GetFEM  5.4.2
getfem_mesh_level_set.cc
1 /*===========================================================================
2 
3  Copyright (C) 2005-2020 Julien Pommier
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
23 
24 
25 namespace getfem {
26 
27  std::ostream &operator<<(std::ostream &os, const mesh_level_set::zone &z) {
28  os << "zone[";
29  for (mesh_level_set::zone::const_iterator it = z.begin();
30  it != z.end(); ++it) {
31  if (it != z.begin()) os << ", ";
32  os << **it;
33  }
34  os << "]";
35  return os;
36  }
37 
38  std::ostream &operator<<(std::ostream &os,const mesh_level_set::zoneset &zs) {
39  os << "zoneset[";
40  for (mesh_level_set::zoneset::const_iterator it = zs.begin();
41  it != zs.end(); ++it) {
42  if (it != zs.begin()) os << "; ";
43  os << **it;
44  }
45  os << "]";
46  return os;
47  }
48 
49 
50 #ifdef DEBUG_LS
51 #include <asm/msr.h>
52 struct Chrono {
53  unsigned long long tprev;
54  unsigned long long acc;
55  unsigned long long tmax;
56  bool running; unsigned cnt;
57  Chrono() : acc(0), tmax(0), running(false), cnt(0) {}
58  void tic() { rdtscll(tprev); running = true; ++cnt; }
59  double toc() {
60  if (!running) return 0.; running = false;
61  unsigned long long t; rdtscll(t);
62  t -= tprev;
63  acc += t; tmax = std::max(tmax, t);
64  return double(t)/2.8e9;
65  }
66  double total() { return double(acc) / 2.8e9; }
67  double max() { return double(tmax) / 2.8e9; }
68  double mean() { return cnt ? total() / cnt : 0.; }
69  unsigned count() { return cnt; }
70 };
71 
72  Chrono interpolate_face_chrono;
73 #endif
74 
75  static bool noisy = false;
76  void getfem_mesh_level_set_noisy(void) { noisy = true; }
77 
78  void mesh_level_set::clear(void) {
79  cut_cv.clear();
80  is_adapted_ = false; touch();
81  }
82 
83  const dal::bit_vector &mesh_level_set::crack_tip_convexes() const {
84  return crack_tip_convexes_;
85  }
86 
87  void mesh_level_set::init_with_mesh(mesh &me) {
88  GMM_ASSERT1(linked_mesh_ == 0, "mesh_level_set already initialized");
89  linked_mesh_ = &me;
90  this->add_dependency(me);
91  is_adapted_ = false;
92  }
93 
94  mesh_level_set::mesh_level_set(mesh &me)
95  { linked_mesh_ = 0; init_with_mesh(me); }
96 
97  mesh_level_set::mesh_level_set(void)
98  { linked_mesh_ = 0; is_adapted_ = false; }
99 
100 
101  mesh_level_set::~mesh_level_set() {}
102 
103  static void interpolate_face(const bgeot::pgeometric_trans &pgt,
104  mesh &m, dal::bit_vector& ptdone,
105  const std::vector<size_type>& ipts,
106  const bgeot::pconvex_structure &cvs,
107  size_type nb_vertices,
108  const std::vector<dal::bit_vector> &constraints,
109  const std::vector<const
110  mesher_signed_distance *> &list_constraints) {
111  if (cvs->dim() == 0) return;
112  else if (cvs->dim() > 1) {
113  std::vector<size_type> fpts;
114  for (short_type f=0; f < cvs->nb_faces(); ++f) {
115  fpts.resize(cvs->nb_points_of_face(f));
116  for (size_type k=0; k < fpts.size(); ++k)
117  fpts[k] = ipts[cvs->ind_points_of_face(f)[k]];
118  interpolate_face(pgt, m,ptdone,fpts,cvs->faces_structure()[f],
119  nb_vertices, constraints, list_constraints);
120  }
121  }
122 
123  if (noisy) {
124  cout << "ipts.size() = " << ipts.size() << endl;
125  cout << " nb_vertices = " << nb_vertices << endl;
126  }
127 
128  dal::bit_vector cts; size_type cnt = 0;
129  for (size_type i=0; i < ipts.size(); ++i) {
130  // cout << "ipts[i] = " << ipts[i] << endl;
131  if (ipts[i] < nb_vertices) {
132  if (noisy)
133  cout << "point " << i << " coordinates "
134  << m.points()[ipts[i]] << " constraints[ipts[i]] = "
135  << constraints[ipts[i]] << endl;
136  if (cnt == 0) cts = constraints[ipts[i]];
137  else cts &= constraints[ipts[i]];
138  ++cnt;
139  }
140  }
141 
142  if (noisy) cout << "cts = " << cts << endl;
143 
144  if (cts.card()) {
145  // dal::bit_vector new_cts;
146  for (size_type i=0; i < ipts.size(); ++i) {
147  if (ipts[i] >= nb_vertices && !ptdone[ipts[i]]) {
148  base_node P = m.points()[ipts[i]];
149  // if (cts.card() > 1)
150  // cout << "WARNING, projection sur " << cts << endl;
151  if (!pure_multi_constraint_projection(list_constraints, P, cts)) {
152  GMM_WARNING1("Pure multi has failed in interpolate_face !! "
153  "Original point " << m.points()[ipts[i]]
154  << " projection " << P);
155  } else {
156  if (pgt->convex_ref()->is_in(P) > 1E-8) {
157  GMM_WARNING1("Projected point outside the reference convex ! "
158  "Projection canceled. P = " << P);
159  } else m.points()[ipts[i]] = P;
160  }
161  ptdone[ipts[i]] = true;
162  // dist(P, new_cts);
163  }
164  }
165  }
166  }
167 
168 
169  struct point_stock {
170 
171  bgeot::node_tab points;
172  std::vector<dal::bit_vector> constraints_;
173  std::vector<scalar_type> radius_;
174  const std::vector<const mesher_signed_distance*> &list_constraints;
175  scalar_type radius_cv;
176 
177  void clear(void) { points.clear(); constraints_.clear();radius_.clear(); }
178  scalar_type radius(size_type i) const { return radius_[i]; }
179  const dal::bit_vector &constraints(size_type i) const
180  { return constraints_[i]; }
181  size_type size(void) const { return points.size(); }
182  const base_node &operator[](size_type i) const { return points[i]; }
183 
184  size_type add(const base_node &pt, const dal::bit_vector &bv,
185  scalar_type r) {
186  size_type j = points.search_node(pt);
187  if (j == size_type(-1)) {
188  j = points.add_node(pt);
189  constraints_.push_back(bv);
190  radius_.push_back(r);
191  }
192  return j;
193  }
194  size_type add(const base_node &pt, scalar_type r) {
195  size_type j = points.search_node(pt);
196  if (j == size_type(-1)) {
197  dal::bit_vector bv;
198  for (size_type i = 0; i < list_constraints.size(); ++i)
199  if (gmm::abs((*(list_constraints[i]))(pt)) < 1E-8*radius_cv)
200  bv.add(i);
201  j = points.add_node(pt);
202  constraints_.push_back(bv);
203  radius_.push_back(r);
204  }
205  return j;
206  }
207  size_type add(const base_node &pt) {
208  size_type j = points.search_node(pt);
209  if (j == size_type(-1)) {
210  dal::bit_vector bv;
211  for (size_type i = 0; i < list_constraints.size(); ++i)
212  if (gmm::abs((*(list_constraints[i]))(pt)) < 1E-8*radius_cv)
213  bv.add(i);
214  j = points.add_node(pt);
215  constraints_.push_back(bv);
216  scalar_type r = min_curvature_radius_estimate(list_constraints,pt,bv);
217  radius_.push_back(r);
218  }
219  return j;
220  }
221 
222  dal::bit_vector decimate(const mesher_signed_distance &ref_element,
223  scalar_type dmin) const {
224  dal::bit_vector retained_points;
225  if (size() != 0) {
226  size_type n = points[0].size();
227 
228  bgeot::kdtree points_tree;
229  points_tree.reserve(size());
230  for (size_type i = 0; i < size(); ++i)
231  points_tree.add_point_with_id(points[i], i);
232 
233  for (size_type nb_co = n; nb_co != size_type(-1); nb_co --) {
234  for (size_type i = 0; i < size(); ++i) {
235  if (!(retained_points.is_in(i)) &&
236  (constraints(i).card() == nb_co ||
237  (nb_co == n && constraints(i).card() > n)) &&
238  ref_element(points[i]) < 1E-8) {
239  bool kept = true;
240  scalar_type d = (nb_co == 0) ? (dmin * 1.5)
241  : std::min(radius(i)*0.25, dmin);
242  base_node min = points[i], max = min;
243  for (size_type m = 0; m < n; ++m) { min[m]-=d; max[m]+=d; }
245  points_tree.points_in_box(inpts, min, max);
246  for (size_type j = 0; j < inpts.size(); ++j)
247  if (retained_points.is_in(inpts[j].i) &&
248  constraints(inpts[j].i).contains(constraints(i))
249  && gmm::vect_dist2(points[i], points[inpts[j].i]) < d)
250  { kept = false; break; }
251  if (kept) {
252  if (noisy) cout << "kept point : " << points[i] << " co = "
253  << constraints(i) << endl;
254  retained_points.add(i);
255  }
256  }
257  }
258  }
259  }
260  return retained_points;
261  }
262 
263  point_stock(const std::vector<const mesher_signed_distance*> &ls,
264  scalar_type rcv)
265  : points(scalar_type(10000000)), list_constraints(ls),
266  radius_cv(rcv) {}
267  };
268 
269 
270  static pmesher_signed_distance new_ref_element(bgeot::pgeometric_trans pgt) {
271  dim_type n = pgt->structure()->dim();
272  size_type nbp = pgt->basic_structure()->nb_points();
273  /* Identifying simplexes. */
274  if (nbp == size_type(n+1) &&
275  pgt->basic_structure() == bgeot::simplex_structure(n)) {
276  return new_mesher_simplex_ref(n);
277  }
278 
279  /* Identifying parallelepiped. */
280  if (nbp == (size_type(1) << n) &&
281  pgt->basic_structure() == bgeot::parallelepiped_structure(n)) {
282  base_node rmin(n), rmax(n);
283  std::fill(rmax.begin(), rmax.end(), scalar_type(1));
284  return new_mesher_rectangle(rmin, rmax);
285  }
286 
287  /* Identifying prisms. */
288  if (nbp == size_type(2 * n) &&
289  pgt->basic_structure() == bgeot::prism_P1_structure(n)) {
290  return new_mesher_prism_ref(n);
291  }
292 
293  GMM_ASSERT1(false, "This element is not taken into account. Contact us");
294  }
295 
296 
297  struct global_mesh_for_mesh_level_set : public mesh {};
298  static mesh& global_mesh() {
300  }
301 
302  void mesh_level_set::run_delaunay(std::vector<base_node> &fixed_points,
303  gmm::dense_matrix<size_type> &simplexes,
304  std::vector<dal::bit_vector> &
305  /* fixed_points_constraints */) {
306  double t0=gmm::uclock_sec();
307  if (noisy) cout << "running delaunay with " << fixed_points.size()
308  << " points.." << std::flush;
309  bgeot::qhull_delaunay(fixed_points, simplexes);
310  if (noisy) cout << " -> " << gmm::mat_ncols(simplexes)
311  << " simplexes [" << gmm::uclock_sec()-t0 << "sec]\n";
312  }
313 
314  static bool intersects(const mesh_level_set::zone &z1,
315  const mesh_level_set::zone &z2) {
316  for (std::set<const std::string *>::const_iterator it = z1.begin(); it != z1.end();
317  ++it)
318  if (z2.find(*it) != z2.end()) return true;
319  return false;
320  }
321 
322  void mesh_level_set::merge_zoneset(zoneset &zones1,
323  const zoneset &zones2) const {
324  for (std::set<const zone *>::const_iterator it2 = zones2.begin();
325  it2 != zones2.end(); ++it2) {
326  zone z = *(*it2);
327  for (std::set<const zone *>::iterator it1 = zones1.begin();
328  it1 != zones1.end(); ) {
329  if (intersects(z, *(*it1))) {
330  z.insert((*it1)->begin(), (*it1)->end());
331  zones1.erase(it1++);
332  }
333  else ++it1;
334  }
335  zones1.insert(&(*(allzones.insert(z).first)));
336  }
337  }
338 
339  /* recursively replace '0' by '+' and '-', and the add the new zones */
340  static void add_sub_zones_no_zero(std::string &s,
341  mesh_level_set::zone &z,
342  std::set<std::string> &allsubzones) {
343  size_t i = s.find('0');
344  if (i != size_t(-1)) {
345  s[i] = '+'; add_sub_zones_no_zero(s, z, allsubzones);
346  s[i] = '-'; add_sub_zones_no_zero(s, z, allsubzones);
347  } else {
348  z.insert(&(*(allsubzones.insert(s).first)));
349  }
350  }
351 
352  void mesh_level_set::merge_zoneset(zoneset &zones1,
353  const std::string &subz) const {
354  // very sub-optimal
355  zone z; std::string s(subz);
356  add_sub_zones_no_zero(s, z, allsubzones);
357  zoneset zs;
358  zs.insert(&(*(allzones.insert(z).first)));
359  merge_zoneset(zones1, zs);
360  }
361 
362  /* prezone was filled for the whole convex by find_crossing_level_set.
363  This information is now refined for each sub-convex.
364  */
365  void mesh_level_set::find_zones_of_element(size_type cv,
366  std::string &prezone,
367  scalar_type radius) {
368  convex_info &cvi = cut_cv[cv];
369  cvi.zones.clear();
370  for (dal::bv_visitor i(cvi.pmsh->convex_index()); !i.finished();++i) {
371  // If the sub element is too small, the zone is not taken into account
372  if (cvi.pmsh->convex_area_estimate(i) > 1e-8) {
373  std::string subz = prezone;
374  //cout << "prezone for convex " << cv << " : " << subz << endl;
375  for (size_type j = 0; j < level_sets.size(); ++j) {
376  if (subz[j] == '*' || subz[j] == '0') {
377  int s = sub_simplex_is_not_crossed_by(cv, level_sets[j], i,radius);
378  // cout << "sub_simplex_is_not_crossed_by = " << s << endl;
379  subz[j] = (s < 0) ? '-' : ((s > 0) ? '+' : '0');
380  }
381  }
382  merge_zoneset(cvi.zones, subz);
383  }
384  }
385  if (noisy) cout << "Number of zones for convex " << cv << " : "
386  << cvi.zones.size() << endl;
387  }
388 
389 
390  void mesh_level_set::cut_element(size_type cv,
391  const dal::bit_vector &primary,
392  const dal::bit_vector &secondary,
393  scalar_type radius_cv) {
394 
395  cut_cv[cv] = convex_info();
396  cut_cv[cv].pmsh = std::make_shared<mesh>();
397  if (noisy) cout << "cutting element " << cv << endl;
398  bgeot::pgeometric_trans pgt = linked_mesh().trans_of_convex(cv);
399  pmesher_signed_distance ref_element = new_ref_element(pgt);
400  std::vector<pmesher_signed_distance> mesher_level_sets;
401 
402  size_type n = pgt->structure()->dim();
403  size_type nbtotls = primary.card() + secondary.card();
404  pintegration_method exactint = classical_exact_im(pgt);
405  GMM_ASSERT1(nbtotls <= 16,
406  "An element is cut by more than 16 level_set, aborting");
407 
408  /*
409  * Step 1 : build the signed distances, estimate the curvature radius
410  * and the degree K.
411  */
412  dim_type K = 0; // Max. degree of the level sets.
413  scalar_type r0 = 1E+10; // min curvature radius
414  std::vector<const mesher_signed_distance*> list_constraints;
415  point_stock mesh_points(list_constraints, radius_cv);
416 
417  ref_element->register_constraints(list_constraints);
418  size_type nbeltconstraints = list_constraints.size();
419  mesher_level_sets.reserve(nbtotls);
420  for (size_type ll = 0; ll < level_sets.size(); ++ll) {
421  if (primary[ll]) {
422  base_node X(n); gmm::fill_random(X);
423  K = std::max(K, (level_sets[ll])->degree());
424  mesher_level_sets.push_back(level_sets[ll]->mls_of_convex(cv, 0));
425  pmesher_signed_distance mls(mesher_level_sets.back());
426  list_constraints.push_back(mesher_level_sets.back().get());
427  r0 = std::min(r0, curvature_radius_estimate(*mls, X, true));
428  GMM_ASSERT1(gmm::abs(r0) >= 1e-13, "Something wrong in your level "
429  "set ... curvature radius = " << r0);
430  if (secondary[ll]) {
431  mesher_level_sets.push_back(level_sets[ll]->mls_of_convex(cv, 1));
432  pmesher_signed_distance mls2(mesher_level_sets.back());
433  list_constraints.push_back(mesher_level_sets.back().get());
434  r0 = std::min(r0, curvature_radius_estimate(*mls2, X, true));
435  }
436  }
437  }
438 
439  /*
440  * Step 2 : projection of points of a grid on each signed distance.
441  */
442  scalar_type h0 = std::min(1./(K+1), 0.2 * r0), dmin = 0.55;
443  bool h0_is_ok = true;
444 
445  do {
446 
447  h0_is_ok = true;
448  mesh_points.clear();
449  scalar_type geps = .001*h0;
450  size_type nbpt = 1;
451  std::vector<size_type> gridnx(n);
452  for (size_type i=0; i < n; ++i)
453  { gridnx[i]=1+(size_type)(1./h0); nbpt *= gridnx[i]; }
454  base_node P(n), Q(n), V(n);
455  dal::bit_vector co;
456  std::vector<size_type> co_v;
457 
458  for (size_type i=0; i < nbpt; ++i) {
459  for (size_type k=0, r = i; k < n; ++k) { // building grid point
460  unsigned p = unsigned(r % gridnx[k]);
461  P[k] = p * scalar_type(1) / scalar_type(gridnx[k]-1);
462  r /= gridnx[k];
463  }
464  co.clear(); co_v.resize(0);
465  if ((*ref_element)(P) < geps) {
466  bool kept = false;
467  for (size_type k=list_constraints.size()-1; k!=size_type(-1); --k) {
468  // for (size_type k=0; k < list_constraints.size(); ++k) {
469  // Rerversed to project on level set before on convex boundaries
470  gmm::copy(P, Q);
471  scalar_type d = list_constraints[k]->grad(P, V);
472  if (gmm::vect_norm2(V)*h0*7 > gmm::abs(d))
473  if (try_projection(*(list_constraints[k]), Q, true)) {
474  if (k >= nbeltconstraints
475  && gmm::vect_dist2(P, Q) < h0 * 3.5) kept = true;
476  if (gmm::vect_dist2(P, Q) < h0 / 1.5) {
477  co.add(k); co_v.push_back(k);
478  if ((*ref_element)(Q) < 1E-8) {
479  scalar_type r =
480  curvature_radius_estimate(*(list_constraints[k]), Q);
481  r0 = std::min(r0, r);
482  if (r0 < 4.*h0) { h0 = 0.2 * r0; h0_is_ok = false; break; }
483  if (k >= nbeltconstraints || kept) mesh_points.add(Q, r);
484  }
485  }
486  }
487  }
488  if (kept) mesh_points.add(P, 1.E10);
489  }
490  size_type nb_co = co.card();
491  dal::bit_vector nn;
492  if (h0_is_ok)
493  for (size_type count = 0; count < (size_type(1) << nb_co); ++count) {
494  nn.clear();
495  for (size_type j = 0; j < nb_co; ++j)
496  if (count & (size_type(1) << j)) nn.add(co_v[j]);
497  if (nn.card() > 1 && nn.card() <= n) {
498  if (noisy) cout << "trying set of constraints" << nn << endl;
499  gmm::copy(P, Q);
500  bool ok=pure_multi_constraint_projection(list_constraints,Q,nn);
501  if (ok && (*ref_element)(Q) < 1E-9) {
502  if (noisy) cout << "Intersection point found " << Q << " with "
503  << nn << endl;
504  mesh_points.add(Q);
505  }
506  }
507  }
508  }
509 
510  if (!h0_is_ok) continue;
511 
512  /*
513  * Step 3 : Delaunay, test if a simplex cut a level set and adapt
514  * the mesh to the curved level sets.
515  */
516 
517  dmin = h0;
518  if (noisy) cout << "dmin = " << dmin << " h0 = " << h0 << endl;
519  if (noisy) cout << "convex " << cv << endl;
520 
521  dal::bit_vector retained_points
522  = mesh_points.decimate(*ref_element, dmin);
523 
524  bool delaunay_again = true;
525 
526  std::vector<base_node> fixed_points;
527  std::vector<dal::bit_vector> fixed_points_constraints;
528  mesh &msh(*(cut_cv[cv].pmsh));
529 
530  mesh_region &ls_border_faces(cut_cv[cv].ls_border_faces);
531  std::vector<base_node> cvpts;
532 
533  size_type nb_delaunay = 0;
534 
535  while (delaunay_again) {
536  if (++nb_delaunay > 15)
537  { h0_is_ok = false; h0 /= 2.0; dmin = 2.*h0; break; }
538 
539  fixed_points.resize(0);
540  fixed_points_constraints.resize(0);
541  // std::vector<base_node> fixed_points;
542  // std::vector<dal::bit_vector> fixed_points_constraints;
543  fixed_points.reserve(retained_points.card());
544  fixed_points_constraints.reserve(retained_points.card());
545  for (dal::bv_visitor i(retained_points); !i.finished(); ++i) {
546  fixed_points.push_back(mesh_points[i]);
547  fixed_points_constraints.push_back(mesh_points.constraints(i));
548  }
549 
550  gmm::dense_matrix<size_type> simplexes;
551  run_delaunay(fixed_points, simplexes, fixed_points_constraints);
552  delaunay_again = false;
553 
554  msh.clear();
555 
556  for (size_type i = 0; i < fixed_points.size(); ++i) {
557  size_type j = msh.add_point(fixed_points[i]);
558  assert(j == i);
559  }
560 
561  cvpts.resize(0);
562  for (size_type i = 0; i < gmm::mat_ncols(simplexes); ++i) {
563  size_type j = msh.add_convex(bgeot::simplex_geotrans(n,1),
564  gmm::vect_const_begin(gmm::mat_col(simplexes, i)));
565  /* remove flat convexes => beware the final mesh may not be conformal ! */
566  if (msh.convex_quality_estimate(j) < 1E-10) msh.sup_convex(j);
567  else {
568  std::vector<scalar_type> signs(list_constraints.size());
569  std::vector<size_type> prev_point(list_constraints.size());
570  for (size_type ii = 0; ii <= n; ++ii) {
571  for (size_type jj = 0; jj < list_constraints.size(); ++jj) {
572  scalar_type dd =
573  (*(list_constraints[jj]))(msh.points_of_convex(j)[ii]);
574  if (gmm::abs(dd) > radius_cv * 1E-7) {
575  if (dd * signs[jj] < 0.0) {
576  if (noisy) cout << "Intersection found ... " << jj
577  << " dd = " << dd << " convex quality = "
578  << msh.convex_quality_estimate(j) << "\n";
579  // calcul d'intersection
580  base_node X = msh.points_of_convex(j)[ii], G;
581  base_node VV = msh.points_of_convex(j)[prev_point[jj]] - X;
582  if (dd > 0.) gmm::scale(VV, -1.);
583  dd = (*(list_constraints[jj])).grad(X, G);
584  size_type nbit = 0;
585  while (gmm::abs(dd) > 1e-16*radius_cv && (++nbit < 20)) { // Newton
586  scalar_type nG = gmm::vect_sp(G, VV);
587  if (gmm::abs(nG) < 1E-8) nG = 1E-8;
588  if (nG < 0) nG = 1.0;
589  gmm::add(gmm::scaled(VV, -dd / nG), X);
590  dd = (*(list_constraints[jj])).grad(X, G);
591  }
592  if (gmm::abs(dd) > 1e-16*radius_cv) { // brute force dychotomy
593  base_node X1 = msh.points_of_convex(j)[ii];
594  base_node X2 = msh.points_of_convex(j)[prev_point[jj]], X3;
595  scalar_type dd1 = (*(list_constraints[jj]))(X1);
596  scalar_type dd2 = (*(list_constraints[jj]))(X2);
597  if (dd1 > dd2) { std::swap(dd1, dd2); std::swap(X1, X2); }
598  while (gmm::abs(dd1) > 1e-15*radius_cv) {
599  X3 = (X1 + X2) / 2.0;
600  scalar_type dd3 = (*(list_constraints[jj]))(X3);
601  if (dd3 == dd1 || dd3 == dd2) break;
602  if (dd3 > 0) { dd2 = dd3; X2 = X3; }
603  else { dd1 = dd3; X1 = X3; }
604  }
605  X = X1;
606  }
607 
608  size_type kk = mesh_points.add(X);
609  if (!(retained_points[kk])) {
610  retained_points.add(kk);
611  delaunay_again = true;
612  break;
613  // goto delaunay_fin;
614  }
615  }
616  if (signs[jj] == 0.0) { signs[jj] = dd; prev_point[jj] = ii; }
617  }
618  }
619  }
620  }
621  }
622 
623  // delaunay_fin :;
624 
625  }
626  if (!h0_is_ok) continue;
627 
628  // Produces an order K mesh for K > 1 (with affine element for the moment)
629  if (K>1) {
630  for (dal::bv_visitor_c j(msh.convex_index()); !j.finished(); ++j) {
631  bgeot::pgeometric_trans pgt2 = bgeot::simplex_geotrans(n, K);
632  cvpts.resize(pgt2->nb_points());
633  for (size_type k=0; k < pgt2->nb_points(); ++k) {
634  cvpts[k] = bgeot::simplex_geotrans(n,1)->transform
635  (pgt2->convex_ref()->points()[k], msh.points_of_convex(j));
636  }
637  msh.sup_convex(j);
638  /* some new point will be added to the mesh, but the initial ones
639  (those on the convex vertices) are kepts */
640  size_type jj = msh.add_convex_by_points(pgt2, cvpts.begin());
641  assert(jj == j);
642  }
643  }
644 
645  if (noisy) {
647  sl.build(msh, getfem::slicer_none(), 1);
648  char s[512]; sprintf(s, "totobefore%d.dx", int(cv));
649  getfem::dx_export exp(s);
650  exp.exporting(sl);
651  exp.exporting_mesh_edges();
652  exp.write_mesh();
653  }
654 
655  /* detect the faces lying on level_set boundaries
656  (not exact, some other faces maybe be found, use this only for vizualistion) */
657  for (dal::bv_visitor i(msh.convex_index()); !i.finished(); ++i) {
658  for (short_type f = 0; f <= n; ++f) {
659  const mesh::ind_pt_face_ct &fpts
660  = msh.ind_points_of_face_of_convex(i, f);
661 
662  dal::bit_vector cts; bool first = true;
663  for (unsigned k=0; k < fpts.size(); ++k) {
664  if (fpts[k] >= fixed_points_constraints.size()) {
665  assert(K>1);
666  continue; /* ignore additional point inserted when K>1 */
667  }
668  if (first) { cts = fixed_points_constraints[fpts[k]]; first=false; }
669  else cts &= fixed_points_constraints[fpts[k]];
670  }
671  for (dal::bv_visitor ii(cts); !ii.finished(); ++ii) {
672  if (ii >= nbeltconstraints)
673  ls_border_faces.add(i, f);
674  }
675  }
676  }
677 
678 
679  if (K > 1) { // To be done only on the level-set ...
680  dal::bit_vector ptdone;
681  std::vector<size_type> ipts;
682  // for (mr_visitor icv(ls_border_faces); !icv.finished(); ++icv) {
683  // size_type i = icv.cv();
684  // short_type f = icv.f();
685  for (dal::bv_visitor i(msh.convex_index()); !i.finished(); ++i) {
686  for (short_type f = 0; f <= n; ++f) {
687  const mesh::ind_pt_face_ct &fpts
688  = msh.ind_points_of_face_of_convex(i, f);
689  ipts.assign(fpts.begin(), fpts.end());
690 #ifdef DEBUG_LS
691  interpolate_face_chrono.tic();
692 #endif
693 
694  interpolate_face(pgt, msh, ptdone, ipts,
695  msh.trans_of_convex(i)->structure()
696  ->faces_structure()[f], fixed_points.size(),
697  fixed_points_constraints, list_constraints);
698 #ifdef DEBUG_LS
699  interpolate_face_chrono.toc();
700 #endif
701  }
702  }
703  }
704 
705  if (noisy) {
707  sl.build(msh, getfem::slicer_none(), 12);
708  char s[512]; sprintf(s, "toto%d.dx", int(cv));
709  getfem::dx_export exp(s);
710  exp.exporting(sl);
711  exp.exporting_mesh_edges();
712  exp.write_mesh();
713  }
714 
715  /*
716  * Step 4 : Building the new integration method.
717  */
718  base_matrix G;
719  bgeot::pgeometric_trans pgt2 = bgeot::simplex_geotrans(n, K);
720  papprox_integration
721  pai = classical_approx_im(pgt2,dim_type(2*K))->approx_method();
722  auto new_approx = std::make_shared<approx_integration>(pgt->convex_ref());
723  base_matrix KK(n,n), CS(n,n);
724  base_matrix pc(pgt2->nb_points(), n);
725  for (dal::bv_visitor i(msh.convex_index()); !i.finished(); ++i) {
726  vectors_to_base_matrix(G, msh.points_of_convex(i));
727  bgeot::geotrans_interpolation_context c(msh.trans_of_convex(i),
728  pai->point(0), G);
729  scalar_type sign = 0.0;
730  for (size_type j = 0; j < pai->nb_points_on_convex(); ++j) {
731  c.set_xref(pai->point(j));
732  pgt2->poly_vector_grad(pai->point(j), pc);
733  gmm::mult(G,pc,KK);
734  scalar_type J = gmm::lu_det(KK);
735  if (noisy && J * sign < 0) {
736  cout << "CAUTION reversal situation in convex " << i
737  << "sign = " << sign << " J = " << J << endl;
738  for (size_type nip = 0; nip < msh.nb_points_of_convex(i); ++nip)
739  cout << "Point " << nip << "=" << msh.points_of_convex(i)[nip]
740  << " ";
741  cout << endl;
742  }
743  if (sign == 0 && gmm::abs(J) > 1E-14) sign = J;
744  new_approx->add_point(c.xreal(), pai->coeff(j) * gmm::abs(J));
745  }
746  }
747 
748  getfem::mesh_region border_faces;
749  getfem::outer_faces_of_mesh(msh, border_faces);
750  for (getfem::mr_visitor it(border_faces); !it.finished(); ++it) {
751  vectors_to_base_matrix(G, msh.points_of_convex(it.cv()));
752  bgeot::geotrans_interpolation_context c(msh.trans_of_convex(it.cv()),
753  pai->point(0), G);
754  for (size_type j = 0; j < pai->nb_points_on_face(it.f()); ++j) {
755  c.set_xref(pai->point_on_face(it.f(), j));
756  new_approx->add_point(c.xreal(), pai->coeff_on_face(it.f(), j)
757  * gmm::abs(c.J()), it.f() /* faux */);
758  }
759  }
760  new_approx->valid_method();
761 
762  /*
763  * Step 5 : Test the validity of produced integration method.
764  */
765 
766  scalar_type error = test_integration_error(new_approx, 1);
767  if (noisy) cout << " max monomial integration err: " << error << "\n";
768  if (error > 1e-5) {
769  if (noisy) cout << "Not Good ! Let us make a finer cut.\n";
770  if (dmin > 3*h0) { dmin /= 2.; }
771  else { h0 /= 2.0; dmin = 2.*h0; }
772  h0_is_ok = false;
773  }
774 
775 
776  if (h0_is_ok && noisy) { // ajout dans global mesh pour visu
777  vectors_to_base_matrix(G, linked_mesh().points_of_convex(cv));
778  std::vector<size_type> pts(msh.nb_points());
779  for (size_type i = 0; i < msh.nb_points(); ++i)
780  pts[i] = global_mesh().add_point(pgt->transform(msh.points()[i], G));
781 
782  for (dal::bv_visitor i(msh.convex_index()); !i.finished(); ++i)
783  global_mesh().add_convex(msh.trans_of_convex(i),
784  gmm::index_ref_iterator(pts.begin(),
785  msh.ind_points_of_convex(i).begin()));
786  }
787 
788  } while (!h0_is_ok);
789 
790 #ifdef DEBUG_LS
791  cout << "Interpolate face: " << interpolate_face_chrono.total()
792  << " moyenne: " << interpolate_face_chrono.mean() << "\n";
793 #endif
794  }
795 
796 
798  m.clear();
799  base_matrix G;
800  for (dal::bv_visitor cv(linked_mesh().convex_index());
801  !cv.finished(); ++cv) {
802  if (is_convex_cut(cv)) {
803  const convex_info &ci = (cut_cv.find(cv))->second;
804  mesh &msh = *(ci.pmsh);
805  bgeot::pgeometric_trans pgt = linked_mesh().trans_of_convex(cv);
806  vectors_to_base_matrix(G, linked_mesh().points_of_convex(cv));
807  std::vector<size_type> pts(msh.nb_points());
808  for (size_type i = 0; i < msh.nb_points(); ++i)
809  pts[i] = m.add_point(pgt->transform(msh.points()[i], G));
810 
811  std::vector<size_type> ic2(msh.nb_allocated_convex());
812 
813  for (dal::bv_visitor i(msh.convex_index()); !i.finished(); ++i) {
814  ic2[i] = m.add_convex(msh.trans_of_convex(i),
815  gmm::index_ref_iterator(pts.begin(),
816  msh.ind_points_of_convex(i).begin()));
817 
818  }
819  for (mr_visitor i(ci.ls_border_faces); !i.finished(); ++i) {
820  m.region(0).add(ic2[i.cv()], i.f());
821  }
822 
823  } else {
824  m.add_convex_by_points(linked_mesh().trans_of_convex(cv),
825  linked_mesh().points_of_convex(cv).begin());
826  }
827  }
828 
829 
830  }
831 
832  void mesh_level_set::update_crack_tip_convexes() {
833  crack_tip_convexes_.clear();
834 
835  for (std::map<size_type, convex_info>::const_iterator it = cut_cv.begin();
836  it != cut_cv.end(); ++it) {
837  size_type cv = it->first;
838  mesh &msh = *(it->second.pmsh);
839  for (unsigned ils = 0; ils < nb_level_sets(); ++ils) {
840  if (get_level_set(ils)->has_secondary()) {
841  pmesher_signed_distance
842  mesherls0 = get_level_set(ils)->mls_of_convex(cv, 0, false);
843  pmesher_signed_distance
844  mesherls1 = get_level_set(ils)->mls_of_convex(cv, 1, false);
845  for (dal::bv_visitor ii(msh.convex_index()); !ii.finished(); ++ii) {
846  for (unsigned ipt = 0; ipt < msh.nb_points_of_convex(ii); ++ipt) {
847  if (gmm::abs((*mesherls0)(msh.points_of_convex(ii)[ipt])) < 1E-10
848  && gmm::abs((*mesherls1)(msh.points_of_convex(ii)[ipt])) < 1E-10) {
849  crack_tip_convexes_.add(cv);
850  goto next_convex;
851  }
852  }
853  }
854  }
855  }
856  next_convex:
857  ;
858  }
859  }
860 
862 
863  // compute the elements touched by each level set
864  // for each element touched, compute the sub mesh
865  // then compute the adapted integration method
866  GMM_ASSERT1(linked_mesh_ != 0, "Uninitialized mesh_level_set");
867  cut_cv.clear();
868  allsubzones.clear();
869  zones_of_convexes.clear();
870  allzones.clear();
871 
872  // noisy = true;
873 
874  std::string z;
875  for (dal::bv_visitor cv(linked_mesh().convex_index());
876  !cv.finished(); ++cv) {
877  scalar_type radius = linked_mesh().convex_radius_estimate(cv);
878  dal::bit_vector prim, sec;
879  find_crossing_level_set(cv, prim, sec, z, radius);
880  zones_of_convexes[cv] = &(*(allsubzones.insert(z).first));
881  if (noisy) cout << "element " << cv << " cut level sets : "
882  << prim << " zone : " << z << endl;
883  if (prim.card()) {
884  cut_element(cv, prim, sec, radius);
885  find_zones_of_element(cv, z, radius);
886  }
887  }
888  if (noisy) {
890  sl.build(global_mesh(), getfem::slicer_none(), 6);
891  getfem::dx_export exp("totoglob.dx");
892  exp.exporting(sl);
893  exp.exporting_mesh_edges();
894  exp.write_mesh();
895  }
896 
897  update_crack_tip_convexes();
898  is_adapted_ = true;
899  }
900 
901  // detect the intersection of two or more level sets and the convexes
902  // where the intersection occurs. To be alled after adapt.
903  void mesh_level_set::find_level_set_potential_intersections
904  (std::vector<size_type> &icv, std::vector<dal::bit_vector> &ils) {
905 
906  GMM_ASSERT1(linked_mesh_ != 0, "Uninitialized mesh_level_set");
907  std::string z;
908  for (dal::bv_visitor cv(linked_mesh().convex_index());
909  !cv.finished(); ++cv)
910  if (is_convex_cut(cv)) {
911  scalar_type radius = linked_mesh().convex_radius_estimate(cv);
912  dal::bit_vector prim, sec;
913  find_crossing_level_set(cv, prim, sec, z, radius);
914  if (prim.card() > 1) {
915  icv.push_back(cv);
916  ils.push_back(prim);
917  }
918  }
919  }
920 
921  // return +1 if the sub-element is on the one side of the level-set
922  // -1 if the sub-element is on the other side of the level-set
923  // 0 if the sub-element is one the positive part of the secundary
924  // level-set if any.
925  int mesh_level_set::sub_simplex_is_not_crossed_by(size_type cv,
926  plevel_set ls,
927  size_type sub_cv,
928  scalar_type radius) {
929  scalar_type EPS = 1e-7 * radius;
930  bgeot::pgeometric_trans pgt = linked_mesh().trans_of_convex(cv);
931  convex_info &cvi = cut_cv[cv];
932  bgeot::pgeometric_trans pgt2 = cvi.pmsh->trans_of_convex(sub_cv);
933 
934  // cout << "cv " << cv << " radius = " << radius << endl;
935 
936  pmesher_signed_distance mls0 = ls->mls_of_convex(cv, 0), mls1(mls0);
937  if (ls->has_secondary()) mls1 = ls->mls_of_convex(cv, 1);
938  int p = 0;
939  bool is_cut = false;
940  scalar_type d2 = 0, d1 = 1, d0 = 0, d0min = 0;
941  for (size_type i = 0; i < pgt2->nb_points(); ++i) {
942  d0 = (*mls0)(cvi.pmsh->points_of_convex(sub_cv)[i]);
943  if (i == 0) d0min = gmm::abs(d0);
944  else d0min = std::min(d0min, gmm::abs(d0));
945  if (ls->has_secondary())
946  d1 = std::min(d1, (*mls1)(cvi.pmsh->points_of_convex(sub_cv)[i]));
947 
948  int p2 = ( (d0 < -EPS) ? -1 : ((d0 > EPS) ? +1 : 0));
949  if (p == 0) p = p2;
950  if (gmm::abs(d0) > gmm::abs(d2)) d2 = d0;
951  if (!p2 || p*p2 < 0) is_cut = true;
952  }
953  // cout << "d0min = " << d0min << " d1 = " << d1 << " iscut = "
954  // << is_cut << endl;
955  if (is_cut && ls->has_secondary() && d1 >= -radius*0.0001) return 0;
956  // if (d0min < EPS && ls->has_secondary() && d1 >= -EPS) return 0;
957  return (d2 < 0.) ? -1 : 1;
958  }
959 
960  int mesh_level_set::is_not_crossed_by(size_type cv, plevel_set ls,
961  unsigned lsnum, scalar_type radius) {
962  const mesh_fem &mf = ls->get_mesh_fem();
963  GMM_ASSERT1(!mf.is_reduced(), "Internal error");
964  const mesh_fem::ind_dof_ct &dofs = mf.ind_basic_dof_of_element(cv);
965  pfem pf = mf.fem_of_element(cv);
966  int p = -2;
967  scalar_type EPS = 1e-8 * radius;
968 
969  /* easy cases:
970  - different sign on the dof nodes => intersection for sure
971  - min value of the levelset greater than the radius of the convex
972  => no intersection
973  */
974  for (mesh_fem::ind_dof_ct::const_iterator it=dofs.begin();
975  it != dofs.end(); ++it) {
976  scalar_type v = ls->values(lsnum)[*it];
977  int p2 = ( (v < -EPS) ? -1 : ((v > EPS) ? +1 : 0));
978  if (p == -2) p=p2;
979  if (!p2 || p*p2 < 0) return 0;
980  }
981 
982  pmesher_signed_distance mls1 = ls->mls_of_convex(cv, lsnum, false);
983  base_node X(pf->dim()), G(pf->dim());
984  gmm::fill_random(X); X *= 1E-2;
985  scalar_type d = mls1->grad(X, G);
986  if (gmm::vect_norm2(G)*2.5 < gmm::abs(d)) return p;
987 
988  bgeot::pgeometric_trans pgt = linked_mesh().trans_of_convex(cv);
989  pmesher_signed_distance ref_element = new_ref_element(pgt);
990 
991  gmm::fill_random(X); X *= 1E-2;
992  mesher_intersection mi1(ref_element, mls1);
993  if (!try_projection(mi1, X)) return p;
994  if ((*ref_element)(X) > 1E-8) return p;
995 
996  gmm::fill_random(X); X *= 1E-2;
997  pmesher_signed_distance mls2 = ls->mls_of_convex(cv, lsnum, true);
998  mesher_intersection mi2(ref_element, mls2);
999  if (!try_projection(mi2, X)) return p;
1000  if ((*ref_element)(X) > 1E-8) return p;
1001 
1002  return 0;
1003  }
1004 
1005  void mesh_level_set::find_crossing_level_set(size_type cv,
1006  dal::bit_vector &prim,
1007  dal::bit_vector &sec,
1008  std::string &z,
1009  scalar_type radius) {
1010  prim.clear(); sec.clear();
1011  z = std::string(level_sets.size(), '*');
1012  unsigned lsnum = 0;
1013  for (size_type k = 0; k < level_sets.size(); ++k, ++lsnum) {
1014  if (noisy) cout << "testing cv " << cv << " with level set "
1015  << k << endl;
1016  int s = is_not_crossed_by(cv, level_sets[k], 0, radius);
1017  if (!s) {
1018  if (noisy) cout << "is cut \n";
1019  if (level_sets[k]->has_secondary()) {
1020  s = is_not_crossed_by(cv, level_sets[k], 1, radius);
1021  if (!s) { sec.add(lsnum); prim.add(lsnum); }
1022  else if (s < 0) prim.add(lsnum); else z[k] = '0';
1023  }
1024  else prim.add(lsnum);
1025  }
1026  else z[k] = (s < 0) ? '-' : '+';
1027  }
1028  }
1029 
1030 
1031 } /* end of namespace getfem. */
1032 
1033 
1034 
dal::singleton::instance
static T & instance()
Instance from the current thread.
Definition: dal_singleton.h:165
getfem::mesh::clear
void clear()
Erase the mesh.
Definition: getfem_mesh.cc:273
bgeot::operator<<
std::ostream & operator<<(std::ostream &o, const convex_structure &cv)
Print the details of the convex structure cvs to the output stream o.
Definition: bgeot_convex_structure.cc:71
getfem::stored_mesh_slice
The output of a getfem::mesh_slicer which has been recorded.
Definition: getfem_mesh_slice.h:47
getfem::mesh_level_set::global_cut_mesh
void global_cut_mesh(mesh &m) const
fill m with the (non-conformal) "cut" mesh.
Definition: getfem_mesh_level_set.cc:797
bgeot::prism_P1_structure
pconvex_structure prism_P1_structure(dim_type nc)
Give a pointer on the structures of a prism of dimension d.
Definition: bgeot_convex_structure.h:206
getfem::mesh::region
const mesh_region region(size_type id) const
Return the region of index 'id'.
Definition: getfem_mesh.h:421
bgeot::kdtree_tab_type
std::vector< index_node_pair > kdtree_tab_type
store a set of points with associated indexes.
Definition: bgeot_kdtree.h:68
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
getfem::mesh_level_set::adapt
void adapt(void)
do all the work (cut the convexes wrt the levelsets)
Definition: getfem_mesh_level_set.cc:861
getfem_mesh_level_set.h
Keep informations about a mesh crossed by level-sets.
bgeot::geotrans_interpolation_context
the geotrans_interpolation_context structure is passed as the argument of geometric transformation in...
Definition: bgeot_geometric_trans.h:411
getfem::mesh_level_set::nb_level_sets
size_type nb_level_sets(void) const
Get number of level-sets referenced in this object.
Definition: getfem_mesh_level_set.h:86
getfem::mesh::add_convex
size_type add_convex(bgeot::pgeometric_trans pgt, ITER ipts)
Add a convex to the mesh.
Definition: getfem_mesh.h:226
bgeot::mesh_structure::ind_points_of_convex
const ind_set & ind_points_of_convex(size_type ic) const
Return a container to the list of points attached to convex ic.
Definition: bgeot_mesh_structure.h:107
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::simplex_structure
pconvex_structure simplex_structure(dim_type nc)
Give a pointer on the structures of a simplex of dimension d.
Definition: bgeot_convex_structure.cc:148
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
bgeot::mesh_structure::convex_index
const dal::bit_vector & convex_index() const
Return the list of valid convex IDs.
Definition: bgeot_mesh_structure.h:91
bgeot::kdtree
Balanced tree over a set of points.
Definition: bgeot_kdtree.h:102
dal::dynamic_array::clear
void clear(void)
Clear and desallocate all the elements.
Definition: dal_basic.h:298
bgeot::node_tab::search_node
size_type search_node(const base_node &pt, const scalar_type radius=0) const
Search a node in the array.
Definition: bgeot_node_tab.cc:56
getfem::mesh::convex_radius_estimate
virtual scalar_type convex_radius_estimate(size_type ic) const
Return an estimate of the convex largest dimension.
Definition: getfem_mesh.cc:461
getfem::mesh_region::visitor
"iterator" class for regions.
Definition: getfem_mesh_region.h:237
getfem::classical_exact_im
pintegration_method classical_exact_im(bgeot::pgeometric_trans pgt)
return an exact integration method for convex type handled by pgt.
Definition: getfem_integration.cc:1235
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
getfem::slicer_none
This slicer does nothing.
Definition: getfem_mesh_slicers.h:245
bgeot::node_tab::clear
void clear(void)
reset the array, remove all points
Definition: bgeot_node_tab.cc:89
getfem::mesh::add_convex_by_points
size_type add_convex_by_points(bgeot::pgeometric_trans pgt, ITER ipts, const scalar_type tol=scalar_type(0))
Add a convex to the mesh, given a geometric transformation and a list of point coordinates.
Definition: getfem_mesh.h:558
gmm::fill_random
void fill_random(L &l, double cfill)
*‍/
Definition: gmm_blas.h:154
bgeot::node_tab
Store a set of points, identifying points that are nearer than a certain very small distance.
Definition: bgeot_node_tab.h:52
getfem::stored_mesh_slice::build
void build(const getfem::mesh &m, const slicer_action &a, size_type nrefine=1)
Build the slice, by applying a slicer_action operation.
Definition: getfem_mesh_slice.h:197
bgeot::node_tab::add_node
size_type add_node(const base_node &pt, const scalar_type radius=0, bool remove_duplicated_nodes=true)
Add a point to the array or use an existing point, located within a distance smaller than radius.
Definition: bgeot_node_tab.cc:96
getfem::classical_approx_im
pintegration_method classical_approx_im(bgeot::pgeometric_trans pgt, dim_type degree)
try to find an approximate integration method for the geometric transformation pgt which is able to i...
Definition: getfem_integration.cc:1280
getfem::mesh::nb_points
size_type nb_points() const
Give the number of geometrical nodes in the mesh.
Definition: getfem_mesh.h:194
bgeot::pconvex_structure
std::shared_ptr< const convex_structure > pconvex_structure
Pointer on a convex structure description.
Definition: bgeot_convex_structure.h:54
bgeot::mesh_structure::nb_allocated_convex
size_type nb_allocated_convex() const
The number of convex indexes from 0 to the index of the last convex.
Definition: bgeot_mesh_structure.h:98
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
bgeot::parallelepiped_structure
pconvex_structure parallelepiped_structure(dim_type nc, dim_type k)
Give a pointer on the structures of a parallelepiped of dimension d.
Definition: bgeot_convex_structure.cc:397
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:55
getfem::outer_faces_of_mesh
void APIDECL outer_faces_of_mesh(const mesh &m, const dal::bit_vector &cvlst, convex_face_ct &flist)
returns a list of "exterior" faces of a mesh (i.e.
Definition: getfem_mesh.cc:822
bgeot::mesh_structure::nb_points_of_convex
short_type nb_points_of_convex(size_type ic) const
Return the number of points of convex ic.
Definition: bgeot_mesh_structure.h:115
getfem::dx_export
A (quite large) class for exportation of data to IBM OpenDX.
Definition: getfem_export.h:318
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1268
getfem::mesh_level_set::linked_mesh
mesh & linked_mesh(void) const
Gives a reference to the linked mesh of type mesh.
Definition: getfem_mesh_level_set.h:99
getfem::dx_export::exporting_mesh_edges
void exporting_mesh_edges(bool with_slice_edge=true)
append edges information (if you want to draw the mesh and are using a refined slice.
Definition: getfem_export.cc:1023
bgeot::kdtree::add_point_with_id
void add_point_with_id(const base_node &n, size_type i)
insert a new point, with an associated number.
Definition: bgeot_kdtree.h:120