branch: fix-project-into-element commit 3e2655bc74191572e01ec7ecc3d949404a82a399 Author: Konstantinos Poulios <logar...@gmail.com> Date: Sat Sep 8 06:41:02 2018 +0200
New implementation of project_into_convex as a geotrans class function --- src/bgeot_convex_ref.cc | 66 ++++++++++++++++++++++++++++++++++++++ src/bgeot_torus.cc | 6 +++- src/getfem/bgeot_convex_ref.h | 7 ++++ src/getfem/bgeot_geometric_trans.h | 2 ++ src/getfem/bgeot_torus.h | 3 +- 5 files changed, 82 insertions(+), 2 deletions(-) diff --git a/src/bgeot_convex_ref.cc b/src/bgeot_convex_ref.cc index 1c94182..dbbbf28 100644 --- a/src/bgeot_convex_ref.cc +++ b/src/bgeot_convex_ref.cc @@ -280,6 +280,22 @@ namespace bgeot { for (; it != ite; e += *it, ++it) {}; return e / sqrt(scalar_type(pt.size())); } + + void project_into(base_node &pt) const { + if (auto_basic) { + GMM_ASSERT1(pt.size() == cvs->dim(), + "K_simplex_of_ref_::project_into: Dimensions mismatch"); + scalar_type sum_coordinates = 0.0; + for (const auto &coord : pt) sum_coordinates += coord; + if (sum_coordinates > 1.0) gmm::scale(pt, 1.0 / sum_coordinates); + for (auto &coord : pt) { + if (coord < 0.0) coord = 0.0; + if (coord > 1.0) coord = 1.0; + } + } else + basic_convex_ref_->project_into(pt); + } + K_simplex_of_ref_(dim_type NN, short_type KK) : convex_of_reference(simplex_structure(NN, KK), (KK == 1) || (NN == 0)) { @@ -443,6 +459,19 @@ namespace bgeot { return r; } + void project_into(base_node &pt) const { + if (auto_basic) { + GMM_ASSERT1(pt.size() == 3, "Dimensions mismatch"); + if (pt[2] < .0) pt[2] = 0.; + for (short_type f = 1; f < 5; ++f) { + scalar_type reldist = gmm::vect_sp(normals_[f], pt)*sqrt(2.); + if (reldist > 1.) + gmm::scale(pt, 1./reldist); + } + } else + basic_convex_ref_->project_into(pt); + } + pyramid_QK_of_ref_(dim_type k) : convex_of_reference(pyramid_QK_structure(k), k == 1) { GMM_ASSERT1(k == 1 || k == 2, "Sorry exist only in degree 1 or 2, not " << k); @@ -651,6 +680,21 @@ namespace bgeot { else return cvr2->is_in_face(short_type(f - cvr1->structure()->nb_faces()), pt2); } + void project_into(base_node &pt) const { + if (auto_basic) { + GMM_ASSERT1(pt.size() == cvs->dim(), "Dimensions mismatch"); + dim_type n1 = cvr1->structure()->dim(), n2 = cvr2->structure()->dim(); + base_node pt1(n1), pt2(n2); + std::copy(pt.begin(), pt.begin()+n1, pt1.begin()); + std::copy(pt.begin()+n1, pt.end(), pt2.begin()); + cvr1->project_into(pt1); + cvr2->project_into(pt2); + std::copy(pt1.begin(), pt1.end(), pt.begin()); + std::copy(pt2.begin(), pt2.end(), pt.begin()+n1); + } else + basic_convex_ref_->project_into(pt); + } + product_ref_(pconvex_ref a, pconvex_ref b) : convex_of_reference( convex_direct_product(*a, *b).structure(), @@ -714,6 +758,7 @@ namespace bgeot { /* equilateral ref convexes are used for estimation of convex quality */ class equilateral_simplex_of_ref_ : public convex_of_reference { + scalar_type r_inscr; public: scalar_type is_in(const base_node &pt) const { GMM_ASSERT1(pt.size() == cvs->dim(), "Dimension does not match"); @@ -733,9 +778,28 @@ namespace bgeot { : convex<base_node>::points().back()); return gmm::vect_sp(pt-x0, normals()[f]); } + + void project_into(base_node &pt) const { + dim_type N = cvs->dim(); + GMM_ASSERT1(pt.size() == N, "Dimension does not match"); + base_node G(N); G.fill(0.); + for (const base_node &x : convex<base_node>::points()) + G += x; + gmm::scale(G, scalar_type(1)/scalar_type(N+1)); + for (size_type f = 0; f < normals().size(); ++f) { + scalar_type r = gmm::vect_sp(pt-G, normals()[f]); + if (r > r_inscr) + pt = G + r_inscr/r*(pt-G); + } + + } + equilateral_simplex_of_ref_(size_type N) : convex_of_reference(simplex_structure(dim_type(N), 1), true) { + //https://math.stackexchange.com/questions/2739915/radius-of-inscribed-sphere-of-n-simplex + r_inscr = scalar_type(1)/sqrt(scalar_type(2*N)*scalar_type(N+1)); + pconvex_ref prev = equilateral_simplex_of_reference(dim_type(N-1)); convex<base_node>::points().resize(N+1); normals_.resize(N+1); @@ -785,6 +849,8 @@ namespace bgeot { { GMM_ASSERT1(false, "Information not available here"); } scalar_type is_in_face(short_type, const base_node &) const { GMM_ASSERT1(false, "Information not available here"); } + void project_into(base_node &) const + { GMM_ASSERT1(false, "Operation not available here"); } generic_dummy_(dim_type d, size_type n, short_type nf) : convex_of_reference(generic_dummy_structure(d, n, nf), true) diff --git a/src/bgeot_torus.cc b/src/bgeot_torus.cc index 9baf2f1..056f501 100644 --- a/src/bgeot_torus.cc +++ b/src/bgeot_torus.cc @@ -174,6 +174,10 @@ namespace bgeot{ GMM_ASSERT1(false, "Sorry, Hessian is not supported in axisymmetric transformation."); } + void torus_geom_trans::project_into_reference_convex(base_node &pt) const { + poriginal_trans_->project_into_reference_convex(pt); + } + torus_geom_trans::torus_geom_trans(pgeometric_trans poriginal_trans) : poriginal_trans_(poriginal_trans){ geometric_trans::is_lin = poriginal_trans_->is_linear(); @@ -210,4 +214,4 @@ namespace bgeot{ return pgt_torus != NULL; } -} \ No newline at end of file +} diff --git a/src/getfem/bgeot_convex_ref.h b/src/getfem/bgeot_convex_ref.h index 7eef34f..37b885f 100644 --- a/src/getfem/bgeot_convex_ref.h +++ b/src/getfem/bgeot_convex_ref.h @@ -108,6 +108,13 @@ namespace bgeot { * positive in the other side. */ virtual scalar_type is_in_face(short_type, const base_node &) const =0; + /** will project any given point lying outside the convex onto the convex + outer surface */ + virtual void project_into(base_node &pt) const { + GMM_ASSERT1(!auto_basic, "This method has to be overloaded in every " + "basic convex"); + basic_convex_ref_->project_into(pt); + } bool is_basic() const { return auto_basic; } /// return the normal vector for each face. const std::vector<base_small_vector> &normals() const diff --git a/src/getfem/bgeot_geometric_trans.h b/src/getfem/bgeot_geometric_trans.h index cc4a7a1..5159b0a 100644 --- a/src/getfem/bgeot_geometric_trans.h +++ b/src/getfem/bgeot_geometric_trans.h @@ -159,6 +159,8 @@ namespace bgeot { template<class CONT> base_node transform(const base_node &pt, const CONT &PTAB) const; base_node transform(const base_node &pt, const base_matrix &G) const; + virtual void project_into_reference_convex(base_node &pt) const + { cvr->project_into(pt); } size_type complexity() const { return complexity_; } virtual ~geometric_trans() { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Geometric transformation"); } diff --git a/src/getfem/bgeot_torus.h b/src/getfem/bgeot_torus.h index 1228792..91d90ae 100644 --- a/src/getfem/bgeot_torus.h +++ b/src/getfem/bgeot_torus.h @@ -57,6 +57,7 @@ struct torus_geom_trans : public geometric_trans{ (const bgeot::base_matrix &, const bgeot::base_matrix &, bgeot::base_matrix &) const; virtual void poly_vector_hess(const base_node &, bgeot::base_matrix &) const; + virtual void project_into_reference_convex(base_node &) const; torus_geom_trans(bgeot::pgeometric_trans poriginal_trans); @@ -76,4 +77,4 @@ bool is_torus_geom_trans(pgeometric_trans pgt); } -#endif /* BGEOT_TORUS_H__ */ \ No newline at end of file +#endif /* BGEOT_TORUS_H__ */