[Getfem-commits] r4777 - in /trunk/getfem: interface/src/gf_mesh.cc src/getfem/getfem_mesh.h src/getfem_mesh.cc

[logari81]

Author: logari81
Date: Tue Oct 14 09:40:09 2014
New Revision: 4777

URL: http://svn.gna.org/viewcvs/getfem?rev=4777&view=rev
Log:
enable higher order mesh extrusion

Modified:
    trunk/getfem/interface/src/gf_mesh.cc
    trunk/getfem/src/getfem/getfem_mesh.h
    trunk/getfem/src/getfem_mesh.cc

Modified: trunk/getfem/interface/src/gf_mesh.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/gf_mesh.cc?rev=4777&r1=4776&r2=4777&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_mesh.cc       (original)
+++ trunk/getfem/interface/src/gf_mesh.cc       Tue Oct 14 09:40:09 2014
@@ -243,8 +243,10 @@
 prismatic_mesh(getfem::mesh *dest_mesh, getfemint::mexargs_in &in)
 {
   const getfem::mesh *src_mesh = in.pop().to_const_mesh();
-  unsigned nblay = in.pop().to_integer(1,2500000);
-  getfem::extrude(*src_mesh, *dest_mesh, nblay);
+  size_type nblay = in.pop().to_integer(1,2500000);
+  short_type degree(1);
+  if (in.remaining()) degree = in.pop().to_integer(1,2500000);
+  getfem::extrude(*src_mesh, *dest_mesh, nblay, degree);
 }
 
 static void
@@ -400,13 +402,16 @@
        );
 
 
-    /address@hidden M = ('prismatic', @tmesh m, @int nl)
+    /address@hidden M = ('prismatic', @tmesh m, @int nl[, @int degree])
       Extrude a prismatic @tmesh `M` from a @tmesh `m`.
 
-      In the additional dimension there are `nl` layers of elements built
-      from ``0`` to address@hidden/
-    sub_command
-      ("prismatic", 2, 2, 0, 1,
+      In the additional dimension there are `nl` layers of elements
+      distributed from ``0`` to ``1``.
+      If the optional parameter `degree` is provided with a value greater
+      than the default value of ``1``, a non-linear transformation of
+      corresponding degree is considered in the extrusion address@hidden/
+    sub_command
+      ("prismatic", 2, 3, 0, 1,
        prismatic_mesh(pmesh, in);
        );
 

Modified: trunk/getfem/src/getfem/getfem_mesh.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_mesh.h?rev=4777&r1=4776&r2=4777&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_mesh.h       (original)
+++ trunk/getfem/src/getfem/getfem_mesh.h       Tue Oct 14 09:40:09 2014
@@ -595,7 +595,8 @@
   /**
    * build a N+1 dimensions mesh from a N-dimensions mesh by extrusion.
    */
-  void APIDECL extrude(const mesh& in, mesh& out, unsigned nb_layers);
+  void APIDECL extrude(const mesh& in, mesh& out, size_type nb_layers,
+                       short_type degree=short_type(1));
 
  template<class ITER>
     size_type mesh::add_convex_by_points(bgeot::pgeometric_trans pgt,

Modified: trunk/getfem/src/getfem_mesh.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_mesh.cc?rev=4777&r1=4776&r2=4777&view=diff
==============================================================================
--- trunk/getfem/src/getfem_mesh.cc     (original)
+++ trunk/getfem/src/getfem_mesh.cc     Tue Oct 14 09:40:09 2014
@@ -244,7 +244,7 @@
     Pmin.clear(); Pmax.clear();
     for (dal::bv_visitor i(pts.index()); !i.finished(); ++i) {
       if (is_first) { Pmin = Pmax = pts[i]; is_first = false; }
-      else for (unsigned j=0; j < dim(); ++j) {
+      else for (dim_type j=0; j < dim(); ++j) {
         Pmin[j] = std::min(Pmin[j], pts[i][j]);
         Pmax[j] = std::max(Pmax[j], pts[i][j]);
       }
@@ -786,33 +786,34 @@
     return mrr;
   }
 
-  void extrude(const mesh& in, mesh& out, unsigned nb_layers) {
-    unsigned dim = in.dim();
+  void extrude(const mesh& in, mesh& out, size_type nb_layers, short_type 
degree) {
+    dim_type dim = in.dim();
     base_node pt(dim+1);
     out.clear();
     size_type nbpt = in.points().index().last()+1;
     GMM_ASSERT1(nbpt == in.points().index().card(),
                 "please optimize the mesh before using "
                 "it as a base for prismatic mesh");
+    size_type nb_layers_total = nb_layers * degree;
     for (size_type i = 0; i < nbpt; ++i) {
       std::copy(in.points()[i].begin(), in.points()[i].end(),pt.begin());
       pt[dim] = 0.0;
-      for (size_type j = 0; j <= nb_layers; ++j, pt[dim] += 1.0 / nb_layers)
+      for (size_type j = 0; j <= nb_layers_total;
+           ++j, pt[dim] += scalar_type(1) / scalar_type(nb_layers_total))
         out.add_point(pt);
     }
 
     std::vector<size_type> tab;
     for (dal::bv_visitor cv(in.convex_index()); !cv.finished(); ++cv) {
       size_type nbp = in.nb_points_of_convex(cv);
-      tab.resize(2*nbp);
+      tab.resize((degree+1)*nbp);
       for (size_type j = 0; j < nb_layers; ++j) {
-        for (size_type k = 0; k < nbp; ++k)
-          tab[k] = (nb_layers+1)*in.ind_points_of_convex(cv)[k] + j;
-        for (size_type k = 0; k < nbp; ++k)
-          tab[k+nbp] = (nb_layers+1)*in.ind_points_of_convex(cv)[k] + j + 1;
+        for (short_type d = 0; d < degree+1; ++d)
+          for (size_type k = 0; k < nbp; ++k)
+            tab[k+nbp*d] = (nb_layers_total+1)*in.ind_points_of_convex(cv)[k] 
+ j*degree + d;
         bgeot::pgeometric_trans pgt =
           bgeot::product_geotrans(in.trans_of_convex(cv),
-                                  bgeot::simplex_geotrans(1,1));
+                                  bgeot::simplex_geotrans(1,degree));
         out.add_convex(pgt, tab.begin());
       }
     }