GetFEM  5.4.2
getfem_mesh_slicers.cc
1 /*===========================================================================
2 
3  Copyright (C) 2004-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 
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  tree.build_tree();
310  }
311 
312  void slicer_mesh_with_mesh::exec(mesh_slicer &ms) {
313  /* identify potientially intersecting convexes of slm */
314  base_node min(ms.nodes[0].pt),max(ms.nodes[0].pt);
315  for (size_type i=1; i < ms.nodes.size(); ++i) {
316  for (size_type k=0; k < min.size(); ++k) {
317  min[k] = std::min(min[k], ms.nodes[i].pt[k]);
318  max[k] = std::max(max[k], ms.nodes[i].pt[k]);
319  }
320  }
321  std::vector<size_type> slmcvs;
322  tree.find_intersecting_boxes(min, max, slmcvs);
323  /* save context */
324  mesh_slicer::cs_simplexes_ct simplexes_final(ms.simplexes);
325  dal::bit_vector splx_in_save(ms.splx_in),
326  simplex_index_save(ms.simplex_index), nodes_index_save(ms.nodes_index);
327  size_type scnt0 = ms.simplexes.size();
328  /* loop over candidate convexes of slm */
329  //cout << "slicer_mesh_with_mesh: convex " << ms.cv << ", " << ms.splx_in.card() << " candidates\n";
330  for (size_type i=0; i < slmcvs.size(); ++i) {
331  size_type slmcv = slmcvs[i];
332  dim_type fcnt_save = dim_type(ms.fcnt);
333  ms.simplexes.resize(scnt0);
334  ms.splx_in = splx_in_save; ms.simplex_index = simplex_index_save; ms.nodes_index = nodes_index_save;
335  //cout << "test intersection of " << ms.cv << " and " << slmcv << "\n";
336  /* loop over the faces and apply slicer_half_space */
337  for (short_type f=0; f<slm.structure_of_convex(slmcv)->nb_faces(); ++f) {
338  base_node x0,n;
339  if (slm.structure_of_convex(slmcv)->dim() == 3 && slm.dim() == 3) {
340  x0 = slm.points_of_face_of_convex(slmcv,f)[0];
341  base_node A = slm.points_of_face_of_convex(slmcv,f)[1] - x0;
342  base_node B = slm.points_of_face_of_convex(slmcv,f)[2] - x0;
343  base_node G = gmm::mean_value(slm.points_of_convex(slmcv).begin(),slm.points_of_convex(slmcv).end());
344  n.resize(3);
345  n[0] = A[1]*B[2] - A[2]*B[1];
346  n[1] = A[2]*B[0] - A[0]*B[2];
347  n[2] = A[0]*B[1] - A[1]*B[0];
348  if (gmm::vect_sp(n,G-x0) > 0) n *= -1.;
349  } else {
350  size_type ip = slm.structure_of_convex(slmcv)->nb_points_of_face(f) / 2;
351  x0 = slm.points_of_face_of_convex(slmcv,f)[ip];
352  n = slm.normal_of_face_of_convex(slmcv,f, x0);
353  }
354  slicer_half_space slf(x0,n,slicer_volume::VOLIN);
355  slf.exec(ms);
356  if (ms.splx_in.card() == 0) break;
357  }
358  if (ms.splx_in.card()) intersection_callback(ms, slmcv);
359  for (dal::bv_visitor is(ms.splx_in); !is.finished(); ++is) {
360  simplexes_final.push_back(ms.simplexes[is]);
361  }
362  ms.fcnt=fcnt_save;
363  }
364  ms.splx_in.clear(); ms.splx_in.add(scnt0, simplexes_final.size()-scnt0);
365  ms.simplexes.swap(simplexes_final);
366  ms.simplex_index = ms.splx_in;
367  ms.update_nodes_index();
368  /*cout << "convex " << ms.cv << "was sliced into " << ms.splx_in.card()
369  << " simplexes, nodes.size=" << ms.nodes.size()
370  << ", used nodes=" << ms.nodes_index.card() << "\n";*/
371  }
372 
373  /* isosurface computations */
374  void slicer_isovalues::prepare(size_type cv,
375  const mesh_slicer::cs_nodes_ct& nodes,
376  const dal::bit_vector& nodes_index) {
377  pt_in.clear(); pt_bin.clear();
378  std::vector<base_node> refpts(nodes.size());
379  Uval.resize(nodes.size());
380  base_vector coeff;
381  base_matrix G;
382  pfem pf = mfU->pmf->fem_of_element(cv);
383  if (pf == 0) return;
384  fem_precomp_pool fprecomp;
385  if (pf->need_G())
386  bgeot::vectors_to_base_matrix
387  (G,mfU->pmf->linked_mesh().points_of_convex(cv));
388  for (size_type i=0; i < nodes.size(); ++i) refpts[i] = nodes[i].pt_ref;
389  pfem_precomp pfp = fprecomp(pf, store_point_tab(refpts));
390  mfU->copy(cv, coeff);
391  //cerr << "cv=" << cv << ", val=" << val << ", coeff=" << coeff << endl;
392  base_vector v(1);
393  fem_interpolation_context ctx(mfU->pmf->linked_mesh().trans_of_convex(cv),
394  pfp, 0, G, cv, short_type(-1));
395  for (dal::bv_visitor i(nodes_index); !i.finished(); ++i) {
396  v[0] = 0;
397  ctx.set_ii(i);
398  pf->interpolation(ctx, coeff, v, mfU->pmf->get_qdim());
399  Uval[i] = v[0];
400  // optimisable -- les bit_vectors sont lents..
401  pt_bin[i] = (gmm::abs(Uval[i] - val) < EPS * val_scaling);
402  pt_in[i] = (Uval[i] - val < 0); if (orient>0) pt_in[i] = !pt_in[i];
403  pt_in[i] = pt_in[i] || pt_bin[i];
404  // cerr << "cv=" << cv << ", node["<< i << "]=" << nodes[i].pt
405  // << ", Uval[i]=" << Uval[i] << ", pt_in[i]=" << pt_in[i]
406  // << ", pt_bin[i]=" << pt_bin[i] << endl;
407  }
408  }
409 
410 
411  void slicer_union::exec(mesh_slicer &ms) {
412  dal::bit_vector splx_in_base = ms.splx_in;
413  size_type c = ms.simplexes.size();
414  dim_type fcnt_0 = dim_type(ms.fcnt);
415  A->exec(ms);
416  dal::bit_vector splx_inA(ms.splx_in);
417  dim_type fcnt_1 = dim_type(ms.fcnt);
418 
419  dal::bit_vector splx_inB = splx_in_base;
420  splx_inB.add(c, ms.simplexes.size()-c);
421  splx_inB.setminus(splx_inA);
422  for (dal::bv_visitor_c i(splx_inB); !i.finished(); ++i) {
423  if (!ms.simplex_index[i]) splx_inB.sup(i);
424  }
425  //splx_inB &= ms.simplex_index;
426  ms.splx_in = splx_inB;
427  B->exec(ms); splx_inB = ms.splx_in;
428  ms.splx_in |= splx_inA;
429 
430  /*
431  the boring part : making sure that the "slice face" node markers
432  are correctly set
433  */
434  for (unsigned f=fcnt_0; f < fcnt_1; ++f) {
435  for (dal::bv_visitor i(splx_inB); !i.finished(); ++i) {
436  for (unsigned j=0; j < ms.simplexes[i].dim()+1; ++j) {
437  bool face_boundA = true;
438  for (unsigned k=0; k < ms.simplexes[i].dim()+1; ++k) {
439  if (j != k && !ms.nodes[ms.simplexes[i].inodes[k]].faces[f]) {
440  face_boundA = false; break;
441  }
442  }
443  if (face_boundA) {
444  /* now we know that the face was on slice A boundary, and
445  that the convex is inside slice B. The conclusion: the
446  face is not on a face of A union B.
447  */
448  for (unsigned k=0; k < ms.simplexes[i].dim()+1; ++k)
449  if (j != k) ms.nodes[ms.simplexes[i].inodes[k]].faces[f] = false;
450  }
451  }
452  }
453  }
454  ms.update_nodes_index();
455  }
456 
457  void slicer_intersect::exec(mesh_slicer& ms) {
458  A->exec(ms);
459  B->exec(ms);
460  }
461 
462  void slicer_complementary::exec(mesh_slicer& ms) {
463  dal::bit_vector splx_inA = ms.splx_in;
464  size_type sz = ms.simplexes.size();
465  A->exec(ms); splx_inA.swap(ms.splx_in);
466  ms.splx_in &= ms.simplex_index;
467  dal::bit_vector bv = ms.splx_in; bv.add(sz, ms.simplexes.size()-sz); bv &= ms.simplex_index;
468  for (dal::bv_visitor_c i(bv); !i.finished(); ++i) {
469  /*cerr << "convex " << cv << ",examining simplex #" << i << ": {";
470  for (size_type j=0; j < simplexes[i].inodes.size(); ++j) cerr << nodes[simplexes[i].inodes[j]].pt << " ";
471  cerr << "}, splx_in=" << splx_in[i] << "splx_inA=" << splx_inA[i] << endl;*/
472  ms.splx_in[i] = !splx_inA.is_in(i);
473  }
474  }
475 
476  void slicer_compute_area::exec(mesh_slicer &ms) {
477  for (dal::bv_visitor is(ms.splx_in); !is.finished(); ++is) {
478  const slice_simplex& s = ms.simplexes[is];
479  base_matrix M(s.dim(),s.dim());
480  for (size_type i=0; i < s.dim(); ++i)
481  for (size_type j=0; j < s.dim(); ++j)
482  M(i,j) = ms.nodes[s.inodes[i+1]].pt[j] - ms.nodes[s.inodes[0]].pt[j];
483  scalar_type v = bgeot::lu_det(&(*(M.begin())), s.dim());
484  for (size_type d=2; d <= s.dim(); ++d) v /= scalar_type(d);
485  a += v;
486  }
487  }
488 
489  void slicer_build_edges_mesh::exec(mesh_slicer &ms) {
490  for (dal::bv_visitor is(ms.splx_in); !is.finished(); ++is) {
491  const slice_simplex& s = ms.simplexes[is];
492  for (size_type i=0; i < s.dim(); ++i) {
493  for (size_type j=i+1; j <= s.dim(); ++j) {
494  const slice_node& A = ms.nodes[s.inodes[i]];
495  const slice_node& B = ms.nodes[s.inodes[j]];
496  /* duplicate with stored_mesh_slice which also
497  builds a list of edges */
498  if ((A.faces & B.faces).count() >= unsigned(ms.cv_dim-1)) {
499  slice_node::faces_ct fmask((1 << ms.cv_nbfaces)-1); fmask.flip();
500  size_type e = edges_m.add_segment_by_points(A.pt,B.pt);
501  if (pslice_edges && (((A.faces & B.faces) & fmask).any())) pslice_edges->add(e);
502  }
503  }
504  }
505  }
506  }
507 
508  void slicer_build_mesh::exec(mesh_slicer &ms) {
509  std::vector<size_type> pid(ms.nodes_index.last_true()+1);
510  for (dal::bv_visitor i(ms.nodes_index); !i.finished(); ++i)
511  pid[i] = m.add_point(ms.nodes[i].pt);
512  for (dal::bv_visitor i(ms.splx_in); !i.finished(); ++i) {
513  for (unsigned j=0; j < ms.simplexes.at(i).inodes.size(); ++j) {
514  assert(m.points_index().is_in(pid.at(ms.simplexes.at(i).inodes[j])));
515  }
516  m.add_convex(bgeot::simplex_geotrans(ms.simplexes[i].dim(),1),
517  gmm::index_ref_iterator(pid.begin(),
518  ms.simplexes[i].inodes.begin()));
519  }
520  }
521 
522  void slicer_explode::exec(mesh_slicer &ms) {
523  if (ms.nodes_index.card() == 0) return;
524 
525  base_node G;
526  if (ms.face < dim_type(-1))
527  G = gmm::mean_value(ms.m.points_of_face_of_convex(ms.cv, ms.face).begin(),
528  ms.m.points_of_face_of_convex(ms.cv, ms.face).end());
529  else
530  G = gmm::mean_value(ms.m.points_of_convex(ms.cv).begin(),
531  ms.m.points_of_convex(ms.cv).end());
532  for (dal::bv_visitor i(ms.nodes_index); !i.finished(); ++i)
533  ms.nodes[i].pt = G + coef*(ms.nodes[i].pt - G);
534 
535  for (dal::bv_visitor cnt(ms.splx_in); !cnt.finished(); ++cnt) {
536  const slice_simplex& s = ms.simplexes[cnt];
537  if (s.dim() == 3) { // keep only faces
538  ms.sup_simplex(cnt);
539  slice_simplex s2(3);
540  for (size_type j=0; j < 4; ++j) {
541  /* usage of s forbidden in this loop since push_back happens .. */
542  static unsigned ord[][3] = {{0,2,1},{1,2,3},{1,3,0},{0,3,2}}; /* keep orientation of faces */
543  for (size_type k=0; k < 3; ++k) { s2.inodes[k] = ms.simplexes[cnt].inodes[ord[j][k]]; } //k + (k<j ? 0 : 1)]; }
544 
545  slice_node::faces_ct f; f.set();
546  for (size_type i=0; i < s2.dim()+1; ++i) {
547  f &= ms.nodes[s2.inodes[i]].faces;
548  }
549  if (f.any()) {
550  ms.add_simplex(s2, true);
551  }
552  }
553  }
554  }
555  }
556 
557  /* -------------------- member functions of mesh_slicer -------------- */
558 
559  mesh_slicer::mesh_slicer(const mesh_level_set &mls_) :
560  m(mls_.linked_mesh()), mls(&mls_), pgt(0), cvr(0) {}
562  m(m_), mls(0), pgt(0), cvr(0) {}
563 
564  void mesh_slicer::using_mesh_level_set(const mesh_level_set &mls_) {
565  mls = &mls_;
566  GMM_ASSERT1(&m == &mls->linked_mesh(), "different meshes");
567  }
568 
569  void mesh_slicer::pack() {
570  std::vector<size_type> new_id(nodes.size());
571  size_type ncnt = 0;
572  for (dal::bv_visitor i(nodes_index); !i.finished(); ++i) {
573  if (i != ncnt) {
574  nodes[i].swap(nodes[ncnt]);
575  }
576  new_id[i] = ncnt++;
577  }
578  nodes.resize(ncnt);
579  size_type scnt = 0;
580  for (dal::bv_visitor j(simplex_index); !j.finished(); ++j) {
581  if (j != scnt) { simplexes[scnt] = simplexes[j]; }
582  for (std::vector<size_type>::iterator it = simplexes[scnt].inodes.begin();
583  it != simplexes[scnt].inodes.end(); ++it) {
584  *it = new_id[*it];
585  }
586  }
587  simplexes.resize(scnt);
588  }
589 
590  void mesh_slicer::update_nodes_index() {
591  nodes_index.clear();
592  for (dal::bv_visitor j(simplex_index); !j.finished(); ++j) {
593  assert(j < simplexes.size());
594  for (std::vector<size_type>::iterator it = simplexes[j].inodes.begin();
595  it != simplexes[j].inodes.end(); ++it) {
596  assert(*it < nodes.size());
597  nodes_index.add(*it);
598  }
599  }
600  }
601 
602  static void flag_points_on_faces(const bgeot::pconvex_ref& cvr,
603  const std::vector<base_node>& pts,
604  std::vector<slice_node::faces_ct>& faces) {
605  GMM_ASSERT1(cvr->structure()->nb_faces() <= 32,
606  "won't work for convexes with more than 32 faces "
607  "(hardcoded limit)");
608  faces.resize(pts.size());
609  for (size_type i=0; i < pts.size(); ++i) {
610  faces[i].reset();
611  for (short_type f=0; f < cvr->structure()->nb_faces(); ++f)
612  faces[i][f] = (gmm::abs(cvr->is_in_face(f, pts[i])) < 1e-10);
613  }
614  }
615 
616  void mesh_slicer::update_cv_data(size_type cv_, short_type f_) {
617  cv = cv_;
618  face = f_;
619  pgt = m.trans_of_convex(cv);
620  prev_cvr = cvr; cvr = pgt->convex_ref();
621  cv_dim = cvr->structure()->dim();
622  cv_nbfaces = cvr->structure()->nb_faces();
623  fcnt = cvr->structure()->nb_faces();
624  discont = (mls && mls->is_convex_cut(cv));
625  }
626 
627  void mesh_slicer::apply_slicers() {
628  simplex_index.clear(); simplex_index.add(0, simplexes.size());
629  splx_in = simplex_index;
630  nodes_index.clear(); nodes_index.add(0, nodes.size());
631  for (size_type i=0; i < action.size(); ++i) {
632  action[i]->exec(*this);
633  //cout << "simplex_index=" << simplex_index << "\n splx_in=" << splx_in << "\n";
634  assert(simplex_index.contains(splx_in));
635  }
636  }
637 
638  void mesh_slicer::simplex_orientation(slice_simplex& s) {
639  size_type N = m.dim();
640  if (s.dim() == N) {
641  base_matrix M(N,N);
642  for (size_type i=1; i < N+1; ++i) {
643  base_small_vector d = nodes[s.inodes[i]].pt - nodes[s.inodes[0]].pt;
644  gmm::copy_n(d.const_begin(), N, M.begin() + (i-1)*N);
645  }
646  scalar_type J = bgeot::lu_det(&(*(M.begin())), N);
647  //cout << " lu_det = " << J << "\n";
648  if (J < 0) {
649  std::swap(s.inodes[1],s.inodes[0]);
650  }
651  }
652  }
653 
654  void mesh_slicer::exec(size_type nrefine, const mesh_region& cvlst) {
655  short_type n = short_type(nrefine);
656  exec_(&n, 0, cvlst);
657  }
658 
659  void mesh_slicer::exec(const std::vector<short_type> &nrefine,
660  const mesh_region& cvlst) {
661  exec_(&nrefine[0], 1, cvlst);
662  }
663 
664  static bool check_orient(size_type cv, bgeot::pgeometric_trans pgt,
665  const mesh& m) {
666  if (pgt->dim() == m.dim() && m.dim()>=2) { /* no orient check for
667  convexes of lower dim */
668  base_matrix G; bgeot::vectors_to_base_matrix(G,m.points_of_convex(cv));
669  base_node g(pgt->dim()); g.fill(.5);
670  base_matrix pc; pgt->poly_vector_grad(g,pc);
671  base_matrix K(pgt->dim(),pgt->dim());
672  gmm::mult(G,pc,K);
673  scalar_type J = bgeot::lu_det(&(*(K.begin())), pgt->dim());
674  // bgeot::geotrans_interpolation_context ctx(pgp,0,G);
675  // scalar_type J = gmm::lu_det(ctx.B()); // pb car inverse K mÍme
676  if (J < 0) return true;
677  //cout << "cv = " << cv << ", J = " << J << "\n";
678  }
679  return false;
680  }
681 
682 #if OLD_MESH_SLICER
683  void mesh_slicer::exec_(const short_type *pnrefine, int nref_stride,
684  const mesh_region& cvlst) {
685  std::vector<base_node> cvm_pts;
686  const bgeot::basic_mesh *cvm = 0;
687  const bgeot::mesh_structure *cvms = 0;
689  bgeot::pgeotrans_precomp pgp = 0;
690  std::vector<slice_node::faces_ct> points_on_faces;
691 
692  cvlst.from_mesh(m);
693  size_type prev_nrefine = 0;
694  for (mr_visitor it(cvlst); !it.finished(); ++it) {
695  size_type nrefine = pnrefine[it.cv()*nref_stride];
696  update_cv_data(it.cv(),it.f());
697  bool revert_orientation = check_orient(cv, pgt,m);
698 
699  /* update structure-dependent data */
700  if (prev_cvr != cvr || nrefine != prev_nrefine) {
701  cvm = bgeot::refined_simplex_mesh_for_convex(cvr, nrefine);
702  cvm_pts.resize(cvm->points().card());
703  std::copy(cvm->points().begin(), cvm->points().end(), cvm_pts.begin());
704  pgp = gppool(pgt, store_point_tab(cvm_pts));
705  flag_points_on_faces(cvr, cvm_pts, points_on_faces);
706  prev_nrefine = nrefine;
707  }
708  if (face < dim_type(-1))
709  cvms = bgeot::refined_simplex_mesh_for_convex_faces(cvr, nrefine)[face].get();
710  else
711  cvms = cvm;
712 
713  /* apply the initial geometric transformation */
714  std::vector<size_type> ptsid(cvm_pts.size()); std::fill(ptsid.begin(), ptsid.end(), size_type(-1));
715  simplexes.resize(cvms->nb_convex());
716  nodes.resize(0);
717  for (size_type snum = 0; snum < cvms->nb_convex(); ++snum) {
718  /* cvms should not contain holes in its convex index.. */
719  simplexes[snum].inodes.resize(cvms->nb_points_of_convex(snum));
720  std::copy(cvms->ind_points_of_convex(snum).begin(),
721  cvms->ind_points_of_convex(snum).end(), simplexes[snum].inodes.begin());
722  if (revert_orientation) std::swap(simplexes[snum].inodes[0],simplexes[snum].inodes[1]);
723  /* store indices of points which are really used , and renumbers them */
724  for (std::vector<size_type>::iterator itp = simplexes[snum].inodes.begin();
725  itp != simplexes[snum].inodes.end(); ++itp) {
726  if (ptsid[*itp] == size_type(-1)) {
727  ptsid[*itp] = nodes.size();
728  nodes.push_back(slice_node());
729  nodes.back().pt_ref = cvm_pts[*itp];
730  nodes.back().faces = points_on_faces[*itp];
731  nodes.back().pt.resize(m.dim()); nodes.back().pt.fill(0.);
732  pgp->transform(m.points_of_convex(cv), *itp, nodes.back().pt);
733  }
734  *itp = ptsid[*itp];
735  }
736  if (0) {
737  static int once = 0;
738  if (once++ < 3) {
739  cout << "check orient cv " << cv << ", snum=" << snum << "/" << cvms->nb_convex();
740  }
741  simplex_orientation(simplexes[snum]);
742  }
743  }
744  apply_slicers();
745  }
746  }
747 #endif // OLD_MESH_SLICER
748 
749  template <typename CONT>
750  static void add_degree1_convex(bgeot::pgeometric_trans pgt, const CONT &pts,
751  mesh &m) {
752  unsigned N = pgt->dim();
753  std::vector<base_node> v; v.reserve(N+1);
754  for (unsigned i=0; i < pgt->nb_points(); ++i) {
755  const base_node P = pgt->convex_ref()->points()[i];
756  scalar_type s = 0;
757  for (unsigned j=0; j < N; ++j) {
758  s += P[j]; if (P[j] == 1) { v.push_back(pts[i]); break; }
759  }
760  if (s == 0) v.push_back(pts[i]);
761  }
762  assert(v.size() == N+1);
763  base_node G = gmm::mean_value(v);
764  /*for (unsigned i=0; i < v.size();++i)
765  v[i] = v[i] + 0.1 * (G - v[i]);*/
766  m.add_convex_by_points(bgeot::simplex_geotrans(N,1), v.begin());
767  }
768 
769  const mesh& mesh_slicer::refined_simplex_mesh_for_convex_cut_by_level_set
770  (const mesh &cvm, unsigned nrefine) {
771  mesh mm; mm.copy_from(cvm);
772  while (nrefine > 1) {
773  mm.Bank_refine(mm.convex_index());
774  nrefine /= 2;
775  }
776 
777  std::vector<size_type> idx;
778  tmp_mesh.clear();
779  //cerr << "nb cv = " << tmp_mesh.convex_index().card() << "\n";
780  for (dal::bv_visitor_c ic(mm.convex_index()); !ic.finished(); ++ic) {
781  add_degree1_convex(mm.trans_of_convex(ic), mm.points_of_convex(ic), tmp_mesh);
782  }
783  /*tmp_mesh.write_to_file(std::cerr);
784  cerr << "\n";*/
785  return tmp_mesh;
786  }
787 
788  const bgeot::mesh_structure &
789  mesh_slicer::refined_simplex_mesh_for_convex_faces_cut_by_level_set
790  (short_type f) {
791  mesh &cvm = tmp_mesh;
792  tmp_mesh_struct.clear();
793  unsigned N = cvm.dim();
794 
795  dal::bit_vector pt_in_face; pt_in_face.sup(0, cvm.points().index().last_true()+1);
796  for (dal::bv_visitor ip(cvm.points().index()); !ip.finished(); ++ip)
797  if (gmm::abs(cvr->is_in_face(short_type(f), cvm.points()[ip]))) pt_in_face.add(ip);
798 
799  for (dal::bv_visitor_c ic(cvm.convex_index()); !ic.finished(); ++ic) {
800  for (short_type ff=0; ff < cvm.nb_faces_of_convex(ic); ++ff) {
801  bool face_ok = true;
802  for (unsigned i=0; i < N; ++i) {
803  if (!pt_in_face.is_in(cvm.ind_points_of_face_of_convex(ic,ff)[i])) {
804  face_ok = false; break;
805  }
806  }
807  if (face_ok) {
808  tmp_mesh_struct.add_convex(bgeot::simplex_structure(dim_type(N-1)),
809  cvm.ind_points_of_face_of_convex(ic, ff).begin());
810  }
811  }
812  }
813  return tmp_mesh_struct;
814  }
815 
816  void mesh_slicer::exec_(const short_type *pnrefine,
817  int nref_stride,
818  const mesh_region& cvlst) {
819  std::vector<base_node> cvm_pts;
820  const bgeot::basic_mesh *cvm = 0;
821  const bgeot::mesh_structure *cvms = 0;
823  bgeot::pgeotrans_precomp pgp = 0;
824  std::vector<slice_node::faces_ct> points_on_faces;
825  bool prev_discont = true;
826 
827  cvlst.from_mesh(m);
828  size_type prev_nrefine = 0;
829  // size_type prev_cv = size_type(-1);
830  for (mr_visitor it(cvlst); !it.finished(); ++it) {
831  size_type nrefine = pnrefine[it.cv()*nref_stride];
832  update_cv_data(it.cv(),it.f());
833  bool revert_orientation = check_orient(cv, pgt,m);
834 
835  /* update structure-dependent data */
836  /* TODO : fix levelset handling when slicing faces .. */
837  if (prev_cvr != cvr || nrefine != prev_nrefine
838  || discont || prev_discont) {
839  if (discont) {
840  cvm = &refined_simplex_mesh_for_convex_cut_by_level_set
841  (mls->mesh_of_convex(cv), unsigned(nrefine));
842  } else {
844  }
845  cvm_pts.resize(cvm->points().card());
846  std::copy(cvm->points().begin(), cvm->points().end(), cvm_pts.begin());
847  pgp = gppool(pgt, store_point_tab(cvm_pts));
848  flag_points_on_faces(cvr, cvm_pts, points_on_faces);
849  prev_nrefine = nrefine;
850  }
851  if (face < dim_type(-1)) {
852  if (!discont) {
854  (cvr, short_type(nrefine))[face].get();
855  } else {
856  cvms = &refined_simplex_mesh_for_convex_faces_cut_by_level_set(face);
857  }
858  } else {
859  cvms = cvm;
860  }
861 
862  /* apply the initial geometric transformation */
863  std::vector<size_type> ptsid(cvm_pts.size()); std::fill(ptsid.begin(), ptsid.end(), size_type(-1));
864  simplexes.resize(cvms->nb_convex());
865  nodes.resize(0);
866 
867  base_node G;
868  for (size_type snum = 0; snum < cvms->nb_convex(); ++snum) {
869  /* cvms should not contain holes in its convex index.. */
870  simplexes[snum].inodes.resize(cvms->nb_points_of_convex(snum));
871  std::copy(cvms->ind_points_of_convex(snum).begin(),
872  cvms->ind_points_of_convex(snum).end(), simplexes[snum].inodes.begin());
873  if (revert_orientation) std::swap(simplexes[snum].inodes[0],simplexes[snum].inodes[1]);
874  /* store indices of points which are really used , and renumbers them */
875  if (discont) {
876  G.resize(m.dim()); G.fill(0.);
877  for (std::vector<size_type>::iterator itp =
878  simplexes[snum].inodes.begin();
879  itp != simplexes[snum].inodes.end(); ++itp) {
880  G += cvm_pts[*itp];
881  }
882  G /= scalar_type(simplexes[snum].inodes.size());
883  }
884 
885  for (std::vector<size_type>::iterator itp =
886  simplexes[snum].inodes.begin();
887  itp != simplexes[snum].inodes.end(); ++itp) {
888  if (discont || ptsid[*itp] == size_type(-1)) {
889  ptsid[*itp] = nodes.size();
890  nodes.push_back(slice_node());
891  if (!discont) {
892  nodes.back().pt_ref = cvm_pts[*itp];
893  } else {
894  /* displace the ref point such that one will not interpolate
895  on the discontinuity (yes this is quite ugly and not
896  robust)
897  */
898  nodes.back().pt_ref = cvm_pts[*itp] + 0.01*(G - cvm_pts[*itp]);
899  }
900  nodes.back().faces = points_on_faces[*itp];
901  nodes.back().pt.resize(m.dim()); nodes.back().pt.fill(0.);
902  pgp->transform(m.points_of_convex(cv), *itp, nodes.back().pt);
903  //nodes.back().pt = pgt->transform(G, m.points_of_convex(cv));
904  //cerr << "G = " << G << " -> pt = " << nodes.back().pt << "\n";
905  }
906  *itp = ptsid[*itp];
907  }
908  }
909  //cerr << "cv = " << cv << ", cvm.nb_points_ = "<< cvm->points().size() << ", nbnodes = " << nodes.size() << ", nb_simpl=" << simplexes.size() << "\n";
910 
911  apply_slicers();
912  // prev_cv = it.cv();
913  prev_discont = discont;
914  }
915  }
916 
917 
918  void mesh_slicer::exec(size_type nrefine) {
919  exec(nrefine,mesh_region(m.convex_index()));
920  }
921 
922  /* apply slice ops to an already stored slice object */
924  GMM_ASSERT1(&sl.linked_mesh() == &m, "wrong mesh");
925  for (stored_mesh_slice::cvlst_ct::const_iterator it = sl.cvlst.begin(); it != sl.cvlst.end(); ++it) {
926  update_cv_data((*it).cv_num);
927  nodes = (*it).nodes;
928  simplexes = (*it).simplexes;
929  apply_slicers();
930  }
931  }
932 
933  /* apply slice ops to a set of nodes */
934  void mesh_slicer::exec(const std::vector<base_node>& pts) {
936  gti.add_points(pts);
939  for (dal::bv_visitor ic(m.convex_index()); !ic.finished(); ++ic) {
940  size_type nb = gti.points_in_convex(m.convex(ic), m.trans_of_convex(ic), ptab, itab);
941  if (!nb) continue;
942  update_cv_data(ic);
943  nodes.resize(0); simplexes.resize(0);
944  for (size_type i=0; i < nb; ++i) {
945  //cerr << "point " << itab[i] << "(" << pts[itab[i]]
946  //<< ") trouve dans le convex " << ic << " [pt_ref=" << ptab[i] << "]\n";
947  nodes.push_back(slice_node(pts[itab[i]],ptab[i])); nodes.back().faces=0;
948  slice_simplex s(1); s.inodes[0] = nodes.size()-1;
949  simplexes.push_back(s);
950  }
951  apply_slicers();
952  }
953  }
954 }
dal::singleton::instance
static T & instance()
Instance from the current thread.
Definition: dal_singleton.h:165
bgeot::geotrans_inv::points_in_convex
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...
Definition: bgeot_geotrans_inv.h:242
getfem::mesh::clear
void clear()
Erase the mesh.
Definition: getfem_mesh.cc:273
bgeot::geotrans_inv
handles the geometric inversion for a given (supposedly quite large) set of points
Definition: bgeot_geotrans_inv.h:178
getfem::stored_mesh_slice
The output of a getfem::mesh_slicer which has been recorded.
Definition: getfem_mesh_slice.h:47
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
dal_singleton.h
A simple singleton implementation.
getfem_mesh_level_set.h
Keep informations about a mesh crossed by level-sets.
getfem::mesh::add_segment_by_points
size_type add_segment_by_points(const base_node &pt1, const base_node &pt2)
Add a segment to the mesh, given the coordinates of its vertices.
Definition: getfem_mesh.h:262
bgeot::geotrans_precomp_pool
The object geotrans_precomp_pool Allow to allocate a certain number of geotrans_precomp and automatic...
Definition: bgeot_geometric_trans.h:383
getfem::mesh::normal_of_face_of_convex
base_small_vector normal_of_face_of_convex(size_type ic, short_type f, const base_node &pt) const
Return the normal of the given convex face, evaluated at the point pt.
Definition: getfem_mesh.cc:384
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::mesh_structure::clear
void clear()
erase the mesh
Definition: bgeot_mesh_structure.cc:205
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
getfem::mesh_slicer::mesh_slicer
mesh_slicer(const mesh &m_)
mesh_slicer constructor.
Definition: getfem_mesh_slicers.cc:561
getfem::mesh_level_set
Keep informations about a mesh crossed by level-sets.
Definition: getfem_mesh_level_set.h:52
getfem::slicer_volume::prepare
virtual void prepare(size_type, const mesh_slicer::cs_nodes_ct &nodes, const dal::bit_vector &nodes_index)
Overload either 'prepare' or 'test_point'.
Definition: getfem_mesh_slicers.h:307
getfem::mesh::points_of_face_of_convex
ref_mesh_face_pt_ct points_of_face_of_convex(size_type ic, short_type f) const
Return a (pseudo)container of points of face of a given convex.
Definition: getfem_mesh.h:182
dal::dynamic_array
Dynamic Array.
Definition: dal_basic.h:47
getfem_mesh_slice.h
Define the class getfem::stored_mesh_slice.
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
bgeot::mesh_structure::nb_convex
size_type nb_convex() const
The total number of convexes in the mesh.
Definition: bgeot_mesh_structure.h:96
bgeot::geotrans_inv::add_points
void add_points(const CONT &c)
Add the points contained in c to the list of points.
Definition: bgeot_geotrans_inv.h:187
bgeot_geotrans_inv.h
Inversion of geometric transformations.
bgeot::refined_simplex_mesh_for_convex_faces
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
Definition: bgeot_poly_composite.cc:565
getfem::slicer_volume::orient
int orient
orient defines the kind of slicing : VOLIN -> keep the inside of the volume, VOLBOUND -> its boundary...
Definition: getfem_mesh_slicers.h:302
getfem::slicer_volume::edge_intersect
virtual scalar_type edge_intersect(size_type, size_type, const mesh_slicer::cs_nodes_ct &) const =0
edge_intersect should always be overloaded
bgeot::refined_simplex_mesh_for_convex
const basic_mesh * refined_simplex_mesh_for_convex(pconvex_ref cvr, short_type k)
simplexify a convex_ref.
Definition: bgeot_poly_composite.cc:558
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:47
bgeot::mesh_structure
Mesh structure definition.
Definition: bgeot_mesh_structure.h:73
getfem_mesh_slicers.h
Define various mesh slicers.
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
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
bgeot::mesh_structure::add_convex
size_type add_convex(pconvex_structure cs, ITER ipts, bool *present=0)
Insert a new convex in the mesh_structure.
Definition: bgeot_mesh_structure.h:306
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
getfem::stored_mesh_slice::linked_mesh
const mesh & linked_mesh() const
return a pointer to the original mesh
Definition: getfem_mesh_slice.h:110
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::mesh_slicer::exec
void exec(size_type nrefine, const mesh_region &cvlst)
build a new mesh_slice.
Definition: getfem_mesh_slicers.cc:654
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
bgeot::mesh_structure::structure_of_convex
pconvex_structure structure_of_convex(size_type ic) const
Return the pconvex_structure of the convex ic.
Definition: bgeot_mesh_structure.h:112