branch: master commit e3372b14e48bb32d6271f6c4658799a922968712 Author: Konstantinos Poulios <logar...@gmail.com> AuthorDate: Fri Mar 8 00:14:06 2024 +0100
Move function implementations from header to source file --- src/getfem/getfem_mesh_slicers.h | 81 +++---------------------------- src/getfem_mesh_slicers.cc | 102 ++++++++++++++++++++++++++++++++++++++- 2 files changed, 109 insertions(+), 74 deletions(-) diff --git a/src/getfem/getfem_mesh_slicers.h b/src/getfem/getfem_mesh_slicers.h index e69e7078..d8f613e3 100644 --- a/src/getfem/getfem_mesh_slicers.h +++ b/src/getfem/getfem_mesh_slicers.h @@ -324,14 +324,7 @@ namespace getfem { slicer_volume(int orient_) : orient(orient_) {} /** Utility function */ - static scalar_type trinom(scalar_type a, scalar_type b, scalar_type c) { - scalar_type delta = b*b - 4*a*c; - if (delta < 0.) return 1./EPS; - delta = sqrt(delta); - scalar_type s1 = (-b - delta) / (2*a); - scalar_type s2 = (-b + delta) / (2*a); - if (gmm::abs(s1-.5) < gmm::abs(s2-.5)) return s1; else return s2; - } + static scalar_type trinom(scalar_type a, scalar_type b, scalar_type c); void split_simplex(mesh_slicer &ms, slice_simplex s, /* s is NOT a reference, it is on * purpose(push_back in the function)*/ @@ -346,22 +339,9 @@ namespace getfem { */ class slicer_half_space : public slicer_volume { const base_node x0, n; /* normal directed from inside toward outside */ - void test_point(const base_node& P, bool& in, bool& bound) const { - scalar_type s = gmm::vect_sp(P-x0,n); - in = (s <= 0); bound = (s*s <= EPS); // gmm::vect_norm2_sqr(P-x0)); No! - // do not try to be smart with the boundary check .. easily broken with - // slicer_mesh_with_mesh - } + void test_point(const base_node& P, bool& in, bool& bound) const; scalar_type edge_intersect(size_type iA, size_type iB, - const mesh_slicer::cs_nodes_ct& nodes) const { - const base_node& A=nodes[iA].pt; - const base_node& B=nodes[iB].pt; - scalar_type s1 = 0., s2 = 0.; - for (unsigned i=0; i < A.size(); ++i) - { s1 += (A[i] - B[i])*n[i]; s2 += (A[i]-x0[i])*n[i]; } - if (gmm::abs(s1) < EPS) return 1./EPS; - else return s2/s1; - } + const mesh_slicer::cs_nodes_ct& nodes) const; public: slicer_half_space(base_node x0_, base_node n_, int orient_) : slicer_volume(orient_), x0(x0_), n(n_/gmm::vect_norm2(n_)) { @@ -375,22 +355,9 @@ namespace getfem { class slicer_sphere : public slicer_volume { base_node x0; scalar_type R; - void test_point(const base_node& P, bool& in, bool& bound) const { - scalar_type R2 = gmm::vect_dist2_sqr(P,x0); - bound = (R2 >= (1-EPS)*R*R && R2 <= (1+EPS)*R*R); - in = R2 <= R*R; - } + void test_point(const base_node& P, bool& in, bool& bound) const; scalar_type edge_intersect(size_type iA, size_type iB, - const mesh_slicer::cs_nodes_ct& nodes) const { - const base_node& A=nodes[iA].pt; - const base_node& B=nodes[iB].pt; - scalar_type a,b,c; // a*x^2 + b*x + c = 0 - a = gmm::vect_norm2_sqr(B-A); - if (a < EPS) return pt_bin.is_in(iA) ? 0. : 1./EPS; - b = 2*gmm::vect_sp(A-x0,B-A); - c = gmm::vect_norm2_sqr(A-x0)-R*R; - return slicer_volume::trinom(a,b,c); - } + const mesh_slicer::cs_nodes_ct& nodes) const; public: /* orient = -1 => select interior, orient = 0 => select boundary @@ -406,35 +373,9 @@ namespace getfem { class slicer_cylinder : public slicer_volume { base_node x0, d; scalar_type R; - void test_point(const base_node& P, bool& in, bool& bound) const { - base_node N = P; - if (2 == N.size()) { /* Add Z coorinate if mesh is 2D */ - N.push_back(0.0); - } - N = N-x0; - scalar_type axpos = gmm::vect_sp(d, N); - scalar_type dist2 = gmm::vect_norm2_sqr(N) - gmm::sqr(axpos); - bound = gmm::abs(dist2-R*R) < EPS; - in = dist2 < R*R; - } + void test_point(const base_node& P, bool& in, bool& bound) const; scalar_type edge_intersect(size_type iA, size_type iB, - const mesh_slicer::cs_nodes_ct& nodes) const { - base_node F=nodes[iA].pt; - base_node D=nodes[iB].pt-nodes[iA].pt; - if (2 == F.size()) { - F.push_back(0.0); - D.push_back(0.0); - } - F = F - x0; - scalar_type Fd = gmm::vect_sp(F,d); - scalar_type Dd = gmm::vect_sp(D,d); - scalar_type a = gmm::vect_norm2_sqr(D) - gmm::sqr(Dd); - if (a < EPS) return pt_bin.is_in(iA) ? 0. : 1./EPS; - assert(a> -EPS); - scalar_type b = 2*(gmm::vect_sp(F,D) - Fd*Dd); - scalar_type c = gmm::vect_norm2_sqr(F) - gmm::sqr(Fd) - gmm::sqr(R); - return slicer_volume::trinom(a,b,c); - } + const mesh_slicer::cs_nodes_ct& nodes) const; public: slicer_cylinder(base_node x0_, base_node x1_,scalar_type R_,int orient_) : slicer_volume(orient_), x0(x0_), d(x1_-x0_), R(R_) { @@ -454,13 +395,7 @@ namespace getfem { void prepare(size_type cv, const mesh_slicer::cs_nodes_ct& nodes, const dal::bit_vector& nodes_index); scalar_type edge_intersect(size_type iA, size_type iB, - const mesh_slicer::cs_nodes_ct&) const { - assert(iA < Uval.size() && iB < Uval.size()); - if (((Uval[iA] < val) && (Uval[iB] > val)) || - ((Uval[iA] > val) && (Uval[iB] < val))) - return (val-Uval[iA])/(Uval[iB]-Uval[iA]); - else return 1./EPS; - } + const mesh_slicer::cs_nodes_ct&) const; public: /* orient = -1: u(x) <= val, 0: u(x) == val, +1: u(x) >= val */ slicer_isovalues(const mesh_slice_cv_dof_data_base& mfU_, diff --git a/src/getfem_mesh_slicers.cc b/src/getfem_mesh_slicers.cc index e636dc4b..34ac04f7 100644 --- a/src/getfem_mesh_slicers.cc +++ b/src/getfem_mesh_slicers.cc @@ -170,7 +170,19 @@ namespace getfem { // return true; //} - /* + scalar_type slicer_volume::trinom(scalar_type a, scalar_type b, scalar_type c) { + scalar_type delta = b * b - 4 * a * c; + if (delta < 0.) return 1. / EPS; + delta = sqrt(delta); + scalar_type s1 = (-b - delta) / (2 * a); + scalar_type s2 = (-b + delta) / (2 * a); + if (gmm::abs(s1 - 0.5) < gmm::abs(s2 - 0.5)) + return s1; + else + return s2; + } + + /* intersects the simplex with the slice, and (recursively) decomposes it into sub-simplices, which are added to the list 'splxs'. If orient == 0, then it is the faces of sub-simplices @@ -407,6 +419,17 @@ namespace getfem { } } + scalar_type + slicer_isovalues::edge_intersect(size_type iA, size_type iB, + const mesh_slicer::cs_nodes_ct&) const { + assert(iA < Uval.size() && iB < Uval.size()); + if (((Uval[iA] < val) && (Uval[iB] > val)) || + ((Uval[iA] > val) && (Uval[iB] < val))) + return (val-Uval[iA])/(Uval[iB]-Uval[iA]); + else + return 1./EPS; + } + void slicer_union::exec(mesh_slicer &ms) { dal::bit_vector splx_in_base = ms.splx_in; @@ -951,4 +974,81 @@ namespace getfem { apply_slicers(); } } + + void + slicer_half_space::test_point(const base_node& P, bool& in, bool& bound) const { + scalar_type s = gmm::vect_sp(P - x0, n); + in = (s <= 0); bound = (s * s <= EPS); // gmm::vect_norm2_sqr(P-x0)); No! + // do not try to be smart with the boundary check .. easily broken with + // slicer_mesh_with_mesh + } + + scalar_type + slicer_half_space::edge_intersect(size_type iA, size_type iB, + const mesh_slicer::cs_nodes_ct& nodes) const { + const base_node& A = nodes[iA].pt; + const base_node& B = nodes[iB].pt; + scalar_type s1 = 0., s2 = 0.; + for (unsigned i = 0; i < A.size(); ++i) { + s1 += (A[i] - B[i]) * n[i]; s2 += (A[i] - x0[i]) * n[i]; + } + if (gmm::abs(s1) < EPS) + return 1. / EPS; + else + return s2 / s1; + } + + void + slicer_sphere::test_point(const base_node& P, bool& in, bool& bound) const { + scalar_type R2 = gmm::vect_dist2_sqr(P,x0); + bound = (R2 >= (1-EPS)*R*R && R2 <= (1+EPS)*R*R); + in = R2 <= R*R; + } + + scalar_type + slicer_sphere::edge_intersect(size_type iA, size_type iB, + const mesh_slicer::cs_nodes_ct& nodes) const { + const base_node& A=nodes[iA].pt; + const base_node& B=nodes[iB].pt; + scalar_type a,b,c; // a*x^2 + b*x + c = 0 + a = gmm::vect_norm2_sqr(B-A); + if (a < EPS) + return pt_bin.is_in(iA) ? 0. : 1./EPS; + b = 2*gmm::vect_sp(A-x0,B-A); + c = gmm::vect_norm2_sqr(A-x0)-R*R; + return slicer_volume::trinom(a,b,c); + } + + void + slicer_cylinder::test_point(const base_node& P, bool& in, bool& bound) const { + base_node N = P; + if (2 == N.size()) /* Add Z coorinate if mesh is 2D */ + N.push_back(0.0); + N = N-x0; + scalar_type axpos = gmm::vect_sp(d, N); + scalar_type dist2 = gmm::vect_norm2_sqr(N) - gmm::sqr(axpos); + bound = gmm::abs(dist2-R*R) < EPS; + in = dist2 < R*R; + } + + scalar_type + slicer_cylinder::edge_intersect(size_type iA, size_type iB, + const mesh_slicer::cs_nodes_ct& nodes) const { + base_node F=nodes[iA].pt; + base_node D=nodes[iB].pt-nodes[iA].pt; + if (2 == F.size()) { + F.push_back(0.0); + D.push_back(0.0); + } + F = F - x0; + scalar_type Fd = gmm::vect_sp(F,d); + scalar_type Dd = gmm::vect_sp(D,d); + scalar_type a = gmm::vect_norm2_sqr(D) - gmm::sqr(Dd); + if (a < EPS) + return pt_bin.is_in(iA) ? 0. : 1./EPS; + assert(a> -EPS); + scalar_type b = 2*(gmm::vect_sp(F,D) - Fd*Dd); + scalar_type c = gmm::vect_norm2_sqr(F) - gmm::sqr(Fd) - gmm::sqr(R); + return slicer_volume::trinom(a,b,c); + } }