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