[Getfem-commits] (no subject)

[Konstantinos Poulios]

branch: master
commit 45c3f1baccf8ff8e804a2d14498532ef23fa1609
Author: Konstantinos Poulios <address@hidden>
Date:   Fri Nov 9 17:29:49 2018 +0100

    Add a ball shell structured mesh generator (tested in 2D and 3D)
---
 src/getfem/getfem_regular_meshes.h |  22 +++++++
 src/getfem_import.cc               |   2 +
 src/getfem_regular_meshes.cc       | 115 +++++++++++++++++++++++++++++++++++++
 3 files changed, 139 insertions(+)

diff --git a/src/getfem/getfem_regular_meshes.h 
b/src/getfem/getfem_regular_meshes.h
index 88270cf..eb007ee 100644
--- a/src/getfem/getfem_regular_meshes.h
+++ b/src/getfem/getfem_regular_meshes.h
@@ -168,6 +168,28 @@ namespace getfem
   */
   void regular_ball_mesh(mesh& m, const std::string &st);
 
+  /**
+     Build a regular mesh on a ball shell, parametrized by the string st.
+     The format of st is similar to getfem::regular_mesh.
+     @see getfem::import_mesh.
+     All parameters except the geometric transformation GT are optional.
+     Here, parameter NSUBDIV has to be a vector of size 2 that holds the
+     number of subdivisions in the circumferential and radial direction
+     of the mesh.
+     SIZES is a two element vector that provides the ball radius and shell
+     thickness (default radius is equal to 1 and default thickness is
+     ewual to 0.5).
+     If NOISED=1 the nodes in the interior of the individual sub-regions
+     of the mesh are randomly "shaken" (default value NOISED=0).
+     An additional integer paramater called SYMMETRIES receiving the
+     values 1, 2 or 3 (in three dimensions) permits to respectively request
+     one half, one quarter or one eighth of the ball to be meshed
+     (default value SYMMETRIES=0).
+
+     @param m the output mesh.
+  */
+  void regular_ball_shell_mesh(mesh& m, const std::string &st);
+
 }  /* end of namespace getfem.                                             */
 
 
diff --git a/src/getfem_import.cc b/src/getfem_import.cc
index 45ae18e..bdfc853 100644
--- a/src/getfem_import.cc
+++ b/src/getfem_import.cc
@@ -1356,6 +1356,8 @@ namespace getfem {
         { regular_mesh(m, filename); return; }
       else if (bgeot::casecmp(format,"structured_ball")==0)
         { regular_ball_mesh(m, filename); return; }
+      else if (bgeot::casecmp(format,"structured_ball_shell")==0)
+        { regular_ball_shell_mesh(m, filename); return; }
 
       std::ifstream f(filename.c_str());
       GMM_ASSERT1(f.good(), "can't open file " << filename);
diff --git a/src/getfem_regular_meshes.cc b/src/getfem_regular_meshes.cc
index d94b8ad..c319277 100644
--- a/src/getfem_regular_meshes.cc
+++ b/src/getfem_regular_meshes.cc
@@ -509,5 +509,120 @@ namespace getfem
     m.translation(org);
 
   }
