[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Mon, 16 Oct 2023 09:10:06 -0400 (EDT) |
branch: master
commit 88a69becf3d33fa2b55d083a76e88a5033bbaab4
Author: Konstantinos Poulios <logari81@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. */