GetFEM  5.4.2
getfem_mesh_im_level_set.cc
1 /*===========================================================================
2 
3  Copyright (C) 2005-2020 Yves Renard
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/getfem_mesh_im_level_set.h"
23 #include "getfem/getfem_mesher.h"
24 #include "getfem/bgeot_kdtree.h"
25 
26 namespace getfem {
27 
28  void mesh_im_level_set::update_from_context(void) const
29  { is_adapted = false; }
30 
31  void mesh_im_level_set::clear_build_methods() {
32  for (size_type i = 0; i < build_methods.size(); ++i)
33  del_stored_object(build_methods[i]);
34  build_methods.clear();
35  cut_im.clear();
36  }
37 
38  void mesh_im_level_set::clear(void) {
39  mesh_im::clear();
40  clear_build_methods();
41  is_adapted = false;
42  }
43 
44  void mesh_im_level_set::init_with_mls(mesh_level_set &me,
45  int integrate_where_,
46  pintegration_method reg,
47  pintegration_method sing) {
48  init_with_mesh(me.linked_mesh());
49  cut_im.init_with_mesh(me.linked_mesh());
50  mls = &me;
51  integrate_where = integrate_where_;
52  set_simplex_im(reg, sing);
53  this->add_dependency(*mls);
54  is_adapted = false;
55  }
56 
57  mesh_im_level_set::mesh_im_level_set(mesh_level_set &me,
58  int integrate_where_,
59  pintegration_method reg,
60  pintegration_method sing) {
61  mls = 0;
62  init_with_mls(me, integrate_where_, reg, sing);
63  }
64 
65  mesh_im_level_set::mesh_im_level_set(void)
66  { mls = 0; is_adapted = false; }
67 
68 
69  pintegration_method
70  mesh_im_level_set::int_method_of_element(size_type cv) const {
71  if (!is_adapted) const_cast<mesh_im_level_set *>(this)->adapt();
72  if (cut_im.convex_index().is_in(cv))
73  return cut_im.int_method_of_element(cv);
74  else {
75  if (ignored_im.is_in(cv)) //integrate_where == INTEGRATE_BOUNDARY)
76  return getfem::im_none();
77 
78  return mesh_im::int_method_of_element(cv);
79  }
80  }
81 
82  DAL_SIMPLE_KEY(special_imls_key, papprox_integration);
83 
84  /* only for INTEGRATE_INSIDE or INTEGRATE_OUTSIDE */
85  mesh_im_level_set::bool2 mesh_im_level_set::is_point_in_selected_area2
86  (const std::vector<pmesher_signed_distance> &mesherls0,
87  const std::vector<pmesher_signed_distance> &mesherls1,
88  const base_node& P) {
89  bool isin = true;
90  int isbin = 0;
91  for (unsigned i = 0; i < mls->nb_level_sets(); ++i) {
92  isin = isin && ((*(mesherls0[i]))(P) < 0);
93  if (gmm::abs((*(mesherls0[i]))(P)) < 1e-7)
94  isbin = i+1;
95  if (mls->get_level_set(i)->has_secondary())
96  isin = isin && ((*(mesherls1[i]))(P) < 0);
97  }
98  bool2 b;
99  b.in = ((integrate_where & INTEGRATE_OUTSIDE)) ? !isin : isin;
100  b.bin = isbin;
101  return b;
102  }
103 
104 
105  /* very rustic set operations evaluator */
106  struct is_in_eval {
107  dal::bit_vector in; // levelsets for which the point is "inside"
108  dal::bit_vector bin; // levelsets for which the point is on the boundary
109  typedef mesh_im_level_set::bool2 bool2;
110  bool2 do_expr(const char *&s) {
111  bool2 r;
112  if (*s == '(') {
113  r = do_expr(++s);
114  GMM_ASSERT1(*s++ == ')',
115  "expecting ')' in csg expression at '" << s-1 << "'");
116  } else if (*s == '!') { // complementary
117  r = do_expr(++s); r.in = !r.in;
118  } else if (*s >= 'a' && *s <= 'z') {
119  unsigned idx = (*s) - 'a';
120  r.in = in.is_in(idx);
121  r.bin = bin.is_in(idx) ? idx+1 : 0;
122  ++s;
123  } else
124  GMM_ASSERT1(false, "parse error in csg expression at '" << s << "'");
125  if (*s == '+') { // Union
126  //cerr << "s = " << s << ", r = " << r << "\n";
127  bool2 a = r, b = do_expr(++s);
128  //cerr << "->b = " << b << "\n";
129  r.in = b.in || a.in;
130  if (b.bin && !a.in) r.bin = b.bin;
131  else if (a.bin && !b.in) r.bin = a.bin;
132  else r.bin = 0;
133  } else if (*s == '-') { // Set difference
134  bool2 a = r, b = do_expr(++s);
135  r.in = a.in && !b.in;
136  if (a.bin && !b.in) r.bin = a.bin;
137  else if (a.in && b.bin) r.bin = b.bin;
138  else r.bin = 0;
139  } else if (*s == '*') { // Intersection
140  bool2 a = r, b = do_expr(++s);
141  r.in = a.in && b.in;
142  if (a.bin && b.in) r.bin = a.bin;
143  else if (a.in && b.bin) r.bin = b.bin;
144  else r.bin = 0;
145  }
146  return r;
147  }
148  bool2 is_in(const char*s) {
149  bool2 b = do_expr(s);
150  GMM_ASSERT1(!(*s), "parse error in CSG expression at " << s);
151  return b;
152  }
153  void check() {
154  const char *s[] = { "a*b*c*d",
155  "a+b+c+d",
156  "(a+b+c+d)",
157  "d*(a+b+c)",
158  "(a+b)-(c+d)",
159  "((a+b)-(c+d))",
160  "!a",
161  0 };
162  for (const char **p = s; *p; ++p)
163  cerr << *p << "\n";
164  for (unsigned c=0; c < 16; ++c) {
165  in[0] = (c&1); bin[0] = 1;
166  in[1] = (c&2); bin[1] = 1;
167  in[2] = (c&4); bin[2] = 1;
168  in[3] = (c&8); bin[3] = 1;
169  cerr << in[0] << in[1] << in[2] << in[3] << ": ";
170  for (const char **p = s; *p; ++p) {
171  bool2 b = is_in(*p);
172  cerr << b.in << "/" << b.bin << " ";
173  }
174  cerr << "\n";
175  }
176  }
177  };
178 
179  mesh_im_level_set::bool2
180  mesh_im_level_set::is_point_in_selected_area
181  (const std::vector<pmesher_signed_distance> &mesherls0,
182  const std::vector<pmesher_signed_distance> &mesherls1,
183  const base_node& P) {
184  is_in_eval ev;
185  for (unsigned i = 0; i < mls->nb_level_sets(); ++i) {
186  bool sec = mls->get_level_set(i)->has_secondary();
187  scalar_type d1 = (*(mesherls0[i]))(P);
188  scalar_type d2 = (sec ? (*(mesherls1[i]))(P) : -1);
189  if (d1 < 0 && d2 < 0) ev.in.add(i);
190  // if ((integrate_where & INTEGRATE_OUTSIDE) /*&& !sec*/)
191  // ev.in[i].flip();
192 
193  if (gmm::abs(d1) < 1e-7 && d2 < 1e-7)
194  ev.bin.add(i);
195  }
196 
197 
198  bool2 r;
199  if (ls_csg_description.size())
200  r = ev.is_in(ls_csg_description.c_str());
201  else {
202  r.in = (ev.in.card() == mls->nb_level_sets());
203  r.bin = (ev.bin.card() >= 1 && ev.in.card() >= mls->nb_level_sets()-1);
204  }
205 
206  if (integrate_where & INTEGRATE_OUTSIDE) r.in = !(r.in);
207 
208 
209 
210  /*bool2 r2 = is_point_in_selected_area2(mesherls0,mesherls1,P);
211  if (r2.in != r.in || r2.bin != r.bin) {
212  cerr << "ev.in = " << ev.in << ", bin=" << ev.bin<<"\n";
213  cerr << "is_point_in_selected_area2("<<P <<"): r="<<r.in<<"/"<<r.bin
214  << ", r2=" << r2.in<<"/"<<r2.bin <<"\n";
215  assert(0);
216  }*/
217 
218  return r;
219  }
220 
221  void mesh_im_level_set::build_method_of_convex(size_type cv) {
222  const mesh &msh(mls->mesh_of_convex(cv));
223  GMM_ASSERT3(msh.convex_index().card() != 0, "Internal error");
224  base_matrix G;
225  base_node B;
226 
227  std::vector<pmesher_signed_distance> mesherls0(mls->nb_level_sets());
228  std::vector<pmesher_signed_distance> mesherls1(mls->nb_level_sets());
229  dal::bit_vector convexes_arein;
230 
231  //std::fstream totof("totof", std::ios::out | std::ios::app);
232  for (unsigned i = 0; i < mls->nb_level_sets(); ++i) {
233  mesherls0[i] = mls->get_level_set(i)->mls_of_convex(cv, 0, false);
234  if (mls->get_level_set(i)->has_secondary())
235  mesherls1[i] = mls->get_level_set(i)->mls_of_convex(cv, 1, false);
236  }
237 
238  if (integrate_where != (INTEGRATE_ALL)) {
239  for (dal::bv_visitor scv(msh.convex_index()); !scv.finished(); ++scv) {
240  B = gmm::mean_value(msh.points_of_convex(scv));
241  convexes_arein[scv] =
242  is_point_in_selected_area(mesherls0, mesherls1, B).in;
243  }
244  }
245 
246  bgeot::pgeometric_trans pgt = linked_mesh().trans_of_convex(cv);
247  bgeot::pgeometric_trans pgt2
248  = msh.trans_of_convex(msh.convex_index().first_true());
249  dim_type n = pgt->dim();
250 
251  if (base_singular_pim) GMM_ASSERT1
252  ((n != 2 ||
253  base_singular_pim->structure()== bgeot::parallelepiped_structure(2))
254  && (n != 3
255  || base_singular_pim->structure() == bgeot::prism_P1_structure(3))
256  && (n >= 2) && (n <= 3),
257  "Base integration method for quasi polar integration not convenient");
258 
259  auto new_approx = std::make_shared<approx_integration>(pgt->convex_ref());
260  new_approx->set_built_on_the_fly();
261  base_matrix KK(n,n), CS(n,n);
262  base_matrix pc(pgt2->nb_points(), n);
263  std::vector<size_type> ptsing;
264 
265  // cout << "testing convex " << cv << ", " << msh.convex_index().card() << " subconvexes\n";
266 
267  for (dal::bv_visitor i(msh.convex_index()); !i.finished(); ++i) {
268  papprox_integration pai = regular_simplex_pim->approx_method();
269 
270  GMM_ASSERT1(regular_simplex_pim->structure() == bgeot::simplex_structure(n), "Base integration method should be defined on a simplex of same dimension than the mesh");
271 
272  if ((integrate_where != INTEGRATE_ALL) &&
273  !convexes_arein[i]) continue;
274 
275  if (base_singular_pim && mls->crack_tip_convexes().is_in(cv)) {
276  ptsing.resize(0);
277  unsigned sing_ls = unsigned(-1);
278 
279  for (unsigned ils = 0; ils < mls->nb_level_sets(); ++ils)
280  if (mls->get_level_set(ils)->has_secondary()) {
281  for (unsigned ipt = 0; ipt <= n; ++ipt) {
282  if (gmm::abs((*(mesherls0[ils]))(msh.points_of_convex(i)[ipt]))
283  < 1E-10
284  && gmm::abs((*(mesherls1[ils]))(msh.points_of_convex(i)[ipt]))
285  < 1E-10) {
286  if (sing_ls == unsigned(-1)) sing_ls = ils;
287  GMM_ASSERT1(sing_ls == ils, "Two singular point in one "
288  "sub element : " << sing_ls << ", " << ils <<
289  ". To be done.");
290  ptsing.push_back(ipt);
291  }
292  }
293  }
294  assert(ptsing.size() < n);
295 
296  if (ptsing.size() > 0) {
297  std::stringstream sts;
298  sts << "IM_QUASI_POLAR(" << name_of_int_method(base_singular_pim)
299  << ", " << ptsing[0];
300  if (ptsing.size() > 1) sts << ", " << ptsing[1];
301  sts << ")";
302  pai = int_method_descriptor(sts.str())->approx_method();
303  }
304  }
305 
306  base_matrix G2;
307  vectors_to_base_matrix(G2, linked_mesh().points_of_convex(cv));
308  bgeot::geotrans_interpolation_context
309  cc(linked_mesh().trans_of_convex(cv), pai->point(0), G2);
310 
311  if (integrate_where & (INTEGRATE_INSIDE | INTEGRATE_OUTSIDE)) {
312 
313  vectors_to_base_matrix(G, msh.points_of_convex(i));
314  bgeot::geotrans_interpolation_context c(msh.trans_of_convex(i),
315  pai->point(0), G);
316 
317  for (size_type j = 0; j < pai->nb_points_on_convex(); ++j) {
318  c.set_xref(pai->point(j));
319  pgt2->poly_vector_grad(pai->point(j), pc);
320  gmm::mult(G,pc,KK);
321  scalar_type J = gmm::lu_det(KK);
322  new_approx->add_point(c.xreal(), pai->coeff(j) * gmm::abs(J));
323 
324  /*if (integrate_where == INTEGRATE_INSIDE) {
325  cc.set_xref(c.xreal());
326  totof << cc.xreal()[0] << "\t" << cc.xreal()[1] << "\t"
327  << pai->coeff(j) * gmm::abs(J) << "\n";
328  }*/
329  }
330  }
331 
332  // pgt2 = msh.trans_of_convex(i);
333 
334  for (short_type f = 0; f < pgt2->structure()->nb_faces(); ++f) {
335  short_type ff = short_type(-1);
336  unsigned isin = unsigned(-1);
337 
338  if (integrate_where == INTEGRATE_BOUNDARY) {
339  bool lisin = true;
340  for (unsigned ipt = 0; ipt <
341  pgt2->structure()->nb_points_of_face(f); ++ipt) {
342  const base_node &P = msh.points_of_face_of_convex(i, f)[ipt];
343  isin = is_point_in_selected_area(mesherls0, mesherls1, P).bin;
344  //cerr << P << ":" << isin << " ";
345  if (!isin) { lisin = false; break; }
346  }
347  if (!lisin) continue;
348  else isin--;
349  } else {
350  B = gmm::mean_value(msh.points_of_face_of_convex(i, f));
351  if (pgt->convex_ref()->is_in(B) < -1E-7) continue;
352  for (short_type fi = 0; fi < pgt->structure()->nb_faces(); ++fi) {
353  if (gmm::abs(pgt->convex_ref()->is_in_face(fi, B)) < 2E-6) ff = fi;
354  }
355 
356  if (ff == short_type(-1)) {
357  cout << "Distance to the element : "
358  << pgt->convex_ref()->is_in(B) << endl;
359  for (short_type fi = 0; fi < pgt->structure()->nb_faces(); ++fi) {
360  cout << "Distance to face " << fi << " : "
361  << gmm::abs(pgt->convex_ref()->is_in_face(fi, B)) << endl;
362  }
363  GMM_ASSERT3(false, "the point is neither in the interior nor "
364  "the boundary of the element");
365  }
366  }
367 
368  vectors_to_base_matrix(G, msh.points_of_convex(i));
369  bgeot::geotrans_interpolation_context c(msh.trans_of_convex(i),
370  pai->point(0), G);
371 
372 
373  for (size_type j = 0; j < pai->nb_points_on_face(f); ++j) {
374  if (gmm::abs(c.J()) > 1E-11) {
375  c.set_xref(pai->point_on_face(f, j));
376  base_small_vector un = pgt2->normals()[f], up(msh.dim());
377  gmm::mult(c.B(), un, up);
378  scalar_type nup = gmm::vect_norm2(up);
379 
380  scalar_type nnup(1);
381  if (integrate_where == INTEGRATE_BOUNDARY) {
382  cc.set_xref(c.xreal());
383  mesherls0[isin]->grad(c.xreal(), un);
384  un /= gmm::vect_norm2(un);
385  gmm::mult(cc.B(), un, up);
386  nnup = gmm::vect_norm2(up);
387  }
388  new_approx->add_point(c.xreal(), pai->coeff_on_face(f, j)
389  * gmm::abs(c.J()) * nup * nnup, ff);
390  }
391  }
392  }
393  }
394 
395  if (new_approx->nb_points()) {
396  new_approx->valid_method();
397  pintegration_method
398  pim = std::make_shared<integration_method>(new_approx);
399  dal::pstatic_stored_object_key
400  pk = std::make_shared<special_imls_key>(new_approx);
401  dal::add_stored_object(pk, pim, new_approx->ref_convex(),
402  new_approx->pintegration_points());
403  build_methods.push_back(pim);
404  cut_im.set_integration_method(cv, pim);
405  }
406  }
407 
408  void mesh_im_level_set::adapt(void) {
409  GMM_ASSERT1(linked_mesh_ != 0, "mesh level set uninitialized");
410  context_check();
411  clear_build_methods();
412  ignored_im.clear();
413  for (dal::bv_visitor cv(linked_mesh().convex_index());
414  !cv.finished(); ++cv) {
415  if (mls->is_convex_cut(cv)) build_method_of_convex(cv);
416 
417  if (!cut_im.convex_index().is_in(cv)) {
418  /* not exclusive with mls->is_convex_cut ... sometimes, cut cv
419  contains no integration points.. */
420 
421  if (integrate_where == INTEGRATE_BOUNDARY) {
422  ignored_im.add(cv);
423  } else if (integrate_where != (INTEGRATE_OUTSIDE|INTEGRATE_INSIDE)) {
424  /* remove convexes that are not in the integration area */
425  std::vector<pmesher_signed_distance> mesherls0(mls->nb_level_sets());
426  std::vector<pmesher_signed_distance> mesherls1(mls->nb_level_sets());
427  for (unsigned i = 0; i < mls->nb_level_sets(); ++i) {
428  mesherls0[i] = mls->get_level_set(i)->mls_of_convex(cv, 0, false);
429  if (mls->get_level_set(i)->has_secondary())
430  mesherls1[i] = mls->get_level_set(i)->mls_of_convex(cv,1, false);
431  }
432 
433  base_node B(gmm::mean_value(linked_mesh().trans_of_convex(cv)
434  ->convex_ref()->points()));
435  if (!is_point_in_selected_area(mesherls0, mesherls1, B).in)
436  ignored_im.add(cv);
437  }
438  }
439  }
440  is_adapted = true; touch();
441  // cout << "Number of built methods : " << build_methods.size() << endl;
442  }
443 
444  void mesh_im_level_set::compute_normal_vector
445  (const fem_interpolation_context &ctx, base_small_vector &vec) const {
446  size_type nb_ls = mls->nb_level_sets(), j = 0;
447  std::vector<pmesher_signed_distance> mesherls0(nb_ls);
448  if (vec.size() != linked_mesh().dim()) vec.resize(linked_mesh().dim());
449  base_small_vector un(ctx.pgt()->dim());
450 
451  if (nb_ls == 0) {
452  gmm::clear(vec); return;
453  } else if (nb_ls == 1) {
454  mesherls0[0]
455  = mls->get_level_set(0)->mls_of_convex(ctx.convex_num(), 0, false);
456  } else {
457  scalar_type d_min(0);
458  for (unsigned i = 0; i < nb_ls; ++i) {
459  mesherls0[i]
460  = mls->get_level_set(i)->mls_of_convex(ctx.convex_num(), 0, false);
461  scalar_type d = gmm::abs((*(mesherls0[i]))(ctx.xref()));
462  if (i == 0 || d < d_min) { d_min = d; j = i; }
463  }
464  }
465  mesherls0[j]->grad(ctx.xref(), un);
466  gmm::mult(ctx.B(), un, vec);
467  scalar_type no = gmm::vect_norm2(vec);
468  if (no != scalar_type(0))
469  gmm::scale(vec, scalar_type(1) / gmm::vect_norm2(vec));
470  }
471 
472 
473 
474  void mesh_im_cross_level_set::update_from_context(void) const
475  { is_adapted = false; }
476 
477  void mesh_im_cross_level_set::clear_build_methods() {
478  for (size_type i = 0; i < build_methods.size(); ++i)
479  if (build_methods[i].get()) del_stored_object(build_methods[i]);
480  build_methods.clear();
481  cut_im.clear();
482  }
483 
484  void mesh_im_cross_level_set::clear(void)
485  { mesh_im::clear(); clear_build_methods(); is_adapted = false; }
486 
487  void mesh_im_cross_level_set::init_with_mls(mesh_level_set &me,
488  size_type ind_ls1_, size_type ind_ls2_,
489  pintegration_method pim) {
490  init_with_mesh(me.linked_mesh());
491  cut_im.init_with_mesh(me.linked_mesh());
492  mls = &me;
493  ind_ls1 = ind_ls1_; ind_ls2 = ind_ls2_;
494  set_segment_im(pim);
495  this->add_dependency(*mls);
496  is_adapted = false;
497  }
498 
499  mesh_im_cross_level_set::mesh_im_cross_level_set(mesh_level_set &me,
500  size_type ind_ls1_, size_type ind_ls2_,
501  pintegration_method pim)
502  { mls = 0; init_with_mls(me, ind_ls1_, ind_ls2_, pim); }
503 
504  mesh_im_cross_level_set::mesh_im_cross_level_set(void)
505  { mls = 0; is_adapted = false; }
506 
507 
508  pintegration_method
509  mesh_im_cross_level_set::int_method_of_element(size_type cv) const {
510  if (!is_adapted) const_cast<mesh_im_cross_level_set *>(this)->adapt();
511  if (cut_im.convex_index().is_in(cv))
512  return cut_im.int_method_of_element(cv);
513  else {
514  if (ignored_im.is_in(cv)) return getfem::im_none();
515 
516  return mesh_im::int_method_of_element(cv);
517  }
518  }
519 
520  static bool is_point_in_intersection
521  (const std::vector<pmesher_signed_distance> &mesherls0,
522  const std::vector<pmesher_signed_distance> &mesherls1,
523  const base_node& P) {
524 
525  bool r = true;
526  for (unsigned i = 0; i < mesherls0.size(); ++i) {
527  bool sec = (dynamic_cast<const mesher_level_set *>(mesherls1[i].get()))->is_initialized();
528  scalar_type d1 = (*(mesherls0[i]))(P);
529  scalar_type d2 = (sec ? (*(mesherls1[i]))(P) : -1);
530  if (!(gmm::abs(d1) < 1e-7 && d2 < 1e-7)) r = false;
531  }
532  return r;
533  }
534 
535  static bool is_edges_intersect(const base_node &PP1, const base_node &PP2,
536  const base_node &PR1, const base_node &PR2) {
537  size_type n = gmm::vect_size(PP1), k1 = 0;
538  scalar_type c1 = scalar_type(0);
539  base_node V = PR2 - PR1;
540  for (size_type k = 0; k < n; ++k)
541  if (gmm::abs(V[k]) > gmm::abs(c1)) { c1 = V[k]; k1 = k; }
542 
543  scalar_type alpha1 = (PP1[k1] - PR1[k1]) / c1;
544  scalar_type alpha2 = (PP2[k1] - PR1[k1]) / c1;
545  base_node W1 = PP1 - PR1 - alpha1 * V;
546  base_node W2 = PP2 - PR1 - alpha2 * V;
547  if (gmm::vect_norm2(W1) > 1e-7*gmm::vect_norm2(V)) return false;
548  if (gmm::vect_norm2(W2) > 1e-7*gmm::vect_norm2(V)) return false;
549  if (alpha1 > 1.-1e-7 && alpha2 > 1.-1e-7) return false;
550  if (alpha1 < 1e-7 && alpha2 < 1e-7) return false;
551  return true;
552  }
553 
554 
555  void mesh_im_cross_level_set::build_method_of_convex
556  (size_type cv, mesh &global_intersection, bgeot::rtree &rtree_seg) {
557  const mesh &msh(mls->mesh_of_convex(cv));
558  GMM_ASSERT3(msh.convex_index().card() != 0, "Internal error");
559  base_matrix G;
560  base_node B;
561 
562  std::vector<pmesher_signed_distance> mesherls0(2);
563  std::vector<pmesher_signed_distance> mesherls1(2);
564  dal::bit_vector convexes_arein;
565 
566  mesherls0[0] = mls->get_level_set(ind_ls1)->mls_of_convex(cv, 0, false);
567  mesherls0[1] = mls->get_level_set(ind_ls2)->mls_of_convex(cv, 0, false);
568  if (mls->get_level_set(ind_ls1)->has_secondary())
569  mesherls1[0] = mls->get_level_set(ind_ls1)->mls_of_convex(cv, 1, false);
570  if (mls->get_level_set(ind_ls2)->has_secondary())
571  mesherls1[1] = mls->get_level_set(ind_ls2)->mls_of_convex(cv, 1, false);
572 
573  bgeot::pgeometric_trans pgt = linked_mesh().trans_of_convex(cv);
574  bgeot::pgeometric_trans pgt2
575  = msh.trans_of_convex(msh.convex_index().first_true());
576  dim_type n = pgt->dim();
577 
578  auto new_approx = std::make_shared<approx_integration>(pgt->convex_ref());
579  new_approx->set_built_on_the_fly();
580  base_matrix KK(n,n), CS(n,n);
581  base_matrix pc(pgt2->nb_points(), n);
582  std::vector<size_type> ptsing;
583 
584  for (dal::bv_visitor i(msh.convex_index()); !i.finished(); ++i) {
585  papprox_integration pai = segment_pim->approx_method();
586  GMM_ASSERT1(gmm::vect_size(pai->point(0)) == 1,
587  "A segment integration method is needed");
588 
589  base_matrix G2;
590  vectors_to_base_matrix(G2, linked_mesh().points_of_convex(cv));
591  bgeot::geotrans_interpolation_context
592  cc(linked_mesh().trans_of_convex(cv), base_node(n), G2);
593 
594  dal::bit_vector ptinter;
595  for (short_type k = 0; k < n; ++k) {
596  size_type ipt = msh.structure_of_convex(i)->ind_dir_points()[k];
597  const base_node &P = msh.points_of_convex(i)[ipt];
598  if (is_point_in_intersection(mesherls0, mesherls1, P))
599  ptinter.add(ipt);
600  }
601 
602  switch (n) {
603  case 2:
604  for (short_type k = 0; k < n; ++k) {
605  size_type ipt = msh.structure_of_convex(i)->ind_dir_points()[k];
606  if (ptinter.is_in(ipt)) {
607 
608  const base_node &P = msh.points_of_convex(i)[ipt];
609  cc.set_xref(P);
610 
611  if (global_intersection.search_point(cc.xreal())
612  == size_type(-1)) {
613  global_intersection.add_point(cc.xreal());
614  new_approx->add_point(msh.points_of_convex(i)[ipt],
615  scalar_type(1));
616  }
617 
618  }
619  }
620  break;
621  case 3:
622  {
623  for (short_type k1 = 1; k1 < n; ++k1) {
624  size_type ipt1 = msh.structure_of_convex(i)->ind_dir_points()[k1];
625  for (short_type k2 = 0; k2 < k1; ++k2) {
626  size_type ipt2=msh.structure_of_convex(i)->ind_dir_points()[k2];
627  if (ptinter.is_in(ipt1) && ptinter.is_in(ipt2)) {
628 
629  const base_node &P1 = msh.points_of_convex(i)[ipt1];
630  const base_node &P2 = msh.points_of_convex(i)[ipt2];
631  cc.set_xref(P1);
632  base_node PR1 = cc.xreal();
633  cc.set_xref(P2);
634  base_node PR2 = cc.xreal();
635 
636  size_type i1 = global_intersection.search_point(PR1);
637  size_type i2 = global_intersection.search_point(PR2);
638 
639  if (i1 == size_type(-1) || i2 == size_type(-1) ||
640  !global_intersection.nb_convex_with_edge(i1, i2)) {
641 
642  base_node min(n), max(n);
643  for (size_type k = 0; k < n; ++k) {
644  min[k] = std::min(PR1[k], PR2[k]);
645  max[k] = std::max(PR1[k], PR2[k]);
646  }
647  bgeot::rtree::pbox_set boxlst;
648  rtree_seg.find_intersecting_boxes(min, max, boxlst);
649 
650  bool found_intersect = false;
651 
652  for (bgeot::rtree::pbox_set::const_iterator
653  it=boxlst.begin(); it != boxlst.end(); ++it) {
654  const base_node &PP1
655  = global_intersection.points_of_convex((*it)->id)[0];
656  const base_node &PP2
657  = global_intersection.points_of_convex((*it)->id)[1];
658  if (is_edges_intersect(PP1, PP2, PR1, PR2))
659  { found_intersect = true; break; }
660  }
661 
662  if (!found_intersect) {
663 
664  i1 = global_intersection.add_point(PR1);
665  i2 = global_intersection.add_point(PR2);
666 
667  size_type is = global_intersection.add_segment(i1, i2);
668 
669  rtree_seg.add_box(min, max, is);
670  rtree_seg.build_tree(); // Not efficient !
671 
672 
673  const base_node &PE1
674  = msh.trans_of_convex(i)->convex_ref()->points()[ipt1];
675  const base_node &PE2
676  = msh.trans_of_convex(i)->convex_ref()->points()[ipt1];
677  base_node V = PE2 - PE1, W1(n), W2(n);
678 
679  base_matrix G3;
680  vectors_to_base_matrix(G3, msh.points_of_convex(i));
681  bgeot::geotrans_interpolation_context
682  ccc(msh.trans_of_convex(i), base_node(n), G3);
683 
684  for (size_type j=0; j < pai->nb_points_on_convex(); ++j) {
685  base_node PE = pai->point(j)[0] * PE2
686  + (scalar_type(1) - pai->point(j)[0]) * PE1;
687  ccc.set_xref(PE);
688  cc.set_xref(ccc.xreal());
689  gmm::mult(ccc.K(), V, W1);
690  gmm::mult(cc.K(), W1, W2);
691  new_approx->add_point(ccc.xreal(),
692  pai->coeff(j) * gmm::vect_norm2(W2));
693  }
694  }
695  }
696  }
697  }
698  }
699  }
700  break;
701  default: GMM_ASSERT1(false, "internal error");
702 
703  }
704  }
705 
706  if (new_approx->nb_points()) {
707  new_approx->valid_method();
708  pintegration_method
709  pim = std::make_shared<integration_method>(new_approx);
710  dal::pstatic_stored_object_key
711  pk = std::make_shared<special_imls_key>(new_approx);
712  dal::add_stored_object(pk, pim, new_approx->ref_convex(),
713  new_approx->pintegration_points());
714  build_methods.push_back(pim);
715  cut_im.set_integration_method(cv, pim);
716  }
717  }
718 
719  void mesh_im_cross_level_set::adapt(void) {
720  GMM_ASSERT1(linked_mesh_ != 0, "mesh level set uninitialized");
721  GMM_ASSERT1(linked_mesh_->dim() > 1 && linked_mesh_->dim() <= 3,
722  "Sorry, works only in dimension 2 or 3");
723 
724  context_check();
725  clear_build_methods();
726  ignored_im.clear();
727  mesh global_intersection;
728  bgeot::rtree rtree_seg;
729 
730  std::vector<size_type> icv;
731  std::vector<dal::bit_vector> ils;
732  mls->find_level_set_potential_intersections(icv, ils);
733 
734  for (size_type i = 0; i < icv.size(); ++i) {
735  if (ils[i].is_in(ind_ls1) && ils[i].is_in(ind_ls2)) {
736  build_method_of_convex(icv[i], global_intersection, rtree_seg);
737  }
738  }
739 
740  for (dal::bv_visitor cv(linked_mesh().convex_index());
741  !cv.finished(); ++cv)
742  if (!cut_im.convex_index().is_in(cv)) ignored_im.add(cv);
743 
744  is_adapted = true; touch();
745  }
746 
747 
748 } /* end of namespace getfem. */
749 
750 
751 
bgeot_kdtree.h
Simple implementation of a KD-tree.
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_im::convex_index
const dal::bit_vector & convex_index(void) const
Get the set of convexes where an integration method has been assigned.
Definition: getfem_mesh_im.h:73
getfem::mesh_im_cross_level_set::int_method_of_element
virtual pintegration_method int_method_of_element(size_type cv) const
return the integration method associated with an element (in no integration is associated,...
Definition: getfem_mesh_im_level_set.cc:509
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::int_method_descriptor
pintegration_method int_method_descriptor(std::string name, bool throw_if_not_found=true)
Get an integration method from its name .
Definition: getfem_integration.cc:1130
getfem::mesh_im_level_set
Describe an adaptable integration method linked to a mesh cut by a level set.
Definition: getfem_mesh_im_level_set.h:56
getfem::mesh_im_cross_level_set::set_segment_im
void set_segment_im(pintegration_method pim)
Set the specific integration methods.
Definition: getfem_mesh_im_level_set.h:219
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_mesher.h
An experimental mesher.
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
bgeot::geotrans_interpolation_context::xref
const base_node & xref() const
coordinates of the current point, in the reference convex.
Definition: bgeot_geometric_trans.cc:291
getfem::fem_interpolation_context::convex_num
size_type convex_num() const
get the current convex number
Definition: getfem_fem.cc:52
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
getfem::name_of_int_method
std::string name_of_int_method(pintegration_method p)
Get the string name of an integration method .
Definition: getfem_integration.cc:1137
bgeot::rtree
Balanced tree of n-dimensional rectangles.
Definition: bgeot_rtree.h:98
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::fem_interpolation_context
structure passed as the argument of fem interpolation functions.
Definition: getfem_fem.h:749
getfem::mesh_im_level_set::int_method_of_element
virtual pintegration_method int_method_of_element(size_type cv) const
return the integration method associated with an element (in no integration is associated,...
Definition: getfem_mesh_im_level_set.cc:70
getfem::mesh_im_cross_level_set
Describe an adaptable integration method linked to a mesh cut by at least two level sets on the inter...
Definition: getfem_mesh_im_level_set.h:186
getfem::mesh_level_set
Keep informations about a mesh crossed by level-sets.
Definition: getfem_mesh_level_set.h:52
getfem::mesh_im_cross_level_set::update_from_context
void update_from_context(void) const
this function has to be defined and should update the object when the context is modified.
Definition: getfem_mesh_im_level_set.cc:474
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:557
getfem::mesh_im_level_set::update_from_context
void update_from_context(void) const
this function has to be defined and should update the object when the context is modified.
Definition: getfem_mesh_im_level_set.cc:28
dal::add_stored_object
void add_stored_object(pstatic_stored_object_key k, pstatic_stored_object o, permanence perm)
Add an object with two optional dependencies.
Definition: dal_static_stored_objects.cc:284
getfem::im_none
pintegration_method im_none(void)
return IM_NONE
Definition: getfem_integration.cc:1293
getfem::mesh_im_cross_level_set::adapt
void adapt(void)
Apply the adequate integration methods.
Definition: getfem_mesh_im_level_set.cc:719
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
getfem::mesh_im_level_set::set_simplex_im
void set_simplex_im(pintegration_method reg, pintegration_method sing=0)
Set the specific integration methods.
Definition: getfem_mesh_im_level_set.h:109
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_im::linked_mesh
const mesh & linked_mesh() const
Give a reference to the linked mesh of type mesh.
Definition: getfem_mesh_im.h:79
getfem::mesh_im::int_method_of_element
virtual pintegration_method int_method_of_element(size_type cv) const
return the integration method associated with an element (in no integration is associated,...
Definition: getfem_mesh_im.h:117
getfem::context_dependencies::context_check
bool context_check() const
return true if update_from_context was called
Definition: getfem_context.h:126
getfem_mesh_im_level_set.h
a subclass of mesh_im which is conformal to a number of level sets.
getfem::mesh_im_level_set::adapt
void adapt(void)
Apply the adequate integration methods.
Definition: getfem_mesh_im_level_set.cc:408
getfem::mesh_im::set_integration_method
void set_integration_method(size_type cv, pintegration_method pim)
Set the integration method of a convex.
Definition: getfem_mesh_im.cc:49