+  void regular_ball_shell_mesh(mesh &m, const std::string &st) {
+    std::stringstream s(st);
+    bgeot::md_param PARAM;
+    PARAM.read_param_file(s);
+
+    std::string GT = PARAM.string_value("GT");
+    GMM_ASSERT1(!GT.empty(), "regular ball mesh : you have at least to "
+                "specify the geometric transformation");
+    bgeot::pgeometric_trans pgt = bgeot::geometric_trans_descriptor(GT);
+
+    size_type N = pgt->dim();
+    base_small_vector org(N);
+    const auto &o = PARAM.array_value("ORG");
+    if (o.size() > 0) {
+      GMM_ASSERT1(o.size() == N,
+                  "ORG parameter should be an array of size " << N);
+      for (size_type i = 0; i < N; ++i) {
+        GMM_ASSERT1(o[i].type_of_param() == bgeot::md_param::REAL_VALUE,
+                    "ORG should be a real array");
+        org[i] = o[i].real();
+      }
+    }
+
+    // "NOISED" applies only to the interior of all merged sub-meshes
+    bool noised = (PARAM.int_value("NOISED") != 0);
+
+    size_type nsubdiv0(3), nsubdiv1(2);
+    const auto &ns = PARAM.array_value("NSUBDIV");
+    if (ns.size() > 0) {
+      GMM_ASSERT1(ns.size() == 2,
+                  "NSUBDIV parameter should be an array of size 2");
+      for (size_type i = 0; i < 2; ++i)
+        GMM_ASSERT1(ns[i].type_of_param() == bgeot::md_param::REAL_VALUE,
+                    "NSUBDIV should be an integer array");
+      nsubdiv0 = size_type(ns[0].real()+0.5);
+      nsubdiv1 = size_type(ns[1].real()+0.5);
+    }
+
+    scalar_type radius(1), thickness(0.5);
+    const auto &si = PARAM.array_value("SIZES");
+    if (si.size() > 0) {
+      GMM_ASSERT1(si.size() == 2,
+                  "SIZES parameter should be an array of size 2");
+      GMM_ASSERT1(si[0].type_of_param() == bgeot::md_param::REAL_VALUE &&
+                  si[1].type_of_param() == bgeot::md_param::REAL_VALUE ,
+                  "SIZES should be a real array");
+      radius = si[0].real();
+      thickness = si[1].real();
+    }
+
+    std::vector<size_type> nsubdiv(N, nsubdiv0);
+    nsubdiv[N-1] = nsubdiv1;
+
+    mesh m0;
+    regular_unit_mesh(m0, nsubdiv, pgt, noised);
+
+    std::vector<base_node> pts(m0.points().card(), base_node(N));
+    size_type j(0);
+    for (dal::bv_visitor pt(m0.points().index()); !pt.finished(); ++pt, ++j) {
+      pts[j] = m0.points()[pt];
+      scalar_type l(1);
+      for (size_type k=0; k < N-1; ++k) {
+        pts[j][k] = tan(pts[j][k]*M_PI_4);
+        l += pts[j][k]*pts[j][k];
+      }
+      l = sqrt(l);
+      scalar_type r(radius-thickness+thickness*pts[j][N-1]);
+      for (size_type k=0; k < N-1; ++k)
+        pts[j][k] *= r/l;
+      pts[j][N-1] = r/l;
+    }
+
+    j = size_type(0);
+    for (dal::bv_visitor pt(m0.points().index()); !pt.finished(); ++pt, ++j)
+      m0.points()[pt] = pts[j];
+    m0.points().resort();
+
+    base_matrix M(N,N);
+    for (size_type i=0; i < N-1; ++i)
+      M(i,i+1) = scalar_type(1);
+    M(N-1,0) = scalar_type(1);
+
+    for (size_type i = 0; i < N; ++i) {
+      m0.transformation(M);
+      for (dal::bv_visitor cv(m0.convex_index()); !cv.finished(); ++cv)
+            m.add_convex_by_points(m0.trans_of_convex(cv),
+                               m0.points_of_convex(cv).begin());
+    }
+
+    size_type symmetries(PARAM.int_value("SYMMETRIES"));
+    symmetries = std::min(symmetries,N);
+
+    for (size_type sym=0; sym < N-symmetries; ++sym) {
+      size_type sym0 = (sym+1) % N;
+      if (sym0 != 0) {
+        gmm::clear(M);
+        M(sym,sym0) = scalar_type(-1);
+        M(sym0,sym) = scalar_type(1);
+        for (size_type i=0; i < N; ++i)
+          if (i != sym && i != sym0) M(i,i) = scalar_type(1);
+      } else {
+        base_matrix M1(M), M2(M);
+        gmm::mult(M1,M2,M);
+      }
+      m0.copy_from(m);
+      m0.transformation(M);
+      for (dal::bv_visitor cv(m0.convex_index()); !cv.finished(); ++cv)
+            m.add_convex_by_points(m0.trans_of_convex(cv),
+                               m0.points_of_convex(cv).begin());
+    }
+
+    m.translation(org);
+
+  }
+
 
 }  /* end of namespace getfem.                                             */