[Getfem-commits] (no subject)

[Konstantinos Poulios]

branch: fix-project-into-element
commit 3e2655bc74191572e01ec7ecc3d949404a82a399
Author: Konstantinos Poulios <address@hidden>
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__  */