GetFEM  5.4.2
getfem_mesher.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2004-2020 Julien Pommier, Yves Renard
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**@file getfem_mesher.h
33  @author Julien Pommier <Julien.Pommier@insa-toulouse.fr>, Yves Renard <Yves.Renard@insa-lyon.fr>
34  @date May 1, 2004.
35  @brief An experimental mesher.
36 */
37 
38 #ifndef GETFEM_MESHER_H__
39 #define GETFEM_MESHER_H__
40 
41 
42 #include "getfem_mesh.h"
43 #include "gmm/gmm_condition_number.h"
44 #include "bgeot_comma_init.h"
45 #include "gmm/gmm_solver_bfgs.h"
46 #include "getfem_export.h"
47 #include "bgeot_kdtree.h"
48 #include <typeinfo>
49 
50 namespace getfem {
51 
52  class mesher_virtual_function : virtual public dal::static_stored_object {
53  public:
54  virtual scalar_type operator()(const base_node &P) const = 0;
55  virtual ~mesher_virtual_function() {}
56  };
57 
58  class mvf_constant : public mesher_virtual_function {
59  scalar_type c;
60  public:
61  mvf_constant(scalar_type c_) : c(c_) {}
62  scalar_type operator()(const base_node &) const { return c; }
63  };
64 
65  // Signed distance definition. Should not be a real signed distance but
66  // should satisfy that dist(P) is less than the euclidean distance from
67  // P to the boundary and more than 1/sqrt(N) time this distance.
68 
69 #define SEPS 1e-8
70 
71  class mesher_signed_distance : public mesher_virtual_function {
72  protected:
73  mutable size_type id;
74  public:
75  mesher_signed_distance() : id(size_type(-1)) {}
76  virtual ~mesher_signed_distance() {}
77  virtual bool bounding_box(base_node &bmin, base_node &bmax) const = 0;
78  virtual scalar_type operator()(const base_node &P,
79  dal::bit_vector &bv) const = 0;
80  virtual scalar_type grad(const base_node &P,
81  base_small_vector &G) const = 0;
82  virtual void hess(const base_node &P, base_matrix &H) const = 0;
83  virtual void register_constraints(std::vector<const
84  mesher_signed_distance*>& list) const=0;
85  virtual scalar_type operator()(const base_node &P) const = 0;
86  };
87 
88  typedef std::shared_ptr<const mesher_signed_distance> pmesher_signed_distance;
89 
90  class mesher_half_space : public mesher_signed_distance {
91  base_node x0; base_small_vector n; scalar_type xon;
92  public:
93  mesher_half_space() = default;
94  mesher_half_space(const base_node &x0_, const base_small_vector &n_)
95  : x0(x0_), n(n_)
96  { n /= gmm::vect_norm2(n); xon = gmm::vect_sp(x0, n); }
97  bool bounding_box(base_node &, base_node &) const
98  { return false; }
99  virtual scalar_type operator()(const base_node &P) const
100  { return xon - gmm::vect_sp(P,n); }
101  virtual scalar_type operator()(const base_node &P,
102  dal::bit_vector &bv) const {
103  scalar_type d = xon - gmm::vect_sp(P,n);
104  bv[id] = (gmm::abs(d) < SEPS);
105  return d;
106  }
107  virtual void register_constraints(std::vector<const
108  mesher_signed_distance*>& list) const {
109  id = list.size(); list.push_back(this);
110  }
111  scalar_type grad(const base_node &P, base_small_vector &G) const {
112  G = n; G *= scalar_type(-1);
113  return xon - gmm::vect_sp(P,n);
114  }
115  void hess(const base_node &P, base_matrix &H) const {
116  gmm::resize(H, P.size(), P.size()); gmm::clear(H);
117  }
118 
119  };
120 
121  inline pmesher_signed_distance new_mesher_half_space
122  (const base_node &x0, const base_small_vector &n)
123  { return std::make_shared<mesher_half_space>(x0, n); }
124 
125  class mesher_level_set : public mesher_signed_distance {
126  bgeot::base_poly base;
127  mutable std::vector<base_poly> gradient;
128  mutable std::vector<base_poly> hessian;
129  const fem<base_poly> *pf;
130  mutable int initialized;
131  scalar_type shift_ls; // for the computation of a gap on a level_set.
132  public:
133  bool is_initialized(void) const { return initialized; }
134  mesher_level_set() : initialized(0) {}
135  template <typename VECT>
136  mesher_level_set(pfem pf_, const VECT &coeff_,
137  scalar_type shift_ls_ = scalar_type(0)) {
138  init_base(pf_, coeff_);
139  set_shift(shift_ls_);
140  }
141  void set_shift(scalar_type shift_ls_) {
142  shift_ls = shift_ls_;
143  if (shift_ls != scalar_type(0)) {
144  base_node P(pf->dim()); base_small_vector G(pf->dim());
145  grad(P, G);
146  shift_ls *= gmm::vect_norm2(G);
147  }
148  }
149  template <typename VECT> void init_base(pfem pf_, const VECT &coeff_);
150  void init_grad(void) const;
151  void init_hess(void) const;
152 
153  bool bounding_box(base_node &, base_node &) const
154  { return false; }
155  virtual scalar_type operator()(const base_node &P) const
156  { return bgeot::to_scalar(base.eval(P.begin())) + shift_ls; }
157  virtual scalar_type operator()(const base_node &P,
158  dal::bit_vector &bv) const
159  { scalar_type d = (*this)(P); bv[id] = (gmm::abs(d) < SEPS); return d; }
160  virtual void register_constraints(std::vector<const
161  mesher_signed_distance*>& list) const {
162  id = list.size(); list.push_back(this);
163  }
164  scalar_type grad(const base_node &P, base_small_vector &G) const;
165  void hess(const base_node &P, base_matrix &H) const;
166  };
167 
168  template <typename VECT>
169  void mesher_level_set::init_base(pfem pf_, const VECT &coeff_) {
170  std::vector<scalar_type> coeff(coeff_.begin(), coeff_.end());
171  GMM_ASSERT1(gmm::vect_norm2(coeff) != 0, "level is zero!");
172  pf = dynamic_cast<const fem<base_poly>*> (pf_.get());
173  GMM_ASSERT1(pf, "A polynomial fem is required for level set (got "
174  << typeid(pf_).name() << ")");
175  base = base_poly(pf->base()[0].dim(), pf->base()[0].degree());
176  for (unsigned i=0; i < pf->nb_base(0); ++i) {
177  base += pf->base()[i] * coeff[i];
178  }
179  initialized = 0;
180  }
181 
182  template <typename VECT>
183  inline pmesher_signed_distance new_mesher_level_set
184  (pfem pf, const VECT &coeff, scalar_type shift_ls = scalar_type(0))
185  { return std::make_shared<mesher_level_set>(pf, coeff, shift_ls); }
186 
187 
188  class mesher_ball : public mesher_signed_distance {
189  base_node x0; scalar_type R;
190  public:
191  mesher_ball(base_node x0_, scalar_type R_) : x0(x0_), R(R_) {}
192  bool bounding_box(base_node &bmin, base_node &bmax) const {
193  bmin = bmax = x0;
194  for (size_type i=0; i < x0.size(); ++i) { bmin[i] -= R; bmax[i] += R; }
195  return true;
196  }
197  virtual scalar_type operator()(const base_node &P,
198  dal::bit_vector &bv) const {
199  scalar_type d = gmm::vect_dist2(P,x0)-R;
200  bv[id] = (gmm::abs(d) < SEPS);
201  return d;
202  }
203  virtual scalar_type operator()(const base_node &P) const
204  { return gmm::vect_dist2(P,x0)-R; }
205  virtual void register_constraints(std::vector<const
206  mesher_signed_distance*>& list) const {
207  id = list.size(); list.push_back(this);
208  }
209  scalar_type grad(const base_node &P, base_small_vector &G) const {
210  G = P; G -= x0;
211  scalar_type e= gmm::vect_norm2(G), d = e - R;
212  while (e == scalar_type(0))
213  { gmm::fill_random(G); e = gmm::vect_norm2(G); }
214  G /= e;
215  return d;
216  }
217  void hess(const base_node &, base_matrix &) const {
218  GMM_ASSERT1(false, "Sorry, to be done");
219  }
220  };
221 
222  inline pmesher_signed_distance new_mesher_ball(base_node x0, scalar_type R)
223  { return std::make_shared<mesher_ball>(x0, R); }
224 
225  class mesher_rectangle : public mesher_signed_distance {
226  // ajouter une rotation rigide et translation ..
227  base_node rmin, rmax;
228  std::vector<mesher_half_space> hfs;
229  public:
230  mesher_rectangle(base_node rmin_, base_node rmax_)
231  : rmin(rmin_), rmax(rmax_), hfs(rmin.size()*2) {
232  base_node n(rmin_.size());
233  for (unsigned k = 0; k < rmin.size(); ++k) {
234  n[k] = 1.0;
235  hfs[k*2] = mesher_half_space(rmin, n);
236  n[k] = -1.0;
237  hfs[k*2+1] = mesher_half_space(rmax, n);
238  n[k] = 0.0;
239  }
240  }
241  bool bounding_box(base_node &bmin, base_node &bmax) const {
242  bmin = rmin; bmax = rmax;
243  return true;
244  }
245  virtual scalar_type operator()(const base_node &P) const {
246  size_type N = rmin.size();
247  scalar_type d = rmin[0] - P[0];
248  for (size_type i=0; i < N; ++i) {
249  d = std::max(d, rmin[i] - P[i]);
250  d = std::max(d, P[i] - rmax[i]);
251  }
252  return d;
253  }
254 
255  virtual scalar_type operator()(const base_node &P, dal::bit_vector &bv)
256  const {
257  scalar_type d = this->operator()(P);
258  if (gmm::abs(d) < SEPS)
259  for (int k = 0; k < 2*rmin.size(); ++k) hfs[k](P, bv);
260  return d;
261  }
262  scalar_type grad(const base_node &P, base_small_vector &G) const {
263  unsigned i = 0; scalar_type di = hfs[i](P);
264  for (int k = 1; k < 2*rmin.size(); ++k) {
265  scalar_type dk = hfs[k](P);
266  if (dk > di) { i = k; di = dk; }
267  }
268  return hfs[i].grad(P, G);
269  }
270  void hess(const base_node &P, base_matrix &H) const {
271  gmm::resize(H, P.size(), P.size()); gmm::clear(H);
272  }
273  virtual void register_constraints(std::vector<const
274  mesher_signed_distance*>& list) const {
275  for (int k = 0; k < 2*rmin.size(); ++k)
276  hfs[k].register_constraints(list);
277  }
278  };
279 
280  inline pmesher_signed_distance new_mesher_rectangle(base_node rmin,
281  base_node rmax)
282  { return std::make_shared<mesher_rectangle>(rmin, rmax); }
283 
284  class mesher_simplex_ref : public mesher_signed_distance {
285  // To be added : rigid motion, dilatation ...
286  std::vector<mesher_half_space> hfs;
287  unsigned N;
288  base_node org;
289  public:
290  mesher_simplex_ref(unsigned N_) : hfs(N_+1), N(N_), org(N_) {
291  base_node no(N);
292  for (unsigned k = 0; k < N; ++k) {
293  no[k] = 1.0;
294  hfs[k] = mesher_half_space(org, no);
295  no[k] = 0.0;
296  }
297  std::fill(org.begin(), org.end(), 1.0/N);
298  no = -org;
299  hfs[N] = mesher_half_space(org, no);
300  }
301  bool bounding_box(base_node &bmin, base_node &bmax) const {
302  bmin.resize(N); bmax.resize(N);
303  std::fill(bmin.begin(), bmin.end(), scalar_type(0));
304  std::fill(bmax.begin(), bmax.end(), scalar_type(1));
305  return true;
306  }
307  virtual scalar_type operator()(const base_node &P) const {
308  scalar_type d = - P[0];
309  for (size_type i=1; i < N; ++i) d = std::max(d, - P[i]);
310  d = std::max(d, gmm::vect_sp(P - org, org) / gmm::vect_norm2(org));
311  return d;
312  }
313 
314  virtual scalar_type operator()(const base_node &P, dal::bit_vector &bv)
315  const {
316  scalar_type d = this->operator()(P);
317  if (gmm::abs(d) < SEPS) for (unsigned k = 0; k < N+1; ++k) hfs[k](P, bv);
318  return d;
319  }
320  scalar_type grad(const base_node &P, base_small_vector &G) const {
321  unsigned i = 0; scalar_type di = hfs[i](P);
322  for (unsigned k = 1; k < N+1; ++k) {
323  scalar_type dk = hfs[k](P);
324  if (dk > di) { i = k; di = dk; }
325  }
326  return hfs[i].grad(P, G);
327  }
328  void hess(const base_node &P, base_matrix &H) const {
329  gmm::resize(H, P.size(), P.size()); gmm::clear(H);
330  }
331  virtual void register_constraints(std::vector<const
332  mesher_signed_distance*>& list) const
333  { for (unsigned k = 0; k < N+1; ++k) hfs[k].register_constraints(list); }
334  };
335 
336  inline pmesher_signed_distance new_mesher_simplex_ref(unsigned N)
337  { return std::make_shared<mesher_simplex_ref>(N); }
338 
339 
340  class mesher_prism_ref : public mesher_signed_distance {
341  // To be added : rigid motion, dilatation ...
342  std::vector<mesher_half_space> hfs;
343  unsigned N;
344  base_node org;
345  public:
346  mesher_prism_ref(unsigned N_) : hfs(N_+2), N(N_) {
347  base_node no(N);
348  org = no;
349  for (unsigned k = 0; k < N; ++k) {
350  no[k] = 1.0;
351  hfs[k] = mesher_half_space(org, no);
352  no[k] = 0.0;
353  }
354  no[N-1] = -1.0;
355  org[N-1] = 1.0;
356  hfs[N] = mesher_half_space(org, no);
357  std::fill(org.begin(), org.end(), 1.0/N);
358  org[N-1] = 0.0;
359  no = -org;
360  hfs[N+1] = mesher_half_space(org, no);
361  }
362  bool bounding_box(base_node &bmin, base_node &bmax) const {
363  bmin.resize(N); bmax.resize(N);
364  std::fill(bmin.begin(), bmin.end(), scalar_type(0));
365  std::fill(bmax.begin(), bmax.end(), scalar_type(1));
366  return true;
367  }
368  virtual scalar_type operator()(const base_node &P) const {
369  scalar_type d = - P[0];
370  for (size_type i=1; i < N; ++i) d = std::max(d, - P[i]);
371  d = std::max(d, P[N-1] - scalar_type(1));
372  d = std::max(d, gmm::vect_sp(P - org, org) / gmm::vect_norm2(org));
373  return d;
374  }
375 
376  virtual scalar_type operator()(const base_node &P, dal::bit_vector &bv)
377  const {
378  scalar_type d = this->operator()(P);
379  if (gmm::abs(d) < SEPS) for (unsigned k = 0; k < N+2; ++k) hfs[k](P, bv);
380  return d;
381  }
382  scalar_type grad(const base_node &P, base_small_vector &G) const {
383  unsigned i = 0; scalar_type di = hfs[i](P);
384  for (unsigned k = 1; k < N+2; ++k) {
385  scalar_type dk = hfs[k](P);
386  if (dk > di) { i = k; di = dk; }
387  }
388  return hfs[i].grad(P, G);
389  }
390  void hess(const base_node &P, base_matrix &H) const {
391  gmm::resize(H, P.size(), P.size()); gmm::clear(H);
392  }
393  virtual void register_constraints(std::vector<const
394  mesher_signed_distance*>& list) const
395  { for (unsigned k = 0; k < N+2; ++k) hfs[k].register_constraints(list); }
396  };
397 
398  inline pmesher_signed_distance new_mesher_prism_ref(unsigned N)
399  { return std::make_shared<mesher_prism_ref>(N); }
400 
401 
402 
403  // Rem : It would be better for the convergence of Newton's methods to
404  // take somthing more smooth than min. for instance tthe square root
405  // of positive parts or negative parts depending on the position of
406  // point (in or out the domain).
407 
408  class mesher_union : public mesher_signed_distance {
409  std::vector<pmesher_signed_distance> dists;
410  mutable std::vector<scalar_type> vd;
411  mutable bool isin;
412  bool with_min;
413  public:
414  mesher_union(const std::vector<pmesher_signed_distance>
415  &dists_) : dists(dists_)
416  { vd.resize(dists.size()); with_min = true; }
417 
418  mesher_union
419  (const pmesher_signed_distance &a,
420  const pmesher_signed_distance &b,
421  const pmesher_signed_distance &c = pmesher_signed_distance(),
422  const pmesher_signed_distance &d = pmesher_signed_distance(),
423  const pmesher_signed_distance &e = pmesher_signed_distance(),
424  const pmesher_signed_distance &f = pmesher_signed_distance(),
425  const pmesher_signed_distance &g = pmesher_signed_distance(),
426  const pmesher_signed_distance &h = pmesher_signed_distance(),
427  const pmesher_signed_distance &i = pmesher_signed_distance(),
428  const pmesher_signed_distance &j = pmesher_signed_distance(),
429  const pmesher_signed_distance &k = pmesher_signed_distance(),
430  const pmesher_signed_distance &l = pmesher_signed_distance(),
431  const pmesher_signed_distance &m = pmesher_signed_distance(),
432  const pmesher_signed_distance &n = pmesher_signed_distance(),
433  const pmesher_signed_distance &o = pmesher_signed_distance(),
434  const pmesher_signed_distance &p = pmesher_signed_distance(),
435  const pmesher_signed_distance &q = pmesher_signed_distance(),
436  const pmesher_signed_distance &r = pmesher_signed_distance(),
437  const pmesher_signed_distance &s = pmesher_signed_distance(),
438  const pmesher_signed_distance &t = pmesher_signed_distance()) {
439  dists.push_back(a); dists.push_back(b);
440  with_min = true;
441  if (c) dists.push_back(c);
442  if (d) dists.push_back(d);
443  if (e) dists.push_back(e);
444  if (f) dists.push_back(f);
445  if (g) dists.push_back(g);
446  if (h) dists.push_back(h);
447  if (i) dists.push_back(i);
448  if (j) dists.push_back(j);
449  if (k) dists.push_back(k);
450  if (l) dists.push_back(l);
451  if (m) dists.push_back(m);
452  if (n) dists.push_back(n);
453  if (o) dists.push_back(o);
454  if (p) dists.push_back(p);
455  if (q) dists.push_back(q);
456  if (r) dists.push_back(r);
457  if (s) dists.push_back(s);
458  if (t) dists.push_back(t);
459  vd.resize(dists.size());
460  }
461 
462  bool bounding_box(base_node &bmin, base_node &bmax) const {
463  base_node bmin2, bmax2;
464  bool b = dists[0]->bounding_box(bmin, bmax);
465  if (!b) return false;
466  for (size_type k = 1; k < dists.size(); ++k) {
467  b = dists[k]->bounding_box(bmin2, bmax2);
468  if (!b) return false;
469  for (unsigned i=0; i < bmin.size(); ++i) {
470  bmin[i] = std::min(bmin[i],bmin2[i]);
471  bmax[i] = std::max(bmax[i],bmax2[i]);
472  }
473  }
474  return true;
475  }
476  virtual scalar_type operator()(const base_node &P) const {
477  scalar_type d, f(0), g(1);
478  if (with_min) {
479  d = (*(dists[0]))(P);
480  for (size_type k = 1; k < dists.size(); ++k)
481  d = std::min(d, (*(dists[k]))(P));
482  }
483  else { // essai rat ...
484  isin = false;
485  for (size_type k = 0; k < dists.size(); ++k) {
486  vd[k] = (*(dists[k]))(P);
487  if (vd[k] <= scalar_type(0)) isin = true;
488  f += gmm::sqr(gmm::neg(vd[k]));
489  g *= gmm::pos(vd[k]);
490  }
491  d = isin ? -gmm::sqrt(f)
492  : pow(g, scalar_type(1) / scalar_type(dists.size()));
493  }
494  return d;
495  }
496  scalar_type operator()(const base_node &P, dal::bit_vector &bv) const {
497  if (with_min) {
498  scalar_type d = vd[0] = (*(dists[0]))(P);
499  bool ok = (d > -SEPS);
500  for (size_type k = 1; k < dists.size(); ++k) {
501  vd[k] = (*(dists[k]))(P); if (vd[k] <= -SEPS) ok = false;
502  d = std::min(d,vd[k]);
503  }
504  for (size_type k = 0; ok && k < dists.size(); ++k) {
505  if (vd[k] < SEPS) (*(dists[k]))(P, bv);
506  }
507  return d;
508  }
509  else { // essai rat ...
510  vd[0] = (*(dists[0]))(P);
511  bool ok = (vd[0] > -SEPS);
512  for (size_type k = 1; k < dists.size(); ++k) {
513  vd[k] = (*(dists[k]))(P); if (vd[k] <= -SEPS) ok = false;
514  }
515  for (size_type k = 0; ok && k < dists.size(); ++k) {
516  if (vd[k] < SEPS) (*(dists[k]))(P, bv);
517  }
518  return operator()(P);
519  }
520  }
521  virtual void register_constraints(std::vector<const
522  mesher_signed_distance*>& list) const {
523  for (size_type k = 0; k < dists.size(); ++k)
524  dists[k]->register_constraints(list);
525  }
526  scalar_type grad(const base_node &P, base_small_vector &G) const {
527  scalar_type d;
528  if (with_min /* || gmm::abs(d) < SEPS*/) {
529  d = (*(dists[0]))(P);
530  size_type i = 0;
531  for (size_type k = 1; k < dists.size(); ++k) {
532  scalar_type d2 = (*(dists[k]))(P);
533  if (d2 < d) { d = d2; i = k; }
534  }
535  return dists[i]->grad(P, G);
536  }
537  else { // essai rat ...
538  d = operator()(P);
539  base_small_vector Gloc;
540  for (size_type k = 0; k < dists.size(); ++k) {
541  dists[k]->grad(P, Gloc);
542  if (isin)
543  Gloc *= -gmm::neg(vd[k]);
544  else
545  Gloc *= pow(d, scalar_type(dists.size())) / vd[k];
546  if (!k) G = Gloc; else G += Gloc;
547  }
548  if (isin)
549  G *= scalar_type(1) / d;
550  else
551  G /= pow(d, scalar_type(dists.size()-1)) * scalar_type(dists.size());
552  return d;
553  }
554  }
555  void hess(const base_node &P, base_matrix &H) const {
556  scalar_type d = (*(dists[0]))(P);
557  if (with_min || gmm::abs(d) < SEPS) {
558  size_type i = 0;
559  for (size_type k = 1; k < dists.size(); ++k) {
560  scalar_type d2 = (*(dists[k]))(P);
561  if (d2 < d) { d = d2; i = k; }
562  }
563  dists[i]->hess(P, H);
564  }
565  else {
566  GMM_ASSERT1(false, "Sorry, to be done");
567  }
568  }
569  };
570 
571  inline pmesher_signed_distance new_mesher_union
572  (const std::vector<pmesher_signed_distance> &dists)
573  { return std::make_shared<mesher_union>(dists); }
574 
575  inline pmesher_signed_distance new_mesher_union
576  (const pmesher_signed_distance &a,
577  const pmesher_signed_distance &b,
578  const pmesher_signed_distance &c = pmesher_signed_distance(),
579  const pmesher_signed_distance &d = pmesher_signed_distance(),
580  const pmesher_signed_distance &e = pmesher_signed_distance(),
581  const pmesher_signed_distance &f = pmesher_signed_distance(),
582  const pmesher_signed_distance &g = pmesher_signed_distance(),
583  const pmesher_signed_distance &h = pmesher_signed_distance(),
584  const pmesher_signed_distance &i = pmesher_signed_distance(),
585  const pmesher_signed_distance &j = pmesher_signed_distance(),
586  const pmesher_signed_distance &k = pmesher_signed_distance(),
587  const pmesher_signed_distance &l = pmesher_signed_distance(),
588  const pmesher_signed_distance &m = pmesher_signed_distance(),
589  const pmesher_signed_distance &n = pmesher_signed_distance(),
590  const pmesher_signed_distance &o = pmesher_signed_distance(),
591  const pmesher_signed_distance &p = pmesher_signed_distance(),
592  const pmesher_signed_distance &q = pmesher_signed_distance(),
593  const pmesher_signed_distance &r = pmesher_signed_distance(),
594  const pmesher_signed_distance &s = pmesher_signed_distance(),
595  const pmesher_signed_distance &t = pmesher_signed_distance()) {
596  return std::make_shared<mesher_union>
597  (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t);
598  }
599 
600 
601 
602  class mesher_intersection : public mesher_signed_distance {
603  std::vector<pmesher_signed_distance> dists;
604  mutable std::vector<scalar_type> vd;
605 
606  // const mesher_signed_distance &a, &b;
607  public:
608 
609  mesher_intersection(const std::vector<pmesher_signed_distance>
610  &dists_) : dists(dists_)
611  { vd.resize(dists.size()); }
612 
613  mesher_intersection
614  (const pmesher_signed_distance &a,
615  const pmesher_signed_distance &b,
616  const pmesher_signed_distance &c = pmesher_signed_distance(),
617  const pmesher_signed_distance &d = pmesher_signed_distance(),
618  const pmesher_signed_distance &e = pmesher_signed_distance(),
619  const pmesher_signed_distance &f = pmesher_signed_distance(),
620  const pmesher_signed_distance &g = pmesher_signed_distance(),
621  const pmesher_signed_distance &h = pmesher_signed_distance(),
622  const pmesher_signed_distance &i = pmesher_signed_distance(),
623  const pmesher_signed_distance &j = pmesher_signed_distance(),
624  const pmesher_signed_distance &k = pmesher_signed_distance(),
625  const pmesher_signed_distance &l = pmesher_signed_distance(),
626  const pmesher_signed_distance &m = pmesher_signed_distance(),
627  const pmesher_signed_distance &n = pmesher_signed_distance(),
628  const pmesher_signed_distance &o = pmesher_signed_distance(),
629  const pmesher_signed_distance &p = pmesher_signed_distance(),
630  const pmesher_signed_distance &q = pmesher_signed_distance(),
631  const pmesher_signed_distance &r = pmesher_signed_distance(),
632  const pmesher_signed_distance &s = pmesher_signed_distance(),
633  const pmesher_signed_distance &t = pmesher_signed_distance()) {
634  dists.push_back(a); dists.push_back(b);
635  if (c) dists.push_back(c);
636  if (d) dists.push_back(d);
637  if (e) dists.push_back(e);
638  if (f) dists.push_back(f);
639  if (g) dists.push_back(g);
640  if (h) dists.push_back(h);
641  if (i) dists.push_back(i);
642  if (j) dists.push_back(j);
643  if (k) dists.push_back(k);
644  if (l) dists.push_back(l);
645  if (m) dists.push_back(m);
646  if (n) dists.push_back(n);
647  if (o) dists.push_back(o);
648  if (p) dists.push_back(p);
649  if (q) dists.push_back(q);
650  if (r) dists.push_back(r);
651  if (s) dists.push_back(s);
652  if (t) dists.push_back(t);
653  vd.resize(dists.size());
654  }
655  bool bounding_box(base_node &bmin, base_node &bmax) const {
656  base_node bmin2, bmax2;
657  bool first;
658  bool b = dists[0]->bounding_box(bmin, bmax); first = !b;
659  for (size_type k = 1; k < dists.size(); ++k) {
660  bool bb = dists[k]->bounding_box(bmin2, bmax2);
661  for (unsigned i=0; i < bmin.size() && bb && !first; ++i) {
662  bmin[i] = std::max(bmin[i],bmin2[i]);
663  bmax[i] = std::max(std::min(bmax[i],bmax2[i]), bmin[i]);
664  }
665  if (first && bb) { bmin = bmin2; bmax = bmax2; first = false; }
666  b = b || bb;
667  }
668  return b;
669  }
670  virtual scalar_type operator()(const base_node &P) const {
671  scalar_type d = (*(dists[0]))(P);
672  for (size_type k = 1; k < dists.size(); ++k)
673  d = std::max(d, (*(dists[k]))(P));
674  return d;
675 
676  }
677  scalar_type operator()(const base_node &P, dal::bit_vector &bv) const {
678  scalar_type d = vd[0] = (*(dists[0]))(P);
679  bool ok = (d < SEPS);
680  for (size_type k = 1; k < dists.size(); ++k) {
681  vd[k] = (*(dists[k]))(P); if (vd[k] >= SEPS) ok = false;
682  d = std::min(d,vd[k]);
683  }
684  for (size_type k = 0; ok && k < dists.size(); ++k) {
685  if (vd[k] > -SEPS) (*(dists[k]))(P, bv);
686  }
687  return d;
688  }
689  virtual void register_constraints(std::vector<const
690  mesher_signed_distance*>& list) const {
691  for (size_type k = 0; k < dists.size(); ++k) {
692  dists[k]->register_constraints(list);
693  }
694  }
695  scalar_type grad(const base_node &P, base_small_vector &G) const {
696  scalar_type d = (*(dists[0]))(P);
697  size_type i = 0;
698  for (size_type k = 1; k < dists.size(); ++k) {
699  scalar_type d2 = (*(dists[k]))(P);
700  if (d2 > d) { d = d2; i = k; }
701  }
702  return dists[i]->grad(P, G);
703  }
704  void hess(const base_node &P, base_matrix &H) const {
705  scalar_type d = (*(dists[0]))(P);
706  size_type i = 0;
707  for (size_type k = 1; k < dists.size(); ++k) {
708  scalar_type d2 = (*(dists[k]))(P);
709  if (d2 > d) { d = d2; i = k; }
710  }
711  dists[i]->hess(P, H);
712  }
713  };
714 
715  inline pmesher_signed_distance new_mesher_intersection
716  (const std::vector<pmesher_signed_distance> &dists)
717  { return std::make_shared<mesher_intersection>(dists); }
718 
719  inline pmesher_signed_distance new_mesher_intersection
720  (const pmesher_signed_distance &a,
721  const pmesher_signed_distance &b,
722  const pmesher_signed_distance &c = pmesher_signed_distance(),
723  const pmesher_signed_distance &d = pmesher_signed_distance(),
724  const pmesher_signed_distance &e = pmesher_signed_distance(),
725  const pmesher_signed_distance &f = pmesher_signed_distance(),
726  const pmesher_signed_distance &g = pmesher_signed_distance(),
727  const pmesher_signed_distance &h = pmesher_signed_distance(),
728  const pmesher_signed_distance &i = pmesher_signed_distance(),
729  const pmesher_signed_distance &j = pmesher_signed_distance(),
730  const pmesher_signed_distance &k = pmesher_signed_distance(),
731  const pmesher_signed_distance &l = pmesher_signed_distance(),
732  const pmesher_signed_distance &m = pmesher_signed_distance(),
733  const pmesher_signed_distance &n = pmesher_signed_distance(),
734  const pmesher_signed_distance &o = pmesher_signed_distance(),
735  const pmesher_signed_distance &p = pmesher_signed_distance(),
736  const pmesher_signed_distance &q = pmesher_signed_distance(),
737  const pmesher_signed_distance &r = pmesher_signed_distance(),
738  const pmesher_signed_distance &s = pmesher_signed_distance(),
739  const pmesher_signed_distance &t = pmesher_signed_distance()) {
740  return std::make_shared<mesher_intersection>
741  (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t);
742  }
743 
744 
745 
746  class mesher_setminus : public mesher_signed_distance {
747  pmesher_signed_distance a, b;
748  public:
749  mesher_setminus(const pmesher_signed_distance& a_,
750  const pmesher_signed_distance &b_) : a(a_), b(b_) {}
751  bool bounding_box(base_node &bmin, base_node &bmax) const
752  { return a->bounding_box(bmin,bmax); }
753  scalar_type operator()(const base_node &P, dal::bit_vector &bv) const {
754  scalar_type da = (*a)(P), db = -(*b)(P);
755  if (da < SEPS && db < SEPS) {
756  if (da > -SEPS) (*a)(P, bv);
757  if (db > -SEPS) (*b)(P, bv);
758  }
759  return std::max(da, db);
760  }
761  scalar_type operator()(const base_node &P) const
762  { return std::max((*a)(P),-(*b)(P)); }
763  virtual void register_constraints(std::vector<const
764  mesher_signed_distance*>& list) const {
765  a->register_constraints(list); b->register_constraints(list);
766  }
767  scalar_type grad(const base_node &P, base_small_vector &G) const {
768  scalar_type da = (*a)(P), db = -(*b)(P);
769  if (da > db) return a->grad(P, G);
770  else { b->grad(P, G); G *= scalar_type(-1); return db; }
771  }
772  void hess(const base_node &P, base_matrix &H) const {
773  scalar_type da = (*a)(P), db = -(*b)(P);
774  if (da > db) a->hess(P, H);
775  else { b->hess(P, H); gmm::scale(H, scalar_type(-1)); }
776  }
777 
778  };
779 
780  inline pmesher_signed_distance new_mesher_setminus
781  (const pmesher_signed_distance &a, const pmesher_signed_distance &b)
782  { return std::make_shared<mesher_setminus>(a, b); }
783 
784 
785  class mesher_tube : public mesher_signed_distance {
786  base_node x0; base_node n; scalar_type R;
787  public:
788  mesher_tube(base_node x0_, base_node n_, scalar_type R_)
789  : x0(x0_), n(n_), R(R_)
790  { n /= gmm::vect_norm2(n); }
791  bool bounding_box(base_node &, base_node &) const
792  { return false; }
793  virtual scalar_type operator()(const base_node &P) const {
794  base_node v(P); v -= x0;
795  gmm::add(gmm::scaled(n, -gmm::vect_sp(v, n)), v);
796  return gmm::vect_norm2(v) - R;
797  }
798  virtual scalar_type operator()(const base_node &P,
799  dal::bit_vector &bv) const {
800  scalar_type d = (*this)(P);
801  bv[id] = (gmm::abs(d) < SEPS);
802  return d;
803  }
804  virtual void register_constraints(std::vector<const
805  mesher_signed_distance*>& list) const {
806  id = list.size(); list.push_back(this);
807  }
808  scalar_type grad(const base_node &P, base_small_vector &G) const {
809  G = P; G -= x0;
810  gmm::add(gmm::scaled(n, -gmm::vect_sp(G, n)), G);
811  scalar_type e = gmm::vect_norm2(G), d = e - R;
812  while (e == scalar_type(0)) {
813  gmm::fill_random(G);
814  gmm::add(gmm::scaled(n, -gmm::vect_sp(G, n)), G);
815  e = gmm::vect_norm2(G);
816  }
817  G /= e;
818  return d;
819  }
820  void hess(const base_node &, base_matrix &) const {
821  GMM_ASSERT1(false, "Sorry, to be done");
822  }
823  };
824 
825  inline pmesher_signed_distance new_mesher_tube(base_node x0, base_node n,
826  scalar_type R)
827  { return std::make_shared<mesher_tube>(x0,n,R); }
828 
829 
830  class mesher_cylinder : public mesher_signed_distance {
831  base_node x0; base_small_vector n;
832  scalar_type L, R;
833  pmesher_signed_distance t, p1, p2, i1;
834  public:
835  mesher_cylinder(const base_node &c, const base_small_vector &no,
836  scalar_type L_, scalar_type R_)
837  : x0(c), n(no/gmm::vect_norm2(no)), L(L_), R(R_),
838  t(new_mesher_tube(x0, n, R)),
839  p1(new_mesher_half_space(x0, n)),
840  p2(new_mesher_half_space(x0+n*L, -1.0 * n)),
841  i1(new_mesher_intersection(p1, p2, t)) {}
842  bool bounding_box(base_node &bmin, base_node &bmax) const {
843  base_node x1(x0+n*L);
844  bmin = bmax = x0;
845  for (unsigned i = 0; i < gmm::vect_size(x0); ++i) {
846  bmin[i] = std::min(x0[i], x1[i]) - R;
847  bmax[i] = std::max(x0[i], x1[i]) + R;
848  }
849  return true;
850  }
851  virtual scalar_type operator()(const base_node &P) const {return (*i1)(P); }
852  virtual scalar_type operator()(const base_node &P,
853  dal::bit_vector& bv) const
854  { return (*i1)(P, bv); }
855  scalar_type grad(const base_node &P, base_small_vector &G) const
856  { return i1->grad(P, G); }
857  void hess(const base_node &, base_matrix &) const {
858  GMM_ASSERT1(false, "Sorry, to be done");
859  }
860  virtual void register_constraints(std::vector<const
861  mesher_signed_distance*>& list) const
862  { i1->register_constraints(list); }
863  };
864 
865  inline pmesher_signed_distance new_mesher_cylinder
866  (const base_node &c, const base_small_vector &no,
867  scalar_type L, scalar_type R)
868  { return std::make_shared<mesher_cylinder>(c,no,L,R); }
869 
870 
871  class mesher_infinite_cone : public mesher_signed_distance {
872  // uses the true distance to a cone.
873  base_node x0; base_node n; scalar_type alpha;
874  public:
875  mesher_infinite_cone(base_node x0_, base_node n_, scalar_type alpha_)
876  : x0(x0_), n(n_), alpha(alpha_)
877  { n /= gmm::vect_norm2(n); }
878  bool bounding_box(base_node &, base_node &) const
879  { return false; }
880  virtual scalar_type operator()(const base_node &P) const {
881  base_node v(P); v -= x0;
882  scalar_type v_n = gmm::vect_sp(v, n);
883  gmm::add(gmm::scaled(n, -v_n), v);
884  return gmm::vect_norm2(v) * cos(alpha) - gmm::abs(v_n) * sin(alpha);
885  }
886  virtual scalar_type operator()(const base_node &P,
887  dal::bit_vector &bv) const {
888  scalar_type d = (*this)(P);
889  bv[id] = (gmm::abs(d) < SEPS);
890  return d;
891  }
892  virtual void register_constraints(std::vector<const
893  mesher_signed_distance*>& list) const {
894  id = list.size(); list.push_back(this);
895  }
896  scalar_type grad(const base_node &P, base_small_vector &v) const {
897  v = P; v -= x0;
898  scalar_type v_n = gmm::vect_sp(v, n);
899  gmm::add(gmm::scaled(n, -v_n), v);
900  scalar_type no = gmm::vect_norm2(v);
901  scalar_type d = no * cos(alpha) - gmm::abs(v_n) * sin(alpha);
902  while (no == scalar_type(0)) {
903  gmm::fill_random(v);
904  gmm::add(gmm::scaled(n, -gmm::vect_sp(v, n)), v);
905  no = gmm::vect_norm2(v);
906  }
907  v *= cos(alpha) / no;
908  v -= (sin(alpha) * gmm::sgn(v_n)) * n;
909  return d;
910  }
911  void hess(const base_node &, base_matrix &) const {
912  GMM_ASSERT1(false, "Sorry, to be done");
913  }
914  };
915 
916  inline pmesher_signed_distance new_mesher_infinite_cone
917  (base_node x0, base_node n, scalar_type alpha)
918  { return std::make_shared<mesher_infinite_cone>(x0,n,alpha); }
919 
920 
921  class mesher_cone : public mesher_signed_distance {
922  base_node x0; base_small_vector n;
923  scalar_type L, alpha;
924  pmesher_signed_distance t, p1, p2, i1;
925  public:
926  mesher_cone(const base_node &c, const base_small_vector &no,
927  scalar_type L_, scalar_type alpha_)
928  : x0(c), n(no/gmm::vect_norm2(no)), L(L_), alpha(alpha_),
929  t(new_mesher_infinite_cone(x0, n, alpha)),
930  p1(new_mesher_half_space(x0, n)),
931  p2(new_mesher_half_space(x0+n*L, -1.0 * n)),
932  i1(new_mesher_intersection(p1, p2, t)) {}
933  bool bounding_box(base_node &bmin, base_node &bmax) const {
934  base_node x1(x0+n*L);
935  scalar_type R = L * sin(alpha);
936  bmin = bmax = x0;
937  for (unsigned i = 0; i < gmm::vect_size(x0); ++i) {
938  bmin[i] = std::min(x0[i], x1[i]) - R;
939  bmax[i] = std::max(x0[i], x1[i]) + R;
940  }
941  return true;
942  }
943  virtual scalar_type operator()(const base_node &P) const {return (*i1)(P); }
944  virtual scalar_type operator()(const base_node &P,
945  dal::bit_vector& bv) const
946  { return (*i1)(P, bv); }
947  scalar_type grad(const base_node &P, base_small_vector &G) const
948  { return i1->grad(P, G); }
949  void hess(const base_node &, base_matrix &) const {
950  GMM_ASSERT1(false, "Sorry, to be done");
951  }
952  virtual void register_constraints(std::vector<const
953  mesher_signed_distance*>& list) const
954  { i1->register_constraints(list); }
955  };
956 
957  inline pmesher_signed_distance new_mesher_cone
958  (const base_node &c, const base_small_vector &no,
959  scalar_type L, scalar_type alpha)
960  { return std::make_shared<mesher_cone>(c,no,L,alpha); }
961 
962 
963  class mesher_ellipse : public mesher_signed_distance {
964  base_node x0; base_small_vector n, t;
965  scalar_type r, R, a;
966  public:
967  mesher_ellipse(const base_node &center, const base_small_vector &no,
968  scalar_type r_, scalar_type R_)
969  : x0(center), n(no/gmm::vect_norm2(no)), r(r_), R(R_) {
970  t[0] = -n[1]; t[1] = n[0];
971  if (R < r) { std::swap(r, R); std::swap(n, t); }
972  a = sqrt(R*R - r*r);
973  }
974  bool bounding_box(base_node &bmin, base_node &bmax) const {
975  bmin = bmax = x0;
976  for (unsigned i = 0; i < 2; ++i) { bmin[i] -= R; bmax[i] += R; }
977  return true;
978  }
979  virtual scalar_type operator()(const base_node &P) const {
980  base_small_vector v(P); v -= x0;
981  scalar_type vt = gmm::vect_sp(v, t);
982  vt = std::max(-a, std::min(a, vt));
983  base_node x1 = x0 + t*vt;
984  base_small_vector v1(P); v1 -= x1;
985  scalar_type v1n = gmm::vect_sp(v1, n), v1t = gmm::vect_sp(v1, t);
986  scalar_type x1n = gmm::vect_sp(x1, n), x1t = gmm::vect_sp(x1, t);
987  scalar_type ea = v1n*v1n / (r*r) + v1t * v1t / (R*R);
988  scalar_type eb = 2. * (x1n*v1n / (r*r) + x1t*v1t / (R*R));
989  scalar_type ec = x1n*x1n / (r*r) + x1t * x1t / (R*R);
990 
991  scalar_type delta = eb*eb - 4 * ea * ec;
992  assert(delta >= 0);
993  scalar_type lambda = (-eb + sqrt(delta)) / (2. * ea);
994  return (1.-lambda)*gmm::vect_norm2(v1);
995  }
996  virtual scalar_type operator()(const base_node &P,
997  dal::bit_vector& bv) const {
998  scalar_type d = this->operator()(P);
999  bv[id] = (gmm::abs(d) < SEPS);
1000  return d;
1001  }
1002  scalar_type grad(const base_node &, base_small_vector &) const
1003  { GMM_ASSERT1(false, "Sorry, to be done"); return 0.; }
1004  void hess(const base_node &, base_matrix &) const {
1005  GMM_ASSERT1(false, "Sorry, to be done");
1006  }
1007  virtual void register_constraints(std::vector<const
1008  mesher_signed_distance*>& list) const
1009  { id = list.size(); list.push_back(this); }
1010  };
1011 
1012  inline pmesher_signed_distance new_mesher_ellipse
1013  (const base_node &center, const base_small_vector &no,
1014  scalar_type r, scalar_type R)
1015  { return std::make_shared<mesher_ellipse>(center,no,r,R); }
1016 
1017 
1018 
1019  class mesher_torus : public mesher_signed_distance { // to be done
1020  // ajouter une rotation rigide et translation ..
1021  scalar_type R, r;
1022  public:
1023  mesher_torus(scalar_type RR, scalar_type rr) : R(RR), r(rr) {}
1024  bool bounding_box(base_node &bmin, base_node &bmax) const {
1025  bmin = base_node(3); bmax = base_node(3);
1026  bmin[0] = bmin[1] = -R-r; bmin[2] = -r;
1027  bmax[0] = bmax[1] = +R+r; bmax[2] = +r;
1028  return true;
1029  }
1030  virtual scalar_type operator()(const base_node &P) const {
1031  scalar_type x = P[0], y = P[1], z = P[2], c = sqrt(x*x + y*y);
1032  return (c == 0.) ? R - r : sqrt(gmm::sqr(c-R) + z*z) - r;
1033  }
1034  virtual scalar_type operator()(const base_node &P, dal::bit_vector&bv)
1035  const {
1036  scalar_type d = this->operator()(P);
1037  bv[id] = (gmm::abs(d) < SEPS);
1038  return d;
1039  }
1040  scalar_type grad(const base_node &P, base_small_vector &G) const {
1041  G.resize(3);
1042  scalar_type x = P[0], y = P[1], z = P[2], c = sqrt(x*x + y*y), d(0);
1043  if (c == 0.) {
1044  d = R - r;
1045  gmm::fill_random(G); G[2] = 0.0; G /= gmm::vect_norm2(G);
1046  }
1047  else {
1048  scalar_type w = 1. - R / c, e = sqrt(gmm::sqr(c-R) + z*z);
1049  d = e - r;
1050  if (e == 0.) {
1051  gmm::fill_random(G); G[0] = x; G[1] = y; G /= gmm::vect_norm2(G);
1052  }
1053  else {
1054  G[0] = x * w / e; G[1] = y * w / e; G[2] = z / e;
1055  }
1056  }
1057  return d;
1058  }
1059  void hess(const base_node &, base_matrix &) const {
1060  GMM_ASSERT1(false, "Sorry, to be done");
1061  }
1062  virtual void register_constraints(std::vector<const
1063  mesher_signed_distance*>&list) const
1064  { id = list.size(); list.push_back(this); }
1065  };
1066 
1067 
1068  inline pmesher_signed_distance new_mesher_torus(scalar_type R, scalar_type r)
1069  { return std::make_shared<mesher_torus>(R,r); }
1070 
1071  // mesher
1072  void build_mesh(mesh &m, const pmesher_signed_distance& dist_,
1073  scalar_type h0, const std::vector<base_node> &fixed_points
1074  = std::vector<base_node>(), size_type K = 1, int noise = -1,
1075  size_type iter_max = 500, int prefind = 1,
1076  scalar_type dist_point_hull = 4,
1077  scalar_type boundary_threshold_flatness = 0.11);
1078 
1079  // exported functions
1080  bool try_projection(const mesher_signed_distance& dist, base_node &X,
1081  bool on_surface = false);
1082  bool pure_multi_constraint_projection
1083  (const std::vector<const mesher_signed_distance*> &list_constraints,
1084  base_node &X, const dal::bit_vector &cts);
1085  scalar_type curvature_radius_estimate(const mesher_signed_distance &dist,
1086  base_node X, bool proj = false);
1087  scalar_type min_curvature_radius_estimate
1088  (const std::vector<const mesher_signed_distance*> &list_constraints,
1089  const base_node &X, const dal::bit_vector &cts, size_type hide_first = 0);
1090 
1091 }
1092 
1093 #endif /* GETFEM_MESHER_H__ */
getfem_export.h
Export solutions to various formats.
bgeot_kdtree.h
Simple implementation of a KD-tree.
gmm_solver_bfgs.h
Implements BFGS (Broyden, Fletcher, Goldfarb, Shanno) algorithm.
bgeot::polynomial::eval
T eval(const ITER &it) const
Evaluate the polynomial.
Definition: bgeot_poly.h:453
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::vect_dist2
number_traits< typename linalg_traits< V1 >::value_type >::magnitude_type vect_dist2(const V1 &v1, const V2 &v2)
Euclidean distance between two vectors.
Definition: gmm_blas.h:597
getfem_mesh.h
Define a getfem::getfem_mesh object.
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
bgeot::polynomial< opt_long_scalar_type >
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:557
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
gmm::fill_random
void fill_random(L &l, double cfill)
*‍/
Definition: gmm_blas.h:154
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:47
dal::static_stored_object
base class for static stored objects
Definition: dal_static_stored_objects.h:206
gmm_condition_number.h
computation of the condition number of dense matrices.
bgeot_comma_init.h
convenient initialization of vectors via overload of "operator,".