branch: master commit 88a69becf3d33fa2b55d083a76e88a5033bbaab4 Author: Konstantinos Poulios <logar...@gmail.com> AuthorDate: Mon Oct 16 14:37:11 2023 +0200
Whitespace changes --- src/getfem_mesh_im_level_set.cc | 694 ++++++++++++++++++++-------------------- 1 file changed, 347 insertions(+), 347 deletions(-) diff --git a/src/getfem_mesh_im_level_set.cc b/src/getfem_mesh_im_level_set.cc index 00edc831..11eae1b6 100644 --- a/src/getfem_mesh_im_level_set.cc +++ b/src/getfem_mesh_im_level_set.cc @@ -39,12 +39,12 @@ namespace getfem { mesh_im::clear(); clear_build_methods(); is_adapted = false; - } + } - void mesh_im_level_set::init_with_mls(mesh_level_set &me, - int integrate_where_, - pintegration_method reg, - pintegration_method sing) { + void mesh_im_level_set::init_with_mls(mesh_level_set &me, + int integrate_where_, + pintegration_method reg, + pintegration_method sing) { init_with_mesh(me.linked_mesh()); cut_im.init_with_mesh(me.linked_mesh()); mls = &me; @@ -55,9 +55,9 @@ namespace getfem { } mesh_im_level_set::mesh_im_level_set(mesh_level_set &me, - int integrate_where_, - pintegration_method reg, - pintegration_method sing) { + int integrate_where_, + pintegration_method reg, + pintegration_method sing) { mls = 0; init_with_mls(me, integrate_where_, reg, sing); } @@ -66,14 +66,14 @@ namespace getfem { { mls = 0; is_adapted = false; } - pintegration_method + pintegration_method mesh_im_level_set::int_method_of_element(size_type cv) const { if (!is_adapted) const_cast<mesh_im_level_set *>(this)->adapt(); - if (cut_im.convex_index().is_in(cv)) - return cut_im.int_method_of_element(cv); + if (cut_im.convex_index().is_in(cv)) + return cut_im.int_method_of_element(cv); else { if (ignored_im.is_in(cv)) //integrate_where == INTEGRATE_BOUNDARY) - return getfem::im_none(); + return getfem::im_none(); return mesh_im::int_method_of_element(cv); } @@ -91,16 +91,16 @@ namespace getfem { for (unsigned i = 0; i < mls->nb_level_sets(); ++i) { isin = isin && ((*(mesherls0[i]))(P) < 0); if (gmm::abs((*(mesherls0[i]))(P)) < 1e-7) - isbin = i+1; + isbin = i+1; if (mls->get_level_set(i)->has_secondary()) - isin = isin && ((*(mesherls1[i]))(P) < 0); + isin = isin && ((*(mesherls1[i]))(P) < 0); } - bool2 b; + bool2 b; b.in = ((integrate_where & INTEGRATE_OUTSIDE)) ? !isin : isin; b.bin = isbin; return b; } - + /* very rustic set operations evaluator */ struct is_in_eval { @@ -110,77 +110,77 @@ namespace getfem { bool2 do_expr(const char *&s) { bool2 r; if (*s == '(') { - r = do_expr(++s); - GMM_ASSERT1(*s++ == ')', - "expecting ')' in csg expression at '" << s-1 << "'"); + r = do_expr(++s); + GMM_ASSERT1(*s++ == ')', + "expecting ')' in csg expression at '" << s-1 << "'"); } else if (*s == '!') { // complementary - r = do_expr(++s); r.in = !r.in; + r = do_expr(++s); r.in = !r.in; } else if (*s >= 'a' && *s <= 'z') { - unsigned idx = (*s) - 'a'; - r.in = in.is_in(idx); - r.bin = bin.is_in(idx) ? idx+1 : 0; - ++s; - } else - GMM_ASSERT1(false, "parse error in csg expression at '" << s << "'"); + unsigned idx = (*s) - 'a'; + r.in = in.is_in(idx); + r.bin = bin.is_in(idx) ? idx+1 : 0; + ++s; + } else + GMM_ASSERT1(false, "parse error in csg expression at '" << s << "'"); if (*s == '+') { // Union - //cerr << "s = " << s << ", r = " << r << "\n"; - bool2 a = r, b = do_expr(++s); - //cerr << "->b = " << b << "\n"; - r.in = b.in || a.in; - if (b.bin && !a.in) r.bin = b.bin; - else if (a.bin && !b.in) r.bin = a.bin; - else r.bin = 0; + //cerr << "s = " << s << ", r = " << r << "\n"; + bool2 a = r, b = do_expr(++s); + //cerr << "->b = " << b << "\n"; + r.in = b.in || a.in; + if (b.bin && !a.in) r.bin = b.bin; + else if (a.bin && !b.in) r.bin = a.bin; + else r.bin = 0; } else if (*s == '-') { // Set difference - bool2 a = r, b = do_expr(++s); - r.in = a.in && !b.in; - if (a.bin && !b.in) r.bin = a.bin; - else if (a.in && b.bin) r.bin = b.bin; - else r.bin = 0; + bool2 a = r, b = do_expr(++s); + r.in = a.in && !b.in; + if (a.bin && !b.in) r.bin = a.bin; + else if (a.in && b.bin) r.bin = b.bin; + else r.bin = 0; } else if (*s == '*') { // Intersection - bool2 a = r, b = do_expr(++s); - r.in = a.in && b.in; - if (a.bin && b.in) r.bin = a.bin; - else if (a.in && b.bin) r.bin = b.bin; - else r.bin = 0; + bool2 a = r, b = do_expr(++s); + r.in = a.in && b.in; + if (a.bin && b.in) r.bin = a.bin; + else if (a.in && b.bin) r.bin = b.bin; + else r.bin = 0; } return r; } - bool2 is_in(const char*s) { - bool2 b = do_expr(s); + bool2 is_in(const char*s) { + bool2 b = do_expr(s); GMM_ASSERT1(!(*s), "parse error in CSG expression at " << s); return b; } void check() { const char *s[] = { "a*b*c*d", - "a+b+c+d", - "(a+b+c+d)", - "d*(a+b+c)", - "(a+b)-(c+d)", - "((a+b)-(c+d))", - "!a", - 0 }; - for (const char **p = s; *p; ++p) - cerr << *p << "\n"; + "a+b+c+d", + "(a+b+c+d)", + "d*(a+b+c)", + "(a+b)-(c+d)", + "((a+b)-(c+d))", + "!a", + 0 }; + for (const char **p = s; *p; ++p) + cerr << *p << "\n"; for (unsigned c=0; c < 16; ++c) { - in[0] = (c&1); bin[0] = 1; - in[1] = (c&2); bin[1] = 1; - in[2] = (c&4); bin[2] = 1; - in[3] = (c&8); bin[3] = 1; - cerr << in[0] << in[1] << in[2] << in[3] << ": "; - for (const char **p = s; *p; ++p) { - bool2 b = is_in(*p); - cerr << b.in << "/" << b.bin << " "; - } - cerr << "\n"; + in[0] = (c&1); bin[0] = 1; + in[1] = (c&2); bin[1] = 1; + in[2] = (c&4); bin[2] = 1; + in[3] = (c&8); bin[3] = 1; + cerr << in[0] << in[1] << in[2] << in[3] << ": "; + for (const char **p = s; *p; ++p) { + bool2 b = is_in(*p); + cerr << b.in << "/" << b.bin << " "; + } + cerr << "\n"; } } }; - - mesh_im_level_set::bool2 + + mesh_im_level_set::bool2 mesh_im_level_set::is_point_in_selected_area (const std::vector<pmesher_signed_distance> &mesherls0, - const std::vector<pmesher_signed_distance> &mesherls1, - const base_node& P) { + const std::vector<pmesher_signed_distance> &mesherls1, + const base_node& P) { is_in_eval ev; for (unsigned i = 0; i < mls->nb_level_sets(); ++i) { bool sec = mls->get_level_set(i)->has_secondary(); @@ -188,12 +188,12 @@ namespace getfem { scalar_type d2 = (sec ? (*(mesherls1[i]))(P) : -1); if (d1 < 0 && d2 < 0) ev.in.add(i); // if ((integrate_where & INTEGRATE_OUTSIDE) /*&& !sec*/) - // ev.in[i].flip(); + // ev.in[i].flip(); - if (gmm::abs(d1) < 1e-7 && d2 < 1e-7) - ev.bin.add(i); + if (gmm::abs(d1) < 1e-7 && d2 < 1e-7) + ev.bin.add(i); } - + bool2 r; if (ls_csg_description.size()) @@ -204,20 +204,20 @@ namespace getfem { } if (integrate_where & INTEGRATE_OUTSIDE) r.in = !(r.in); - + /*bool2 r2 = is_point_in_selected_area2(mesherls0,mesherls1,P); if (r2.in != r.in || r2.bin != r.bin) { cerr << "ev.in = " << ev.in << ", bin=" << ev.bin<<"\n"; cerr << "is_point_in_selected_area2("<<P <<"): r="<<r.in<<"/"<<r.bin - << ", r2=" << r2.in<<"/"<<r2.bin <<"\n"; + << ", r2=" << r2.in<<"/"<<r2.bin <<"\n"; assert(0); }*/ return r; } - + void mesh_im_level_set::build_method_of_convex(size_type cv) { const mesh &msh(mls->mesh_of_convex(cv)); GMM_ASSERT3(msh.convex_index().card() != 0, "Internal error"); @@ -232,17 +232,17 @@ namespace getfem { for (unsigned i = 0; i < mls->nb_level_sets(); ++i) { mesherls0[i] = mls->get_level_set(i)->mls_of_convex(cv, 0, false); if (mls->get_level_set(i)->has_secondary()) - mesherls1[i] = mls->get_level_set(i)->mls_of_convex(cv, 1, false); + mesherls1[i] = mls->get_level_set(i)->mls_of_convex(cv, 1, false); } if (integrate_where != (INTEGRATE_ALL)) { for (dal::bv_visitor scv(msh.convex_index()); !scv.finished(); ++scv) { - B = gmm::mean_value(msh.points_of_convex(scv)); - convexes_arein[scv] = - is_point_in_selected_area(mesherls0, mesherls1, B).in; + B = gmm::mean_value(msh.points_of_convex(scv)); + convexes_arein[scv] = + is_point_in_selected_area(mesherls0, mesherls1, B).in; } } - + bgeot::pgeometric_trans pgt = linked_mesh().trans_of_convex(cv); bgeot::pgeometric_trans pgt2 = msh.trans_of_convex(msh.convex_index().first_true()); @@ -250,9 +250,9 @@ namespace getfem { if (base_singular_pim) GMM_ASSERT1 ((n != 2 || - base_singular_pim->structure()== bgeot::parallelepiped_structure(2)) + base_singular_pim->structure()== bgeot::parallelepiped_structure(2)) && (n != 3 - || base_singular_pim->structure() == bgeot::prism_P1_structure(3)) + || base_singular_pim->structure() == bgeot::prism_P1_structure(3)) && (n >= 2) && (n <= 3), "Base integration method for quasi polar integration not convenient"); @@ -266,140 +266,140 @@ namespace getfem { for (dal::bv_visitor i(msh.convex_index()); !i.finished(); ++i) { papprox_integration pai = regular_simplex_pim->approx_method(); - + 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"); - + if ((integrate_where != INTEGRATE_ALL) && - !convexes_arein[i]) continue; - + !convexes_arein[i]) continue; + if (base_singular_pim && mls->crack_tip_convexes().is_in(cv)) { - ptsing.resize(0); - unsigned sing_ls = unsigned(-1); - - for (unsigned ils = 0; ils < mls->nb_level_sets(); ++ils) - if (mls->get_level_set(ils)->has_secondary()) { - for (unsigned ipt = 0; ipt <= n; ++ipt) { - if (gmm::abs((*(mesherls0[ils]))(msh.points_of_convex(i)[ipt])) - < 1E-10 - && gmm::abs((*(mesherls1[ils]))(msh.points_of_convex(i)[ipt])) - < 1E-10) { - if (sing_ls == unsigned(-1)) sing_ls = ils; - GMM_ASSERT1(sing_ls == ils, "Two singular point in one " - "sub element : " << sing_ls << ", " << ils << - ". To be done."); - ptsing.push_back(ipt); - } - } - } - assert(ptsing.size() < n); - - if (ptsing.size() > 0) { - std::stringstream sts; - sts << "IM_QUASI_POLAR(" << name_of_int_method(base_singular_pim) - << ", " << ptsing[0]; - if (ptsing.size() > 1) sts << ", " << ptsing[1]; - sts << ")"; - pai = int_method_descriptor(sts.str())->approx_method(); - } + ptsing.resize(0); + unsigned sing_ls = unsigned(-1); + + for (unsigned ils = 0; ils < mls->nb_level_sets(); ++ils) + if (mls->get_level_set(ils)->has_secondary()) { + for (unsigned ipt = 0; ipt <= n; ++ipt) { + if (gmm::abs((*(mesherls0[ils]))(msh.points_of_convex(i)[ipt])) + < 1E-10 + && gmm::abs((*(mesherls1[ils]))(msh.points_of_convex(i)[ipt])) + < 1E-10) { + if (sing_ls == unsigned(-1)) sing_ls = ils; + GMM_ASSERT1(sing_ls == ils, "Two singular point in one " + "sub element : " << sing_ls << ", " << ils << + ". To be done."); + ptsing.push_back(ipt); + } + } + } + assert(ptsing.size() < n); + + if (ptsing.size() > 0) { + std::stringstream sts; + sts << "IM_QUASI_POLAR(" << name_of_int_method(base_singular_pim) + << ", " << ptsing[0]; + if (ptsing.size() > 1) sts << ", " << ptsing[1]; + sts << ")"; + pai = int_method_descriptor(sts.str())->approx_method(); + } } base_matrix G2; vectors_to_base_matrix(G2, linked_mesh().points_of_convex(cv)); bgeot::geotrans_interpolation_context - cc(linked_mesh().trans_of_convex(cv), pai->point(0), G2); + cc(linked_mesh().trans_of_convex(cv), pai->point(0), G2); if (integrate_where & (INTEGRATE_INSIDE | INTEGRATE_OUTSIDE)) { - - vectors_to_base_matrix(G, msh.points_of_convex(i)); - bgeot::geotrans_interpolation_context c(msh.trans_of_convex(i), - pai->point(0), G); - - for (size_type j = 0; j < pai->nb_points_on_convex(); ++j) { - c.set_xref(pai->point(j)); - pgt2->poly_vector_grad(pai->point(j), pc); - gmm::mult(G,pc,KK); - scalar_type J = gmm::lu_det(KK); - new_approx->add_point(c.xreal(), pai->coeff(j) * gmm::abs(J)); - - /*if (integrate_where == INTEGRATE_INSIDE) { - cc.set_xref(c.xreal()); - totof << cc.xreal()[0] << "\t" << cc.xreal()[1] << "\t" - << pai->coeff(j) * gmm::abs(J) << "\n"; - }*/ - } + + vectors_to_base_matrix(G, msh.points_of_convex(i)); + bgeot::geotrans_interpolation_context c(msh.trans_of_convex(i), + pai->point(0), G); + + for (size_type j = 0; j < pai->nb_points_on_convex(); ++j) { + c.set_xref(pai->point(j)); + pgt2->poly_vector_grad(pai->point(j), pc); + gmm::mult(G,pc,KK); + scalar_type J = gmm::lu_det(KK); + new_approx->add_point(c.xreal(), pai->coeff(j) * gmm::abs(J)); + + /*if (integrate_where == INTEGRATE_INSIDE) { + cc.set_xref(c.xreal()); + totof << cc.xreal()[0] << "\t" << cc.xreal()[1] << "\t" + << pai->coeff(j) * gmm::abs(J) << "\n"; + }*/ + } } // pgt2 = msh.trans_of_convex(i); for (short_type f = 0; f < pgt2->structure()->nb_faces(); ++f) { - short_type ff = short_type(-1); - unsigned isin = unsigned(-1); - - if (integrate_where == INTEGRATE_BOUNDARY) { - bool lisin = true; - for (unsigned ipt = 0; ipt < - pgt2->structure()->nb_points_of_face(f); ++ipt) { - const base_node &P = msh.points_of_face_of_convex(i, f)[ipt]; - isin = is_point_in_selected_area(mesherls0, mesherls1, P).bin; - //cerr << P << ":" << isin << " "; - if (!isin) { lisin = false; break; } - } - if (!lisin) continue; - else isin--; - } else { - B = gmm::mean_value(msh.points_of_face_of_convex(i, f)); - if (pgt->convex_ref()->is_in(B) < -1E-7) continue; - for (short_type fi = 0; fi < pgt->structure()->nb_faces(); ++fi) { - if (gmm::abs(pgt->convex_ref()->is_in_face(fi, B)) < 2E-6) ff = fi; - } - - if (ff == short_type(-1)) { - cout << "Distance to the element : " - << pgt->convex_ref()->is_in(B) << endl; - for (short_type fi = 0; fi < pgt->structure()->nb_faces(); ++fi) { - cout << "Distance to face " << fi << " : " - << gmm::abs(pgt->convex_ref()->is_in_face(fi, B)) << endl; - } - GMM_ASSERT3(false, "the point is neither in the interior nor " - "the boundary of the element"); - } - } - - vectors_to_base_matrix(G, msh.points_of_convex(i)); - bgeot::geotrans_interpolation_context c(msh.trans_of_convex(i), - pai->point(0), G); - - - for (size_type j = 0; j < pai->nb_points_on_face(f); ++j) { - if (gmm::abs(c.J()) > 1E-11) { - c.set_xref(pai->point_on_face(f, j)); - base_small_vector un = pgt2->normals()[f], up(msh.dim()); - gmm::mult(c.B(), un, up); - scalar_type nup = gmm::vect_norm2(up); - - scalar_type nnup(1); - if (integrate_where == INTEGRATE_BOUNDARY) { - cc.set_xref(c.xreal()); - mesherls0[isin]->grad(c.xreal(), un); - un /= gmm::vect_norm2(un); - gmm::mult(cc.B(), un, up); - nnup = gmm::vect_norm2(up); - } - new_approx->add_point(c.xreal(), pai->coeff_on_face(f, j) - * gmm::abs(c.J()) * nup * nnup, ff); - } - } + short_type ff = short_type(-1); + unsigned isin = unsigned(-1); + + if (integrate_where == INTEGRATE_BOUNDARY) { + bool lisin = true; + for (unsigned ipt = 0; ipt < + pgt2->structure()->nb_points_of_face(f); ++ipt) { + const base_node &P = msh.points_of_face_of_convex(i, f)[ipt]; + isin = is_point_in_selected_area(mesherls0, mesherls1, P).bin; + //cerr << P << ":" << isin << " "; + if (!isin) { lisin = false; break; } + } + if (!lisin) continue; + else isin--; + } else { + B = gmm::mean_value(msh.points_of_face_of_convex(i, f)); + if (pgt->convex_ref()->is_in(B) < -1E-7) continue; + for (short_type fi = 0; fi < pgt->structure()->nb_faces(); ++fi) { + if (gmm::abs(pgt->convex_ref()->is_in_face(fi, B)) < 2E-6) ff = fi; + } + + if (ff == short_type(-1)) { + cout << "Distance to the element : " + << pgt->convex_ref()->is_in(B) << endl; + for (short_type fi = 0; fi < pgt->structure()->nb_faces(); ++fi) { + cout << "Distance to face " << fi << " : " + << gmm::abs(pgt->convex_ref()->is_in_face(fi, B)) << endl; + } + GMM_ASSERT3(false, "the point is neither in the interior nor " + "the boundary of the element"); + } + } + + vectors_to_base_matrix(G, msh.points_of_convex(i)); + bgeot::geotrans_interpolation_context c(msh.trans_of_convex(i), + pai->point(0), G); + + + for (size_type j = 0; j < pai->nb_points_on_face(f); ++j) { + if (gmm::abs(c.J()) > 1E-11) { + c.set_xref(pai->point_on_face(f, j)); + base_small_vector un = pgt2->normals()[f], up(msh.dim()); + gmm::mult(c.B(), un, up); + scalar_type nup = gmm::vect_norm2(up); + + scalar_type nnup(1); + if (integrate_where == INTEGRATE_BOUNDARY) { + cc.set_xref(c.xreal()); + mesherls0[isin]->grad(c.xreal(), un); + un /= gmm::vect_norm2(un); + gmm::mult(cc.B(), un, up); + nnup = gmm::vect_norm2(up); + } + new_approx->add_point(c.xreal(), pai->coeff_on_face(f, j) + * gmm::abs(c.J()) * nup * nnup, ff); + } + } } } if (new_approx->nb_points()) { new_approx->valid_method(); pintegration_method - pim = std::make_shared<integration_method>(new_approx); + pim = std::make_shared<integration_method>(new_approx); dal::pstatic_stored_object_key - pk = std::make_shared<special_imls_key>(new_approx); + pk = std::make_shared<special_imls_key>(new_approx); dal::add_stored_object(pk, pim, new_approx->ref_convex(), - new_approx->pintegration_points()); + new_approx->pintegration_points()); build_methods.push_back(pim); cut_im.set_integration_method(cv, pim); } @@ -410,31 +410,31 @@ namespace getfem { context_check(); clear_build_methods(); ignored_im.clear(); - for (dal::bv_visitor cv(linked_mesh().convex_index()); - !cv.finished(); ++cv) { + for (dal::bv_visitor cv(linked_mesh().convex_index()); + !cv.finished(); ++cv) { if (mls->is_convex_cut(cv)) build_method_of_convex(cv); if (!cut_im.convex_index().is_in(cv)) { - /* not exclusive with mls->is_convex_cut ... sometimes, cut cv - contains no integration points.. */ - - if (integrate_where == INTEGRATE_BOUNDARY) { - ignored_im.add(cv); - } else if (integrate_where != (INTEGRATE_OUTSIDE|INTEGRATE_INSIDE)) { - /* remove convexes that are not in the integration area */ - std::vector<pmesher_signed_distance> mesherls0(mls->nb_level_sets()); - std::vector<pmesher_signed_distance> mesherls1(mls->nb_level_sets()); - for (unsigned i = 0; i < mls->nb_level_sets(); ++i) { - mesherls0[i] = mls->get_level_set(i)->mls_of_convex(cv, 0, false); - if (mls->get_level_set(i)->has_secondary()) - mesherls1[i] = mls->get_level_set(i)->mls_of_convex(cv,1, false); - } - - base_node B(gmm::mean_value(linked_mesh().trans_of_convex(cv) - ->convex_ref()->points())); - if (!is_point_in_selected_area(mesherls0, mesherls1, B).in) - ignored_im.add(cv); - } + /* not exclusive with mls->is_convex_cut ... sometimes, cut cv + contains no integration points.. */ + + if (integrate_where == INTEGRATE_BOUNDARY) { + ignored_im.add(cv); + } else if (integrate_where != (INTEGRATE_OUTSIDE|INTEGRATE_INSIDE)) { + /* remove convexes that are not in the integration area */ + std::vector<pmesher_signed_distance> mesherls0(mls->nb_level_sets()); + std::vector<pmesher_signed_distance> mesherls1(mls->nb_level_sets()); + for (unsigned i = 0; i < mls->nb_level_sets(); ++i) { + mesherls0[i] = mls->get_level_set(i)->mls_of_convex(cv, 0, false); + if (mls->get_level_set(i)->has_secondary()) + mesherls1[i] = mls->get_level_set(i)->mls_of_convex(cv,1, false); + } + + base_node B(gmm::mean_value(linked_mesh().trans_of_convex(cv) + ->convex_ref()->points())); + if (!is_point_in_selected_area(mesherls0, mesherls1, B).in) + ignored_im.add(cv); + } } } is_adapted = true; touch(); @@ -447,19 +447,19 @@ namespace getfem { std::vector<pmesher_signed_distance> mesherls0(nb_ls); if (vec.size() != linked_mesh().dim()) vec.resize(linked_mesh().dim()); base_small_vector un(ctx.pgt()->dim()); - + if (nb_ls == 0) { - gmm::clear(vec); return; + gmm::clear(vec); return; } else if (nb_ls == 1) { mesherls0[0] - = mls->get_level_set(0)->mls_of_convex(ctx.convex_num(), 0, false); + = mls->get_level_set(0)->mls_of_convex(ctx.convex_num(), 0, false); } else { scalar_type d_min(0); for (unsigned i = 0; i < nb_ls; ++i) { - mesherls0[i] - = mls->get_level_set(i)->mls_of_convex(ctx.convex_num(), 0, false); - scalar_type d = gmm::abs((*(mesherls0[i]))(ctx.xref())); - if (i == 0 || d < d_min) { d_min = d; j = i; } + mesherls0[i] + = mls->get_level_set(i)->mls_of_convex(ctx.convex_num(), 0, false); + scalar_type d = gmm::abs((*(mesherls0[i]))(ctx.xref())); + if (i == 0 || d < d_min) { d_min = d; j = i; } } } mesherls0[j]->grad(ctx.xref(), un); @@ -482,11 +482,11 @@ namespace getfem { } void mesh_im_cross_level_set::clear(void) - { mesh_im::clear(); clear_build_methods(); is_adapted = false; } + { mesh_im::clear(); clear_build_methods(); is_adapted = false; } - void mesh_im_cross_level_set::init_with_mls(mesh_level_set &me, - size_type ind_ls1_, size_type ind_ls2_, - pintegration_method pim) { + void mesh_im_cross_level_set::init_with_mls(mesh_level_set &me, + size_type ind_ls1_, size_type ind_ls2_, + pintegration_method pim) { init_with_mesh(me.linked_mesh()); cut_im.init_with_mesh(me.linked_mesh()); mls = &me; @@ -497,31 +497,31 @@ namespace getfem { } mesh_im_cross_level_set::mesh_im_cross_level_set(mesh_level_set &me, - size_type ind_ls1_, size_type ind_ls2_, - pintegration_method pim) + size_type ind_ls1_, size_type ind_ls2_, + pintegration_method pim) { mls = 0; init_with_mls(me, ind_ls1_, ind_ls2_, pim); } mesh_im_cross_level_set::mesh_im_cross_level_set(void) { mls = 0; is_adapted = false; } - pintegration_method + pintegration_method mesh_im_cross_level_set::int_method_of_element(size_type cv) const { if (!is_adapted) const_cast<mesh_im_cross_level_set *>(this)->adapt(); - if (cut_im.convex_index().is_in(cv)) - return cut_im.int_method_of_element(cv); + if (cut_im.convex_index().is_in(cv)) + return cut_im.int_method_of_element(cv); else { if (ignored_im.is_in(cv)) return getfem::im_none(); return mesh_im::int_method_of_element(cv); } } - + static bool is_point_in_intersection (const std::vector<pmesher_signed_distance> &mesherls0, const std::vector<pmesher_signed_distance> &mesherls1, const base_node& P) { - + bool r = true; for (unsigned i = 0; i < mesherls0.size(); ++i) { bool sec = (dynamic_cast<const mesher_level_set *>(mesherls1[i].get()))->is_initialized(); @@ -533,13 +533,13 @@ namespace getfem { } static bool is_edges_intersect(const base_node &PP1, const base_node &PP2, - const base_node &PR1, const base_node &PR2) { + const base_node &PR1, const base_node &PR2) { size_type n = gmm::vect_size(PP1), k1 = 0; scalar_type c1 = scalar_type(0); base_node V = PR2 - PR1; for (size_type k = 0; k < n; ++k) if (gmm::abs(V[k]) > gmm::abs(c1)) { c1 = V[k]; k1 = k; } - + scalar_type alpha1 = (PP1[k1] - PR1[k1]) / c1; scalar_type alpha2 = (PP2[k1] - PR1[k1]) / c1; base_node W1 = PP1 - PR1 - alpha1 * V; @@ -551,7 +551,7 @@ namespace getfem { return true; } - + void mesh_im_cross_level_set::build_method_of_convex (size_type cv, mesh &global_intersection, bgeot::rtree &rtree_seg) { const mesh &msh(mls->mesh_of_convex(cv)); @@ -569,7 +569,7 @@ namespace getfem { mesherls1[0] = mls->get_level_set(ind_ls1)->mls_of_convex(cv, 1, false); if (mls->get_level_set(ind_ls2)->has_secondary()) mesherls1[1] = mls->get_level_set(ind_ls2)->mls_of_convex(cv, 1, false); - + bgeot::pgeometric_trans pgt = linked_mesh().trans_of_convex(cv); bgeot::pgeometric_trans pgt2 = msh.trans_of_convex(msh.convex_index().first_true()); @@ -584,133 +584,133 @@ namespace getfem { for (dal::bv_visitor i(msh.convex_index()); !i.finished(); ++i) { papprox_integration pai = segment_pim->approx_method(); GMM_ASSERT1(gmm::vect_size(pai->point(0)) == 1, - "A segment integration method is needed"); + "A segment integration method is needed"); base_matrix G2; vectors_to_base_matrix(G2, linked_mesh().points_of_convex(cv)); bgeot::geotrans_interpolation_context - cc(linked_mesh().trans_of_convex(cv), base_node(n), G2); - + cc(linked_mesh().trans_of_convex(cv), base_node(n), G2); + dal::bit_vector ptinter; for (short_type k = 0; k < n; ++k) { - size_type ipt = msh.structure_of_convex(i)->ind_dir_points()[k]; - const base_node &P = msh.points_of_convex(i)[ipt]; - if (is_point_in_intersection(mesherls0, mesherls1, P)) - ptinter.add(ipt); + size_type ipt = msh.structure_of_convex(i)->ind_dir_points()[k]; + const base_node &P = msh.points_of_convex(i)[ipt]; + if (is_point_in_intersection(mesherls0, mesherls1, P)) + ptinter.add(ipt); } switch (n) { case 2: - for (short_type k = 0; k < n; ++k) { - size_type ipt = msh.structure_of_convex(i)->ind_dir_points()[k]; - if (ptinter.is_in(ipt)) { - - const base_node &P = msh.points_of_convex(i)[ipt]; - cc.set_xref(P); - - if (global_intersection.search_point(cc.xreal()) - == size_type(-1)) { - global_intersection.add_point(cc.xreal()); - new_approx->add_point(msh.points_of_convex(i)[ipt], - scalar_type(1)); - } - - } - } + for (short_type k = 0; k < n; ++k) { + size_type ipt = msh.structure_of_convex(i)->ind_dir_points()[k]; + if (ptinter.is_in(ipt)) { + + const base_node &P = msh.points_of_convex(i)[ipt]; + cc.set_xref(P); + + if (global_intersection.search_point(cc.xreal()) + == size_type(-1)) { + global_intersection.add_point(cc.xreal()); + new_approx->add_point(msh.points_of_convex(i)[ipt], + scalar_type(1)); + } + + } + } break; case 3: - { - for (short_type k1 = 1; k1 < n; ++k1) { - size_type ipt1 = msh.structure_of_convex(i)->ind_dir_points()[k1]; - for (short_type k2 = 0; k2 < k1; ++k2) { - size_type ipt2=msh.structure_of_convex(i)->ind_dir_points()[k2]; - if (ptinter.is_in(ipt1) && ptinter.is_in(ipt2)) { - - const base_node &P1 = msh.points_of_convex(i)[ipt1]; - const base_node &P2 = msh.points_of_convex(i)[ipt2]; - cc.set_xref(P1); - base_node PR1 = cc.xreal(); - cc.set_xref(P2); - base_node PR2 = cc.xreal(); - - size_type i1 = global_intersection.search_point(PR1); - size_type i2 = global_intersection.search_point(PR2); - - if (i1 == size_type(-1) || i2 == size_type(-1) || - !global_intersection.nb_convex_with_edge(i1, i2)) { - - base_node min(n), max(n); - for (size_type k = 0; k < n; ++k) { - min[k] = std::min(PR1[k], PR2[k]); - max[k] = std::max(PR1[k], PR2[k]); - } - bgeot::rtree::pbox_set boxlst; - rtree_seg.find_intersecting_boxes(min, max, boxlst); - - bool found_intersect = false; - - for (bgeot::rtree::pbox_set::const_iterator - it=boxlst.begin(); it != boxlst.end(); ++it) { - const base_node &PP1 - = global_intersection.points_of_convex((*it)->id)[0]; - const base_node &PP2 - = global_intersection.points_of_convex((*it)->id)[1]; - if (is_edges_intersect(PP1, PP2, PR1, PR2)) - { found_intersect = true; break; } - } - - if (!found_intersect) { - - i1 = global_intersection.add_point(PR1); - i2 = global_intersection.add_point(PR2); - - size_type is = global_intersection.add_segment(i1, i2); - - rtree_seg.add_box(min, max, is); + { + for (short_type k1 = 1; k1 < n; ++k1) { + size_type ipt1 = msh.structure_of_convex(i)->ind_dir_points()[k1]; + for (short_type k2 = 0; k2 < k1; ++k2) { + size_type ipt2=msh.structure_of_convex(i)->ind_dir_points()[k2]; + if (ptinter.is_in(ipt1) && ptinter.is_in(ipt2)) { + + const base_node &P1 = msh.points_of_convex(i)[ipt1]; + const base_node &P2 = msh.points_of_convex(i)[ipt2]; + cc.set_xref(P1); + base_node PR1 = cc.xreal(); + cc.set_xref(P2); + base_node PR2 = cc.xreal(); + + size_type i1 = global_intersection.search_point(PR1); + size_type i2 = global_intersection.search_point(PR2); + + if (i1 == size_type(-1) || i2 == size_type(-1) || + !global_intersection.nb_convex_with_edge(i1, i2)) { + + base_node min(n), max(n); + for (size_type k = 0; k < n; ++k) { + min[k] = std::min(PR1[k], PR2[k]); + max[k] = std::max(PR1[k], PR2[k]); + } + bgeot::rtree::pbox_set boxlst; + rtree_seg.find_intersecting_boxes(min, max, boxlst); + + bool found_intersect = false; + + for (bgeot::rtree::pbox_set::const_iterator + it=boxlst.begin(); it != boxlst.end(); ++it) { + const base_node &PP1 + = global_intersection.points_of_convex((*it)->id)[0]; + const base_node &PP2 + = global_intersection.points_of_convex((*it)->id)[1]; + if (is_edges_intersect(PP1, PP2, PR1, PR2)) + { found_intersect = true; break; } + } + + if (!found_intersect) { + + i1 = global_intersection.add_point(PR1); + i2 = global_intersection.add_point(PR2); + + size_type is = global_intersection.add_segment(i1, i2); + + rtree_seg.add_box(min, max, is); rtree_seg.build_tree(); // Not efficient ! - - const base_node &PE1 - = msh.trans_of_convex(i)->convex_ref()->points()[ipt1]; - const base_node &PE2 - = msh.trans_of_convex(i)->convex_ref()->points()[ipt1]; - base_node V = PE2 - PE1, W1(n), W2(n); - - base_matrix G3; - vectors_to_base_matrix(G3, msh.points_of_convex(i)); - bgeot::geotrans_interpolation_context - ccc(msh.trans_of_convex(i), base_node(n), G3); - - for (size_type j=0; j < pai->nb_points_on_convex(); ++j) { - base_node PE = pai->point(j)[0] * PE2 - + (scalar_type(1) - pai->point(j)[0]) * PE1; - ccc.set_xref(PE); - cc.set_xref(ccc.xreal()); - gmm::mult(ccc.K(), V, W1); - gmm::mult(cc.K(), W1, W2); - new_approx->add_point(ccc.xreal(), - pai->coeff(j) * gmm::vect_norm2(W2)); - } - } - } - } - } - } - } - break; + + const base_node &PE1 + = msh.trans_of_convex(i)->convex_ref()->points()[ipt1]; + const base_node &PE2 + = msh.trans_of_convex(i)->convex_ref()->points()[ipt1]; + base_node V = PE2 - PE1, W1(n), W2(n); + + base_matrix G3; + vectors_to_base_matrix(G3, msh.points_of_convex(i)); + bgeot::geotrans_interpolation_context + ccc(msh.trans_of_convex(i), base_node(n), G3); + + for (size_type j=0; j < pai->nb_points_on_convex(); ++j) { + base_node PE = pai->point(j)[0] * PE2 + + (scalar_type(1) - pai->point(j)[0]) * PE1; + ccc.set_xref(PE); + cc.set_xref(ccc.xreal()); + gmm::mult(ccc.K(), V, W1); + gmm::mult(cc.K(), W1, W2); + new_approx->add_point(ccc.xreal(), + pai->coeff(j) * gmm::vect_norm2(W2)); + } + } + } + } + } + } + } + break; default: GMM_ASSERT1(false, "internal error"); - + } } if (new_approx->nb_points()) { new_approx->valid_method(); pintegration_method - pim = std::make_shared<integration_method>(new_approx); + pim = std::make_shared<integration_method>(new_approx); dal::pstatic_stored_object_key - pk = std::make_shared<special_imls_key>(new_approx); + pk = std::make_shared<special_imls_key>(new_approx); dal::add_stored_object(pk, pim, new_approx->ref_convex(), - new_approx->pintegration_points()); + new_approx->pintegration_points()); build_methods.push_back(pim); cut_im.set_integration_method(cv, pim); } @@ -719,8 +719,8 @@ namespace getfem { void mesh_im_cross_level_set::adapt(void) { GMM_ASSERT1(linked_mesh_ != 0, "mesh level set uninitialized"); GMM_ASSERT1(linked_mesh_->dim() > 1 && linked_mesh_->dim() <= 3, - "Sorry, works only in dimension 2 or 3"); - + "Sorry, works only in dimension 2 or 3"); + context_check(); clear_build_methods(); ignored_im.clear(); @@ -733,18 +733,18 @@ namespace getfem { for (size_type i = 0; i < icv.size(); ++i) { if (ils[i].is_in(ind_ls1) && ils[i].is_in(ind_ls2)) { - build_method_of_convex(icv[i], global_intersection, rtree_seg); + build_method_of_convex(icv[i], global_intersection, rtree_seg); } } - for (dal::bv_visitor cv(linked_mesh().convex_index()); - !cv.finished(); ++cv) + for (dal::bv_visitor cv(linked_mesh().convex_index()); + !cv.finished(); ++cv) if (!cut_im.convex_index().is_in(cv)) ignored_im.add(cv); is_adapted = true; touch(); } - + } /* end of namespace getfem. */