GetFEM++  5.3
getfem_mesh_slicers.cc
1 /*===========================================================================
2 
3  Copyright (C) 2004-2017 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 
22 #include "getfem/dal_singleton.h"
27 
28 namespace getfem {
29  const float slicer_action::EPS = float(1e-13);
30 
31  /* ------------------------------ slicers ------------------------------- */
32 
33  slicer_none& slicer_none::static_instance() {
35  }
36 
37  /* boundary extraction */
38  slicer_boundary::slicer_boundary(const mesh& m, slicer_action &sA,
39  const mesh_region& cvflst) : A(&sA) {
40  build_from(m,cvflst);
41  }
42 
43  slicer_boundary::slicer_boundary(const mesh& m, slicer_action &sA) : A(&sA) {
44  mesh_region cvflist;
45  outer_faces_of_mesh(m, cvflist);
46  build_from(m,cvflist);
47  }
48 
49  void slicer_boundary::build_from(const mesh& m, const mesh_region& cvflst) {
50  if (m.convex_index().card()==0) return;
51  convex_faces.resize(m.convex_index().last()+1, slice_node::faces_ct(0L));
52  for (mr_visitor i(cvflst); !i.finished(); ++i)
53  if (i.is_face()) convex_faces[i.cv()][i.f()]=1;
54  else convex_faces[i.cv()].set();
55  /* set the mask to 1 for all other possible faces of the convexes, which may
56  appear after slicing the convex, hence they will be part of the "boundary" */
57  for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
58  for (short_type f=m.structure_of_convex(cv)->nb_faces(); f < convex_faces[cv].size(); ++f)
59  convex_faces[cv][f]=1;
60  }
61  }
62 
63  bool slicer_boundary::test_bound(const slice_simplex& s, slice_node::faces_ct& fmask, const mesh_slicer::cs_nodes_ct& nodes) const {
64  slice_node::faces_ct f; f.set();
65  for (size_type i=0; i < s.dim()+1; ++i) {
66  f &= nodes[s.inodes[i]].faces;
67  }
68  f &= fmask;
69  return (f.any());
70  }
71 
72  void slicer_boundary::exec(mesh_slicer& ms) {
73  if (A) A->exec(ms);
74  if (ms.splx_in.card() == 0) return;
75  slice_node::faces_ct fmask(ms.cv < convex_faces.size() ? convex_faces[ms.cv] : 0);
76  /* quickly check if the convex have any chance to be part of the boundary */
77  if (!convex_faces[ms.cv].any()) { ms.splx_in.clear(); return; }
78 
79  for (dal::bv_visitor_c cnt(ms.splx_in); !cnt.finished(); ++cnt) {
80  const slice_simplex& s = ms.simplexes[cnt];
81  if (s.dim() < ms.nodes[0].pt.size()) {
82  if (!test_bound(s, fmask, ms.nodes)) ms.splx_in.sup(cnt);
83  } else if (s.dim() == 2) {
84  ms.sup_simplex(cnt);
85  slice_simplex s2(2);
86  for (size_type j=0; j < 3; ++j) {
87  /* usage of s forbidden in this loop since push_back happens .. */
88  static unsigned ord[][2] = {{0,1},{1,2},{2,0}}; /* keep orientation of faces */
89  for (size_type k=0; k < 2; ++k) { s2.inodes[k] = ms.simplexes[cnt].inodes[ord[j][k]]; }
90  if (test_bound(s2, fmask, ms.nodes)) {
91  ms.add_simplex(s2, true);
92  }
93  }
94  } else if (s.dim() == 3) {
95  //ms.simplex_orientation(ms.simplexes[cnt]);
96  ms.sup_simplex(cnt);
97  slice_simplex s2(3);
98  for (size_type j=0; j < 4; ++j) {
99  /* usage of s forbidden in this loop since push_back happens .. */
100  static unsigned ord[][3] = {{0,2,1},{1,2,3},{1,3,0},{0,3,2}}; /* keep orientation of faces */
101  for (size_type k=0; k < 3; ++k) { s2.inodes[k] = ms.simplexes[cnt].inodes[ord[j][k]]; } //k + (k<j ? 0 : 1)]; }
102  /*cerr << " -> testing "; for (size_type iA=0; iA < s2.dim()+1; ++iA) cerr << s2.inodes[iA] << " ";
103  cerr << " : " << test_bound(s2, fmask, nodes) << endl;*/
104  if (test_bound(s2, fmask, ms.nodes)) {
105  ms.add_simplex(s2, true);
106  }
107  }
108  } /* simplexes of higher dimension are ignored... */
109  }
110  ms.update_nodes_index();
111  }
112 
113  /* apply deformation from a mesh_fem to the nodes */
114  void slicer_apply_deformation::exec(mesh_slicer& ms) {
115  base_vector coeff;
116  base_matrix G;
117  bool ref_pts_changed = false;
118  if (ms.cvr != ms.prev_cvr
119  || defdata->pmf->fem_of_element(ms.cv) != pf) {
120  pf = defdata->pmf->fem_of_element(ms.cv);
121  if (pf->need_G())
122  bgeot::vectors_to_base_matrix
123  (G, defdata->pmf->linked_mesh().points_of_convex(ms.cv));
124  }
125  /* check that the points are still the same
126  * -- or recompute the fem_precomp */
127  std::vector<base_node> ref_pts2; ref_pts2.reserve(ms.nodes_index.card());
128  for (dal::bv_visitor i(ms.nodes_index); !i.finished(); ++i) {
129  ref_pts2.push_back(ms.nodes[i].pt_ref);
130  if (ref_pts2.size() > ref_pts.size()
131  || gmm::vect_dist2_sqr(ref_pts2[i],ref_pts[i])>1e-20)
132  ref_pts_changed = true;
133  }
134  if (ref_pts2.size() != ref_pts.size()) ref_pts_changed = true;
135  if (ref_pts_changed) {
136  ref_pts.swap(ref_pts2);
137  fprecomp.clear();
138  }
139  bgeot::pstored_point_tab pspt = store_point_tab(ref_pts);
140  pfem_precomp pfp = fprecomp(pf, pspt);
141  defdata->copy(ms.cv, coeff);
142 
143  base_vector val(ms.m.dim());
144  size_type cnt = 0;
145  fem_interpolation_context ctx(ms.pgt, pfp, 0, G, ms.cv, short_type(-1));
146  for (dal::bv_visitor i(ms.nodes_index); !i.finished(); ++i, ++cnt) {
147  ms.nodes[i].pt.resize(defdata->pmf->get_qdim());
148  ctx.set_ii(cnt);
149  pf->interpolation(ctx, coeff, val, defdata->pmf->get_qdim());
150  gmm::add(val, ms.nodes[cnt].pt);
151  }
152  }
153 
154 
155  //static bool check_flat_simplex(mesh_slicer::cs_nodes_ct& /*nodes*/, const slice_simplex /*s*/) {
156  /*base_matrix M(s.dim(),s.dim());
157  for (size_type i=0; i < s.dim(); ++i) {
158  for (size_type j=0; j < s.dim(); ++j) {
159  M(i,j) = nodes[s.inodes[i+1]].pt_ref[j] - nodes[s.inodes[0]].pt_ref[j];
160  }
161  }
162  scalar_type d = gmm::lu_det(M);
163  if (gmm::abs(d) < pow(1e-6,s.dim())) {
164  cout.precision(10);
165  cout << "!!Flat (" << d << ") :";
166  for (size_type i=0; i < s.dim()+1; ++i) cout << " " << nodes[s.inodes[i]].pt;
167  cout << "\n";
168  return false;
169  }*/
170  // return true;
171  //}
172 
173  /*
174  intersects the simplex with the slice, and (recursively)
175  decomposes it into sub-simplices, which are added to the list
176  'splxs'. If orient == 0, then it is the faces of sub-simplices
177  which are added
178 
179  assertion: when called, it will always push *at least* one new
180  simplex on the stack.
181 
182  remark: s is not reference: on purpose.
183  */
184  void slicer_volume::split_simplex(mesh_slicer& ms,
185  slice_simplex s, size_type sstart,
186  std::bitset<32> spin, std::bitset<32> spbin) {
187  scalar_type alpha = 0; size_type iA=0, iB = 0;
188  bool intersection = false;
189  static int level = 0;
190 
191  level++;
192  /*
193  cerr << "split_simplex: level=" << level << " simplex: ";
194  for (iA=0; iA < s.dim()+1 && !intersection; ++iA)
195  cerr << "node#" << s.inodes[iA] << "=" << nodes[s.inodes[iA]].pt
196  << ", in=" << pt_in[s.inodes[iA]] << ", bin=" << pt_bin[s.inodes[iA]] << "; "; cerr << endl;
197  */
198  assert(level < 100);
199  for (iA=0; iA < s.dim(); ++iA) {
200  if (spbin[iA]) continue;
201  for (iB=iA+1; iB < s.dim()+1; ++iB) {
202  if (!spbin[iB] && spin[iA] != spin[iB]) {
203  alpha=edge_intersect(s.inodes[iA],s.inodes[iB],ms.nodes);
204  if (alpha >= 1e-8 && alpha <= 1-1e-8) { intersection = true; break; }
205  }
206  }
207  if (intersection) break;
208  }
209  if (intersection) {
210  /* will call split_simplex recursively */
211  const slice_node& A = ms.nodes[s.inodes[iA]];
212  const slice_node& B = ms.nodes[s.inodes[iB]];
213  slice_node n(A.pt + alpha*(B.pt-A.pt), A.pt_ref + alpha*(B.pt_ref-A.pt_ref));
214  n.faces = A.faces & B.faces;
215  size_type nn = ms.nodes.size();
216  ms.nodes.push_back(n); /* invalidate A and B.. */
217  pt_bin.add(nn); pt_in.add(nn);
218 
219  std::bitset<32> spin2(spin), spbin2(spbin);
220  std::swap(s.inodes[iA],nn);
221  spin2.set(iA); spbin2.set(iA);
222  split_simplex(ms, s, sstart, spin2, spbin2);
223 
224  std::swap(s.inodes[iA],nn); std::swap(s.inodes[iB],nn);
225  spin2 = spin; spbin2 = spbin; spin2.set(iB); spbin2.set(iB);
226  split_simplex(ms, s, sstart, spin2, spbin2);
227 
228  } else {
229  /* end of recursion .. */
230  bool all_in = true;
231  for (size_type i=0; i < s.dim()+1; ++i) if (!spin[i]) { all_in = false; break; }
232  //check_flat_simplex(ms.nodes,s);
233 
234  // even simplexes "outside" are pushed, in case of a slicer_complementary op
235  ms.add_simplex(s,(all_in && orient != VOLBOUND) || orient == VOLSPLIT);
236  if (orient==0) { /* reduce dimension */
237  slice_simplex face(s.dim());
238  for (size_type f=0; f < s.dim()+1; ++f) {
239  all_in = true;
240  for (size_type i=0; i < s.dim(); ++i) {
241  size_type p = i + (i<f?0:1);
242  if (!spbin[p]) { all_in = false; break; }
243  else face.inodes[i] = s.inodes[p];
244  }
245  if (all_in) {
246  /* prevent addition of a face twice */
247  std::sort(face.inodes.begin(), face.inodes.end());
248  if (std::find(ms.simplexes.begin()+sstart, ms.simplexes.end(), face) == ms.simplexes.end()) {
249  ms.add_simplex(face,true);
250  }
251  }
252  }
253  }
254  }
255  level--;
256  }
257 
258  /* nodes : list of nodes (new nodes may be added)
259  splxs : list of simplexes (new simplexes may be added)
260  splx_in : input: simplexes to take into account, output: list of simplexes inside the slice
261 
262  note that the simplexes in the list may have different dimensions
263  */
264  void slicer_volume::exec(mesh_slicer& ms) {
265  //cerr << "\n----\nslicer_volume::slice : entree, splx_in=" << splx_in << endl;
266  if (ms.splx_in.card() == 0) return;
267  prepare(ms.cv,ms.nodes,ms.nodes_index);
268  for (dal::bv_visitor_c cnt(ms.splx_in); !cnt.finished(); ++cnt) {
269  slice_simplex& s = ms.simplexes[cnt];
270  /*cerr << "\n--------slicer_volume::slice : slicing convex " << cnt << endl;
271  for (size_type i=0; i < s.dim()+1; ++i)
272  cerr << " * pt[" << i << "]=" << nodes[s.inodes[i]].pt << ", is_in=" <<
273  is_in(nodes[s.inodes[i]].pt) << ", is_bin=" << is_in(nodes[s.inodes[i]].pt,true) << endl;
274  */
275  size_type in_cnt = 0, in_bcnt = 0;
276  std::bitset<32> spin, spbin;
277  for (size_type i=0; i < s.dim()+1; ++i) {
278  if (pt_in.is_in(s.inodes[i])) { ++in_cnt; spin.set(i); }
279  if (pt_bin.is_in(s.inodes[i])) { ++in_bcnt; spbin.set(i); }
280  }
281 
282  if (in_cnt == 0) {
283  if (orient != VOLSPLIT) ms.splx_in.sup(cnt);
284  } else if (in_cnt != s.dim()+1 || orient==VOLBOUND) { /* the simplex crosses the slice boundary */
285  ms.sup_simplex(cnt);
286  split_simplex(ms, s, ms.simplexes.size(), spin, spbin);
287  }
288  }
289 
290  /* signalement des points qui se trouvent pile-poil sur la bordure */
291  if (pt_bin.card()) {
292  GMM_ASSERT1(ms.fcnt != dim_type(-1),
293  "too much {faces}/{slices faces} in the convex " << ms.cv
294  << " (nbfaces=" << ms.fcnt << ")");
295  for (dal::bv_visitor cnt(pt_bin); !cnt.finished(); ++cnt) {
296  ms.nodes[cnt].faces.set(ms.fcnt);
297  }
298  ms.fcnt++;
299  }
300  ms.update_nodes_index();
301  }
302 
303  slicer_mesh_with_mesh::slicer_mesh_with_mesh(const mesh& slm_) : slm(slm_) {
304  base_node min,max;
305  for (dal::bv_visitor cv(slm.convex_index()); !cv.finished(); ++cv) {
306  bgeot::bounding_box(min,max,slm.points_of_convex(cv),slm.trans_of_convex(cv));
307  tree.add_box(min, max, cv);
308  }
309  }
310 
311  void slicer_mesh_with_mesh::exec(mesh_slicer &ms) {
312  /* identify potientially intersecting convexes of slm */
313  base_node min(ms.nodes[0].pt),max(ms.nodes[0].pt);
314  for (size_type i=1; i < ms.nodes.size(); ++i) {
315  for (size_type k=0; k < min.size(); ++k) {
316  min[k] = std::min(min[k], ms.nodes[i].pt[k]);
317  max[k] = std::max(max[k], ms.nodes[i].pt[k]);
318  }
319  }
320  std::vector<size_type> slmcvs;
321  tree.find_intersecting_boxes(min, max, slmcvs);
322  /* save context */
323  mesh_slicer::cs_simplexes_ct simplexes_final(ms.simplexes);
324  dal::bit_vector splx_in_save(ms.splx_in),
325  simplex_index_save(ms.simplex_index), nodes_index_save(ms.nodes_index);
326  size_type scnt0 = ms.simplexes.size();
327  /* loop over candidate convexes of slm */
328  //cout << "slicer_mesh_with_mesh: convex " << ms.cv << ", " << ms.splx_in.card() << " candidates\n";
329  for (size_type i=0; i < slmcvs.size(); ++i) {
330  size_type slmcv = slmcvs[i];
331  dim_type fcnt_save = dim_type(ms.fcnt);
332  ms.simplexes.resize(scnt0);
333  ms.splx_in = splx_in_save; ms.simplex_index = simplex_index_save; ms.nodes_index = nodes_index_save;
334  //cout << "test intersection of " << ms.cv << " and " << slmcv << "\n";
335  /* loop over the faces and apply slicer_half_space */
336  for (short_type f=0; f<slm.structure_of_convex(slmcv)->nb_faces(); ++f) {
337  base_node x0,n;
338  if (slm.structure_of_convex(slmcv)->dim() == 3 && slm.dim() == 3) {
339  x0 = slm.points_of_face_of_convex(slmcv,f)[0];
340  base_node A = slm.points_of_face_of_convex(slmcv,f)[1] - x0;
341  base_node B = slm.points_of_face_of_convex(slmcv,f)[2] - x0;
342  base_node G = gmm::mean_value(slm.points_of_convex(slmcv).begin(),slm.points_of_convex(slmcv).end());
343  n.resize(3);
344  n[0] = A[1]*B[2] - A[2]*B[1];
345  n[1] = A[2]*B[0] - A[0]*B[2];
346  n[2] = A[0]*B[1] - A[1]*B[0];
347  if (gmm::vect_sp(n,G-x0) > 0) n *= -1.;
348  } else {
349  size_type ip = slm.structure_of_convex(slmcv)->nb_points_of_face(f) / 2;
350  x0 = slm.points_of_face_of_convex(slmcv,f)[ip];
351  n = slm.normal_of_face_of_convex(slmcv,f, x0);
352  }
353  slicer_half_space slf(x0,n,slicer_volume::VOLIN);
354  slf.exec(ms);
355  if (ms.splx_in.card() == 0) break;
356  }
357  if (ms.splx_in.card()) intersection_callback(ms, slmcv);
358  for (dal::bv_visitor is(ms.splx_in); !is.finished(); ++is) {
359  simplexes_final.push_back(ms.simplexes[is]);
360  }
361  ms.fcnt=fcnt_save;
362  }
363  ms.splx_in.clear(); ms.splx_in.add(scnt0, simplexes_final.size()-scnt0);
364  ms.simplexes.swap(simplexes_final);
365  ms.simplex_index = ms.splx_in;
366  ms.update_nodes_index();
367  /*cout << "convex " << ms.cv << "was sliced into " << ms.splx_in.card()
368  << " simplexes, nodes.size=" << ms.nodes.size()
369  << ", used nodes=" << ms.nodes_index.card() << "\n";*/
370  }
371 
372  /* isosurface computations */
373  void slicer_isovalues::prepare(size_type cv,
374  const mesh_slicer::cs_nodes_ct& nodes,
375  const dal::bit_vector& nodes_index) {
376  pt_in.clear(); pt_bin.clear();
377  std::vector<base_node> refpts(nodes.size());
378  Uval.resize(nodes.size());
379  base_vector coeff;
380  base_matrix G;
381  pfem pf = mfU->pmf->fem_of_element(cv);
382  if (pf == 0) return;
383  fem_precomp_pool fprecomp;
384  if (pf->need_G())
385  bgeot::vectors_to_base_matrix
386  (G,mfU->pmf->linked_mesh().points_of_convex(cv));
387  for (size_type i=0; i < nodes.size(); ++i) refpts[i] = nodes[i].pt_ref;
388  pfem_precomp pfp = fprecomp(pf, store_point_tab(refpts));
389  mfU->copy(cv, coeff);
390  //cerr << "cv=" << cv << ", val=" << val << ", coeff=" << coeff << endl;
391  base_vector v(1);
392  fem_interpolation_context ctx(mfU->pmf->linked_mesh().trans_of_convex(cv),
393  pfp, 0, G, cv, short_type(-1));
394  for (dal::bv_visitor i(nodes_index); !i.finished(); ++i) {
395  v[0] = 0;
396  ctx.set_ii(i);
397  pf->interpolation(ctx, coeff, v, mfU->pmf->get_qdim());
398  Uval[i] = v[0];
399  // optimisable -- les bit_vectors sont lents..
400  pt_bin[i] = (gmm::abs(Uval[i] - val) < EPS * val_scaling);
401  pt_in[i] = (Uval[i] - val < 0); if (orient>0) pt_in[i] = !pt_in[i];
402  pt_in[i] = pt_in[i] || pt_bin[i];
403  // cerr << "cv=" << cv << ", node["<< i << "]=" << nodes[i].pt
404  // << ", Uval[i]=" << Uval[i] << ", pt_in[i]=" << pt_in[i]
405  // << ", pt_bin[i]=" << pt_bin[i] << endl;
406  }
407  }
408 
409 
410  void slicer_union::exec(mesh_slicer &ms) {
411  dal::bit_vector splx_in_base = ms.splx_in;
412  size_type c = ms.simplexes.size();
413  dim_type fcnt_0 = dim_type(ms.fcnt);
414  A->exec(ms);
415  dal::bit_vector splx_inA(ms.splx_in);
416  dim_type fcnt_1 = dim_type(ms.fcnt);
417 
418  dal::bit_vector splx_inB = splx_in_base;
419  splx_inB.add(c, ms.simplexes.size()-c);
420  splx_inB.setminus(splx_inA);
421  for (dal::bv_visitor_c i(splx_inB); !i.finished(); ++i) {
422  if (!ms.simplex_index[i]) splx_inB.sup(i);
423  }
424  //splx_inB &= ms.simplex_index;
425  ms.splx_in = splx_inB;
426  B->exec(ms); splx_inB = ms.splx_in;
427  ms.splx_in |= splx_inA;
428 
429  /*
430  the boring part : making sure that the "slice face" node markers
431  are correctly set
432  */
433  for (unsigned f=fcnt_0; f < fcnt_1; ++f) {
434  for (dal::bv_visitor i(splx_inB); !i.finished(); ++i) {
435  for (unsigned j=0; j < ms.simplexes[i].dim()+1; ++j) {
436  bool face_boundA = true;
437  for (unsigned k=0; k < ms.simplexes[i].dim()+1; ++k) {
438  if (j != k && !ms.nodes[ms.simplexes[i].inodes[k]].faces[f]) {
439  face_boundA = false; break;
440  }
441  }
442  if (face_boundA) {
443  /* now we know that the face was on slice A boundary, and
444  that the convex is inside slice B. The conclusion: the
445  face is not on a face of A union B.
446  */
447  for (unsigned k=0; k < ms.simplexes[i].dim()+1; ++k)
448  if (j != k) ms.nodes[ms.simplexes[i].inodes[k]].faces[f] = false;
449  }
450  }
451  }
452  }
453  ms.update_nodes_index();
454  }
455 
456  void slicer_intersect::exec(mesh_slicer& ms) {
457  A->exec(ms);
458  B->exec(ms);
459  }
460 
461  void slicer_complementary::exec(mesh_slicer& ms) {
462  dal::bit_vector splx_inA = ms.splx_in;
463  size_type sz = ms.simplexes.size();
464  A->exec(ms); splx_inA.swap(ms.splx_in);
465  ms.splx_in &= ms.simplex_index;
466  dal::bit_vector bv = ms.splx_in; bv.add(sz, ms.simplexes.size()-sz); bv &= ms.simplex_index;
467  for (dal::bv_visitor_c i(bv); !i.finished(); ++i) {
468  /*cerr << "convex " << cv << ",examining simplex #" << i << ": {";
469  for (size_type j=0; j < simplexes[i].inodes.size(); ++j) cerr << nodes[simplexes[i].inodes[j]].pt << " ";
470  cerr << "}, splx_in=" << splx_in[i] << "splx_inA=" << splx_inA[i] << endl;*/
471  ms.splx_in[i] = !splx_inA.is_in(i);
472  }
473  }
474 
475  void slicer_compute_area::exec(mesh_slicer &ms) {
476  for (dal::bv_visitor is(ms.splx_in); !is.finished(); ++is) {
477  const slice_simplex& s = ms.simplexes[is];
478  base_matrix M(s.dim(),s.dim());
479  for (size_type i=0; i < s.dim(); ++i)
480  for (size_type j=0; j < s.dim(); ++j)
481  M(i,j) = ms.nodes[s.inodes[i+1]].pt[j] - ms.nodes[s.inodes[0]].pt[j];
482  scalar_type v = bgeot::lu_det(&(*(M.begin())), s.dim());
483  for (size_type d=2; d <= s.dim(); ++d) v /= scalar_type(d);
484  a += v;
485  }
486  }
487 
488  void slicer_build_edges_mesh::exec(mesh_slicer &ms) {
489  for (dal::bv_visitor is(ms.splx_in); !is.finished(); ++is) {
490  const slice_simplex& s = ms.simplexes[is];
491  for (size_type i=0; i < s.dim(); ++i) {
492  for (size_type j=i+1; j <= s.dim(); ++j) {
493  const slice_node& A = ms.nodes[s.inodes[i]];
494  const slice_node& B = ms.nodes[s.inodes[j]];
495  /* duplicate with stored_mesh_slice which also
496  builds a list of edges */
497  if ((A.faces & B.faces).count() >= unsigned(ms.cv_dim-1)) {
498  slice_node::faces_ct fmask((1 << ms.cv_nbfaces)-1); fmask.flip();
499  size_type e = edges_m.add_segment_by_points(A.pt,B.pt);
500  if (pslice_edges && (((A.faces & B.faces) & fmask).any())) pslice_edges->add(e);
501  }
502  }
503  }
504  }
505  }
506 
507  void slicer_build_mesh::exec(mesh_slicer &ms) {
508  std::vector<size_type> pid(ms.nodes_index.last_true()+1);
509  for (dal::bv_visitor i(ms.nodes_index); !i.finished(); ++i)
510  pid[i] = m.add_point(ms.nodes[i].pt);
511  for (dal::bv_visitor i(ms.splx_in); !i.finished(); ++i) {
512  for (unsigned j=0; j < ms.simplexes.at(i).inodes.size(); ++j) {
513  assert(m.points_index().is_in(pid.at(ms.simplexes.at(i).inodes[j])));
514  }
515  m.add_convex(bgeot::simplex_geotrans(ms.simplexes[i].dim(),1),
516  gmm::index_ref_iterator(pid.begin(),
517  ms.simplexes[i].inodes.begin()));
518  }
519  }
520 
521  void slicer_explode::exec(mesh_slicer &ms) {
522  if (ms.nodes_index.card() == 0) return;
523 
524  base_node G;
525  if (ms.face < dim_type(-1))
526  G = gmm::mean_value(ms.m.points_of_face_of_convex(ms.cv, ms.face).begin(),
527  ms.m.points_of_face_of_convex(ms.cv, ms.face).end());
528  else
529  G = gmm::mean_value(ms.m.points_of_convex(ms.cv).begin(),
530  ms.m.points_of_convex(ms.cv).end());
531  for (dal::bv_visitor i(ms.nodes_index); !i.finished(); ++i)
532  ms.nodes[i].pt = G + coef*(ms.nodes[i].pt - G);
533 
534  for (dal::bv_visitor cnt(ms.splx_in); !cnt.finished(); ++cnt) {
535  const slice_simplex& s = ms.simplexes[cnt];
536  if (s.dim() == 3) { // keep only faces
537  ms.sup_simplex(cnt);
538  slice_simplex s2(3);
539  for (size_type j=0; j < 4; ++j) {
540  /* usage of s forbidden in this loop since push_back happens .. */
541  static unsigned ord[][3] = {{0,2,1},{1,2,3},{1,3,0},{0,3,2}}; /* keep orientation of faces */
542  for (size_type k=0; k < 3; ++k) { s2.inodes[k] = ms.simplexes[cnt].inodes[ord[j][k]]; } //k + (k<j ? 0 : 1)]; }
543 
544  slice_node::faces_ct f; f.set();
545  for (size_type i=0; i < s2.dim()+1; ++i) {
546  f &= ms.nodes[s2.inodes[i]].faces;
547  }
548  if (f.any()) {
549  ms.add_simplex(s2, true);
550  }
551  }
552  }
553  }
554  }
555 
556  /* -------------------- member functions of mesh_slicer -------------- */
557 
558  mesh_slicer::mesh_slicer(const mesh_level_set &mls_) :
559  m(mls_.linked_mesh()), mls(&mls_), pgt(0), cvr(0) {}
561  m(m_), mls(0), pgt(0), cvr(0) {}
562 
563  void mesh_slicer::using_mesh_level_set(const mesh_level_set &mls_) {
564  mls = &mls_;
565  GMM_ASSERT1(&m == &mls->linked_mesh(), "different meshes");
566  }
567 
568  void mesh_slicer::pack() {
569  std::vector<size_type> new_id(nodes.size());
570  size_type ncnt = 0;
571  for (dal::bv_visitor i(nodes_index); !i.finished(); ++i) {
572  if (i != ncnt) {
573  nodes[i].swap(nodes[ncnt]);
574  }
575  new_id[i] = ncnt++;
576  }
577  nodes.resize(ncnt);
578  size_type scnt = 0;
579  for (dal::bv_visitor j(simplex_index); !j.finished(); ++j) {
580  if (j != scnt) { simplexes[scnt] = simplexes[j]; }
581  for (std::vector<size_type>::iterator it = simplexes[scnt].inodes.begin();
582  it != simplexes[scnt].inodes.end(); ++it) {
583  *it = new_id[*it];
584  }
585  }
586  simplexes.resize(scnt);
587  }
588 
589  void mesh_slicer::update_nodes_index() {
590  nodes_index.clear();
591  for (dal::bv_visitor j(simplex_index); !j.finished(); ++j) {
592  assert(j < simplexes.size());
593  for (std::vector<size_type>::iterator it = simplexes[j].inodes.begin();
594  it != simplexes[j].inodes.end(); ++it) {
595  assert(*it < nodes.size());
596  nodes_index.add(*it);
597  }
598  }
599  }
600 
601  static void flag_points_on_faces(const bgeot::pconvex_ref& cvr,
602  const std::vector<base_node>& pts,
603  std::vector<slice_node::faces_ct>& faces) {
604  GMM_ASSERT1(cvr->structure()->nb_faces() <= 32,
605  "won't work for convexes with more than 32 faces "
606  "(hardcoded limit)");
607  faces.resize(pts.size());
608  for (size_type i=0; i < pts.size(); ++i) {
609  faces[i].reset();
610  for (short_type f=0; f < cvr->structure()->nb_faces(); ++f)
611  faces[i][f] = (gmm::abs(cvr->is_in_face(f, pts[i])) < 1e-10);
612  }
613  }
614 
615  void mesh_slicer::update_cv_data(size_type cv_, short_type f_) {
616  cv = cv_;
617  face = f_;
618  pgt = m.trans_of_convex(cv);
619  prev_cvr = cvr; cvr = pgt->convex_ref();
620  cv_dim = cvr->structure()->dim();
621  cv_nbfaces = cvr->structure()->nb_faces();
622  fcnt = cvr->structure()->nb_faces();
623  discont = (mls && mls->is_convex_cut(cv));
624  }
625 
626  void mesh_slicer::apply_slicers() {
627  simplex_index.clear(); simplex_index.add(0, simplexes.size());
628  splx_in = simplex_index;
629  nodes_index.clear(); nodes_index.add(0, nodes.size());
630  for (size_type i=0; i < action.size(); ++i) {
631  action[i]->exec(*this);
632  //cout << "simplex_index=" << simplex_index << "\n splx_in=" << splx_in << "\n";
633  assert(simplex_index.contains(splx_in));
634  }
635  }
636 
637  void mesh_slicer::simplex_orientation(slice_simplex& s) {
638  size_type N = m.dim();
639  if (s.dim() == N) {
640  base_matrix M(N,N);
641  for (size_type i=1; i < N+1; ++i) {
642  base_small_vector d = nodes[s.inodes[i]].pt - nodes[s.inodes[0]].pt;
643  gmm::copy_n(d.const_begin(), N, M.begin() + (i-1)*N);
644  }
645  scalar_type J = bgeot::lu_det(&(*(M.begin())), N);
646  //cout << " lu_det = " << J << "\n";
647  if (J < 0) {
648  std::swap(s.inodes[1],s.inodes[0]);
649  }
650  }
651  }
652 
653  void mesh_slicer::exec(size_type nrefine, const mesh_region& cvlst) {
654  short_type n = short_type(nrefine);
655  exec_(&n, 0, cvlst);
656  }
657 
658  void mesh_slicer::exec(const std::vector<short_type> &nrefine,
659  const mesh_region& cvlst) {
660  exec_(&nrefine[0], 1, cvlst);
661  }
662 
663  static bool check_orient(size_type cv, bgeot::pgeometric_trans pgt,
664  const mesh& m) {
665  if (pgt->dim() == m.dim() && m.dim()>=2) { /* no orient check for
666  convexes of lower dim */
667  base_matrix G; bgeot::vectors_to_base_matrix(G,m.points_of_convex(cv));
668  base_node g(pgt->dim()); g.fill(.5);
669  base_matrix pc; pgt->poly_vector_grad(g,pc);
670  base_matrix K(pgt->dim(),pgt->dim());
671  gmm::mult(G,pc,K);
672  scalar_type J = bgeot::lu_det(&(*(K.begin())), pgt->dim());
673  // bgeot::geotrans_interpolation_context ctx(pgp,0,G);
674  // scalar_type J = gmm::lu_det(ctx.B()); // pb car inverse K mÍme
675  if (J < 0) return true;
676  //cout << "cv = " << cv << ", J = " << J << "\n";
677  }
678  return false;
679  }
680 
681 #if OLD_MESH_SLICER
682  void mesh_slicer::exec_(const short_type *pnrefine, int nref_stride,
683  const mesh_region& cvlst) {
684  std::vector<base_node> cvm_pts;
685  const bgeot::basic_mesh *cvm = 0;
686  const bgeot::mesh_structure *cvms = 0;
688  bgeot::pgeotrans_precomp pgp = 0;
689  std::vector<slice_node::faces_ct> points_on_faces;
690 
691  cvlst.from_mesh(m);
692  size_type prev_nrefine = 0;
693  for (mr_visitor it(cvlst); !it.finished(); ++it) {
694  size_type nrefine = pnrefine[it.cv()*nref_stride];
695  update_cv_data(it.cv(),it.f());
696  bool revert_orientation = check_orient(cv, pgt,m);
697 
698  /* update structure-dependent data */
699  if (prev_cvr != cvr || nrefine != prev_nrefine) {
700  cvm = bgeot::refined_simplex_mesh_for_convex(cvr, nrefine);
701  cvm_pts.resize(cvm->points().card());
702  std::copy(cvm->points().begin(), cvm->points().end(), cvm_pts.begin());
703  pgp = gppool(pgt, store_point_tab(cvm_pts));
704  flag_points_on_faces(cvr, cvm_pts, points_on_faces);
705  prev_nrefine = nrefine;
706  }
707  if (face < dim_type(-1))
708  cvms = bgeot::refined_simplex_mesh_for_convex_faces(cvr, nrefine)[face].get();
709  else
710  cvms = cvm;
711 
712  /* apply the initial geometric transformation */
713  std::vector<size_type> ptsid(cvm_pts.size()); std::fill(ptsid.begin(), ptsid.end(), size_type(-1));
714  simplexes.resize(cvms->nb_convex());
715  nodes.resize(0);
716  for (size_type snum = 0; snum < cvms->nb_convex(); ++snum) {
717  /* cvms should not contain holes in its convex index.. */
718  simplexes[snum].inodes.resize(cvms->nb_points_of_convex(snum));
719  std::copy(cvms->ind_points_of_convex(snum).begin(),
720  cvms->ind_points_of_convex(snum).end(), simplexes[snum].inodes.begin());
721  if (revert_orientation) std::swap(simplexes[snum].inodes[0],simplexes[snum].inodes[1]);
722  /* store indices of points which are really used , and renumbers them */
723  for (std::vector<size_type>::iterator itp = simplexes[snum].inodes.begin();
724  itp != simplexes[snum].inodes.end(); ++itp) {
725  if (ptsid[*itp] == size_type(-1)) {
726  ptsid[*itp] = nodes.size();
727  nodes.push_back(slice_node());
728  nodes.back().pt_ref = cvm_pts[*itp];
729  nodes.back().faces = points_on_faces[*itp];
730  nodes.back().pt.resize(m.dim()); nodes.back().pt.fill(0.);
731  pgp->transform(m.points_of_convex(cv), *itp, nodes.back().pt);
732  }
733  *itp = ptsid[*itp];
734  }
735  if (0) {
736  static int once = 0;
737  if (once++ < 3) {
738  cout << "check orient cv " << cv << ", snum=" << snum << "/" << cvms->nb_convex();
739  }
740  simplex_orientation(simplexes[snum]);
741  }
742  }
743  apply_slicers();
744  }
745  }
746 #endif // OLD_MESH_SLICER
747 
748  template <typename CONT>
749  static void add_degree1_convex(bgeot::pgeometric_trans pgt, const CONT &pts,
750  mesh &m) {
751  unsigned N = pgt->dim();
752  std::vector<base_node> v; v.reserve(N+1);
753  for (unsigned i=0; i < pgt->nb_points(); ++i) {
754  const base_node P = pgt->convex_ref()->points()[i];
755  scalar_type s = 0;
756  for (unsigned j=0; j < N; ++j) {
757  s += P[j]; if (P[j] == 1) { v.push_back(pts[i]); break; }
758  }
759  if (s == 0) v.push_back(pts[i]);
760  }
761  assert(v.size() == N+1);
762  base_node G = gmm::mean_value(v);
763  /*for (unsigned i=0; i < v.size();++i)
764  v[i] = v[i] + 0.1 * (G - v[i]);*/
765  m.add_convex_by_points(bgeot::simplex_geotrans(N,1), v.begin());
766  }
767 
768  const mesh& mesh_slicer::refined_simplex_mesh_for_convex_cut_by_level_set
769  (const mesh &cvm, unsigned nrefine) {
770  mesh mm; mm.copy_from(cvm);
771  while (nrefine > 1) {
772  mm.Bank_refine(mm.convex_index());
773  nrefine /= 2;
774  }
775 
776  std::vector<size_type> idx;
777  tmp_mesh.clear();
778  //cerr << "nb cv = " << tmp_mesh.convex_index().card() << "\n";
779  for (dal::bv_visitor_c ic(mm.convex_index()); !ic.finished(); ++ic) {
780  add_degree1_convex(mm.trans_of_convex(ic), mm.points_of_convex(ic), tmp_mesh);
781  }
782  /*tmp_mesh.write_to_file(std::cerr);
783  cerr << "\n";*/
784  return tmp_mesh;
785  }
786 
787  const bgeot::mesh_structure &
788  mesh_slicer::refined_simplex_mesh_for_convex_faces_cut_by_level_set
789  (short_type f) {
790  mesh &cvm = tmp_mesh;
791  tmp_mesh_struct.clear();
792  unsigned N = cvm.dim();
793 
794  dal::bit_vector pt_in_face; pt_in_face.sup(0, cvm.points().index().last_true()+1);
795  for (dal::bv_visitor ip(cvm.points().index()); !ip.finished(); ++ip)
796  if (gmm::abs(cvr->is_in_face(short_type(f), cvm.points()[ip]))) pt_in_face.add(ip);
797 
798  for (dal::bv_visitor_c ic(cvm.convex_index()); !ic.finished(); ++ic) {
799  for (short_type ff=0; ff < cvm.nb_faces_of_convex(ic); ++ff) {
800  bool face_ok = true;
801  for (unsigned i=0; i < N; ++i) {
802  if (!pt_in_face.is_in(cvm.ind_points_of_face_of_convex(ic,ff)[i])) {
803  face_ok = false; break;
804  }
805  }
806  if (face_ok) {
807  tmp_mesh_struct.add_convex(bgeot::simplex_structure(dim_type(N-1)),
808  cvm.ind_points_of_face_of_convex(ic, ff).begin());
809  }
810  }
811  }
812  return tmp_mesh_struct;
813  }
814 
815  void mesh_slicer::exec_(const short_type *pnrefine,
816  int nref_stride,
817  const mesh_region& cvlst) {
818  std::vector<base_node> cvm_pts;
819  const bgeot::basic_mesh *cvm = 0;
820  const bgeot::mesh_structure *cvms = 0;
822  bgeot::pgeotrans_precomp pgp = 0;
823  std::vector<slice_node::faces_ct> points_on_faces;
824  bool prev_discont = true;
825 
826  cvlst.from_mesh(m);
827  size_type prev_nrefine = 0;
828  // size_type prev_cv = size_type(-1);
829  for (mr_visitor it(cvlst); !it.finished(); ++it) {
830  size_type nrefine = pnrefine[it.cv()*nref_stride];
831  update_cv_data(it.cv(),it.f());
832  bool revert_orientation = check_orient(cv, pgt,m);
833 
834  /* update structure-dependent data */
835  /* TODO : fix levelset handling when slicing faces .. */
836  if (prev_cvr != cvr || nrefine != prev_nrefine
837  || discont || prev_discont) {
838  if (discont) {
839  cvm = &refined_simplex_mesh_for_convex_cut_by_level_set
840  (mls->mesh_of_convex(cv), unsigned(nrefine));
841  } else {
843  }
844  cvm_pts.resize(cvm->points().card());
845  std::copy(cvm->points().begin(), cvm->points().end(), cvm_pts.begin());
846  pgp = gppool(pgt, store_point_tab(cvm_pts));
847  flag_points_on_faces(cvr, cvm_pts, points_on_faces);
848  prev_nrefine = nrefine;
849  }
850  if (face < dim_type(-1)) {
851  if (!discont) {
853  (cvr, short_type(nrefine))[face].get();
854  } else {
855  cvms = &refined_simplex_mesh_for_convex_faces_cut_by_level_set(face);
856  }
857  } else {
858  cvms = cvm;
859  }
860 
861  /* apply the initial geometric transformation */
862  std::vector<size_type> ptsid(cvm_pts.size()); std::fill(ptsid.begin(), ptsid.end(), size_type(-1));
863  simplexes.resize(cvms->nb_convex());
864  nodes.resize(0);
865 
866  base_node G;
867  for (size_type snum = 0; snum < cvms->nb_convex(); ++snum) {
868  /* cvms should not contain holes in its convex index.. */
869  simplexes[snum].inodes.resize(cvms->nb_points_of_convex(snum));
870  std::copy(cvms->ind_points_of_convex(snum).begin(),
871  cvms->ind_points_of_convex(snum).end(), simplexes[snum].inodes.begin());
872  if (revert_orientation) std::swap(simplexes[snum].inodes[0],simplexes[snum].inodes[1]);
873  /* store indices of points which are really used , and renumbers them */
874  if (discont) {
875  G.resize(m.dim()); G.fill(0.);
876  for (std::vector<size_type>::iterator itp =
877  simplexes[snum].inodes.begin();
878  itp != simplexes[snum].inodes.end(); ++itp) {
879  G += cvm_pts[*itp];
880  }
881  G /= scalar_type(simplexes[snum].inodes.size());
882  }
883 
884  for (std::vector<size_type>::iterator itp =
885  simplexes[snum].inodes.begin();
886  itp != simplexes[snum].inodes.end(); ++itp) {
887  if (discont || ptsid[*itp] == size_type(-1)) {
888  ptsid[*itp] = nodes.size();
889  nodes.push_back(slice_node());
890  if (!discont) {
891  nodes.back().pt_ref = cvm_pts[*itp];
892  } else {
893  /* displace the ref point such that one will not interpolate
894  on the discontinuity (yes this is quite ugly and not
895  robust)
896  */
897  nodes.back().pt_ref = cvm_pts[*itp] + 0.01*(G - cvm_pts[*itp]);
898  }
899  nodes.back().faces = points_on_faces[*itp];
900  nodes.back().pt.resize(m.dim()); nodes.back().pt.fill(0.);
901  pgp->transform(m.points_of_convex(cv), *itp, nodes.back().pt);
902  //nodes.back().pt = pgt->transform(G, m.points_of_convex(cv));
903  //cerr << "G = " << G << " -> pt = " << nodes.back().pt << "\n";
904  }
905  *itp = ptsid[*itp];
906  }
907  }
908  //cerr << "cv = " << cv << ", cvm.nb_points_ = "<< cvm->points().size() << ", nbnodes = " << nodes.size() << ", nb_simpl=" << simplexes.size() << "\n";
909 
910  apply_slicers();
911  // prev_cv = it.cv();
912  prev_discont = discont;
913  }
914  }
915 
916 
917  void mesh_slicer::exec(size_type nrefine) {
918  exec(nrefine,mesh_region(m.convex_index()));
919  }
920 
921  /* apply slice ops to an already stored slice object */
923  GMM_ASSERT1(&sl.linked_mesh() == &m, "wrong mesh");
924  for (stored_mesh_slice::cvlst_ct::const_iterator it = sl.cvlst.begin(); it != sl.cvlst.end(); ++it) {
925  update_cv_data((*it).cv_num);
926  nodes = (*it).nodes;
927  simplexes = (*it).simplexes;
928  apply_slicers();
929  }
930  }
931 
932  /* apply slice ops to a set of nodes */
933  void mesh_slicer::exec(const std::vector<base_node>& pts) {
935  gti.add_points(pts);
938  for (dal::bv_visitor ic(m.convex_index()); !ic.finished(); ++ic) {
939  size_type nb = gti.points_in_convex(m.convex(ic), m.trans_of_convex(ic), ptab, itab);
940  if (!nb) continue;
941  update_cv_data(ic);
942  nodes.resize(0); simplexes.resize(0);
943  for (size_type i=0; i < nb; ++i) {
944  //cerr << "point " << itab[i] << "(" << pts[itab[i]]
945  //<< ") trouve dans le convex " << ic << " [pt_ref=" << ptab[i] << "]\n";
946  nodes.push_back(slice_node(pts[itab[i]],ptab[i])); nodes.back().faces=0;
947  slice_simplex s(1); s.inodes[0] = nodes.size()-1;
948  simplexes.push_back(s);
949  }
950  apply_slicers();
951  }
952  }
953 }
const mesh & linked_mesh() const
return a pointer to the original mesh
void add_points(const CONT &c)
Add the points contained in c to the list of points.
structure used to hold a set of convexes and/or convex faces.
short_type nb_points_of_convex(size_type ic) const
Return the number of points of convex ic.
const mesh_slicer::cs_nodes_ct & nodes(size_type ic) const
Return the list of nodes for the &#39;ic&#39;th convex of the slice.
Keep informations about a mesh crossed by level-sets.
The object geotrans_precomp_pool Allow to allocate a certain number of geotrans_precomp and automatic...
ind_pt_face_ct ind_points_of_face_of_convex(size_type ic, short_type f) const
Return a container of the (global) point number for face f or convex ic.
const basic_mesh * refined_simplex_mesh_for_convex(pconvex_ref cvr, short_type k)
simplexify a convex_ref.
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
ref_convex convex(size_type ic) const
return a bgeot::convex object for the convex number ic.
Definition: getfem_mesh.h:189
handles the geometric inversion for a given (supposedly quite large) set of points ...
void clear()
Erase the mesh.
Definition: getfem_mesh.cc:273
void clear()
erase the mesh
const std::vector< std::unique_ptr< mesh_structure > > & refined_simplex_mesh_for_convex_faces(pconvex_ref cvr, short_type k)
simplexify the faces of a convex_ref
static T & instance()
Instance from the current thread.
const mesh_region & from_mesh(const mesh &m) const
For regions which have been built with just a number &#39;id&#39;, from_mesh(m) sets the current region to &#39;m...
size_type nb_convex() const
The total number of convexes in the mesh.
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
Inversion of geometric transformations.
void copy_from(const mesh &m)
Clone a mesh.
Definition: getfem_mesh.cc:446
void exec(size_type nrefine, const mesh_region &cvlst)
build a new mesh_slice.
Dynamic Array.
Definition: dal_basic.h:47
A simple singleton implementation.
void Bank_refine(dal::bit_vector)
Use the Bank strategy to refine some convexes.
size_type points_in_convex(const convex< base_node, TAB > &cv, pgeometric_trans pgt, CONT1 &pftab, CONT2 &itab, bool bruteforce=false)
Search all the points in the convex cv, which is the transformation of the convex cref via the geomet...
"iterator" class for regions.
The output of a getfem::mesh_slicer which has been recorded.
Define various mesh slicers.
GEneric Tool for Finite Element Methods.
number_traits< typename linalg_traits< V1 >::value_type >::magnitude_type vect_dist2_sqr(const V1 &v1, const V2 &v2)
squared Euclidean distance between two vectors
Definition: gmm_blas.h:565
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:239
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:79
size_type add_convex(pconvex_structure cs, ITER ipts, bool *present=0)
Insert a new convex in the mesh_structure.
Define the class getfem::stored_mesh_slice.
pconvex_structure simplex_structure(dim_type nc)
Give a pointer on the structures of a simplex of dimension d.
Mesh structure definition.
mesh_slicer(const mesh &m_)
mesh_slicer constructor.
const ind_cv_ct & ind_points_of_convex(size_type ic) const
Return a container to the list of points attached to convex ic.
short_type nb_faces_of_convex(size_type ic) const
Return the number of faces of convex ic.
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:551
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
mesh & linked_mesh(void) const
Gives a reference to the linked mesh of type mesh.
const dal::bit_vector & convex_index() const
Return the list of valid convex IDs.
Keep informations about a mesh crossed by level-sets.
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
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:780