Merge branch 'bndlayer_surfacegeom' into 'master'

Bndlayer surfacegeom

See merge request ngsolve/netgen!575

libsrc/meshing/python_mesh.cpp

@@ -1023,6 +1023,11 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
             return self.AddFaceDescriptor (fd);
           })
 
+    .def ("AddSingularity", [](Mesh & self, PointIndex pi, double factor)
+         {
+	   self[pi].Singularity(factor);
+         })
+
     .def ("AddPoints", [](Mesh & self, py::buffer b1)
           {
             static Timer timer("Mesh::AddPoints");
@@ -1683,14 +1688,14 @@ project_boundaries : Optional[str] = None
                   }), py::arg("mapping"))
     .def(NGSPickle<SurfaceGeometry>())
     .def("GenerateMesh", [](shared_ptr<SurfaceGeometry> geo,
-                            bool quads, int nx, int ny, bool flip_triangles, py::list py_bbbpts, py::list py_bbbnames, py::list py_hppts, py::list py_hpbnd)
+                            bool quads, int nx, int ny, bool flip_triangles, py::list py_bbbpts, py::list py_bbbnames, py::list py_hppnts, py::dict/*list*/ py_hpbnd, py::dict py_layers)
            {
              if (py::len(py_bbbpts) != py::len(py_bbbnames))
                throw Exception("In SurfaceGeometry::GenerateMesh bbbpts and bbbnames do not have same lengths.");
              Array<Point<3>> bbbpts(py::len(py_bbbpts));
              Array<string> bbbname(py::len(py_bbbpts));
-             Array<Point<3>> hppts(py::len(py_hppts));
+             Array<Point<3>> hppnts(py::len(py_hppnts));
-             Array<float> hpptsfac(py::len(py_hppts));
+             Array<float> hppntsfac(py::len(py_hppnts));
              Array<string> hpbnd(py::len(py_hpbnd));
              Array<float> hpbndfac(py::len(py_hpbnd));
              for(int i = 0; i<py::len(py_bbbpts);i++)
@@ -1699,30 +1704,82 @@ project_boundaries : Optional[str] = None
                    bbbpts[i] = Point<3>(py::extract<double>(pnt[0])(),py::extract<double>(pnt[1])(),py::extract<double>(pnt[2])());
                    bbbname[i] = py::extract<string>(py_bbbnames[i])();
                  }
-             for(int i = 0; i<py::len(py_hppts);i++)
+             for(int i = 0; i<py::len(py_hppnts);i++)
 		 {
-                   py::tuple pnt = py::extract<py::tuple>(py_hppts[i])();
+                   py::tuple pnt = py::extract<py::tuple>(py_hppnts[i])();
-                   hppts[i] = Point<3>(py::extract<double>(pnt[0])(),py::extract<double>(pnt[1])(),py::extract<double>(pnt[2])());
+                   hppnts[i] = Point<3>(py::extract<double>(pnt[0])(),py::extract<double>(pnt[1])(),py::extract<double>(pnt[2])());
-                   //hpptsfac[i] = py::len(pnt) > 3 ? py::extract<double>(pnt[3])() : 0.0;
+                   hppntsfac[i] = py::extract<double>(pnt[3])();
-                   hpptsfac[i] = py::extract<double>(pnt[3])();
                  }
 
-             for(int i = 0; i<py::len(py_hpbnd);i++)
+	     int ii=0;
+             for(auto val : py_hpbnd)
                {
-                   py::tuple bnd = py::extract<py::tuple>(py_hpbnd[i])();
+                 hpbnd[ii] = py::cast<string>(val.first);
-                   hpbnd[i] = py::extract<string>(bnd[0])();
+		 hpbndfac[ii] = py::cast<float>(val.second);
-                   hpbndfac[i] = py::extract<double>(bnd[1])();
+		 ii++;
 	       }
+
+             
+             Array<double> layer_thickness[4];
+             bool layer_quad = false;
+
+             for(auto val : py_layers)
+               {
+		 int index = -1;
+                 if (py::cast<string>(val.first) == "left") index = 0;
+                 else if (py::cast<string>(val.first) == "top") index = 3;
+                 else if (py::cast<string>(val.first) == "right") index = 2;
+                 else if (py::cast<string>(val.first) == "bottom") index = 1;
+		 else if (py::cast<string>(val.first) == "quads") layer_quad = py::cast<bool>(val.second);
+		 else throw Exception("Unknown parameter " + string(py::cast<string>(val.first)));
+		 if (index < 0) continue;
+
+		 auto list = py::cast<py::list>(val.second);
+		 layer_thickness[index] = Array<double>(py::len(list));
+		 for (size_t i = 0; i < py::len(list); i++)
+		   layer_thickness[index][i] = py::cast<double>(list[i]);
+               }
+                   
              auto mesh = make_shared<Mesh>();
              SetGlobalMesh (mesh);
              mesh->SetGeometry(geo);
 	     ng_geometry = geo;
-             auto result = geo->GenerateStructuredMesh (mesh, quads, nx, ny, flip_triangles, bbbpts, bbbname, hppts, hpptsfac, hpbnd, hpbndfac);
+             auto result = geo->GenerateStructuredMesh (mesh, quads, nx, ny, flip_triangles, bbbpts, bbbname, hppnts, hppntsfac, hpbnd, hpbndfac, layer_thickness, layer_quad);
              if(result != 0)
                throw Exception("SurfaceGeometry: Meshing failed!");
              return mesh;
-           }, py::arg("quads")=true, py::arg("nx")=10, py::arg("ny")=10, py::arg("flip_triangles")=false, py::arg("bbbpts")=py::list(), py::arg("bbbnames")=py::list(), py::arg("hppts")=py::list(), py::arg("hpbnd")=py::list())
+           }, py::arg("quads")=true, py::arg("nx")=10, py::arg("ny")=10, py::arg("flip_triangles")=false, py::arg("bbbpts")=py::list(), py::arg("bbbnames")=py::list(), py::arg("hppnts")=py::list(), py::arg("hpbnd")=py::dict(), py::arg("boundarylayer")=py::dict());/*, R"raw_string(
-      ;
+      Generate a structured 2D surface mesh
+
+    Parameters:
+    
+    quads : bool
+      If True, a quadrilateral mesh is generated. If False, the quads are split to triangles.
+
+    nx : int
+      Number of cells in x-direction.
+
+    ny : int
+      Number of cells in y-direction.
+
+    flip_triangles : bool
+      If set to True together with quads=False the quads are cut the other way round
+
+    bbbpts : list
+      List of points which should be handled as BBBND and are named with bbbnames. The mesh must be constructed in such a way that the bbbpts coincide with generated points.
+
+    bbbnames : list
+      List of bbbnd names as strings. Size must coincide with size of bbbpts.
+
+    hppnts : list
+      If not None it expects a list of the form [ (px1,py1,pz1, hpref1), (px2,py2,pz2, hpref2), ... ] where px,py,pz are the point coordinates which have to be resolved in the mesh and hpref the refinement factor.
+
+    hpbnd : dict
+      If not None it expects a dictionary of the form {"boundaryname" : hpref } where boundaryname in [left, right, top, bottom] and hpref the refinement factor.
+
+    boundarylayer : dict
+      If not None it expects a dictionary of the form { "boundaryname" : [t1,...,tn], "quads" : False } where ti denote the thickness of layer i. The number of layers are included in nx/ny. After the layers are placed the remaining number of cells are used to divide the remaining grid uniformly. If quads are set to True quadrilaterals are used inside the boundarylayer. If set False the value of "quads" of the function call is used.
+      )raw_string");*/
     ;
 
     py::class_<ClearSolutionClass> (m, "ClearSolutionClass")

libsrc/meshing/surfacegeom.cpp

@@ -197,6 +197,7 @@ namespace netgen
     return false;
   }
 
+  
   void  SurfaceGeometry :: PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint, int surfi1, int surfi2, const EdgePointGeomInfo & ap1, const EdgePointGeomInfo & ap2, Point<3> & newp, EdgePointGeomInfo & newgi) const
   {
     newgi.u = ap1.u+secpoint*(ap2.u-ap1.u);
@@ -208,6 +209,32 @@ namespace netgen
     newp = Point<3>(func(Point<2>(newgi.u, newgi.v)));
   }
 
+    void CheckForBBBPnt(const Array<Point<3>>& bbbpts, const Point<3>& pnt, Array<bool>& found, Array<PointIndex>& indbbbpts, const Array<PointIndex>& pids)
+  {
+    for (int k = 0; k < bbbpts.Size(); k++)
+      {
+        auto diff = pnt - bbbpts[k];
+        if(diff.Length2() < 1e-14)
+          {
+            found[k] = true;
+            indbbbpts[k] = pids[pids.Size()-1];
+          }
+      }
+  }
+
+  void CheckForSingularity(const Array<Point<3>>& hppoints, const Point<3>& pnt, const Array<float>& hpptsfac, shared_ptr<Mesh> & mesh, const Array<PointIndex>& pids)
+  {
+    for (int k = 0; k < hppoints.Size(); k++)
+      {
+        auto diff = pnt - hppoints[k];
+        if(diff.Length2() < 1e-14)
+          {
+            (*mesh)[pids[pids.Size()-1]].Singularity(hpptsfac[k]);
+          }
+        
+      }
+  }
+  
   void SurfaceGeometry :: PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
                                        int surfi, 
                                        const PointGeomInfo & gi1, 
@@ -223,7 +250,7 @@ namespace netgen
     //ProjectPointGI(surfi, newp, newgi);
   }
 
-  int SurfaceGeometry :: GenerateStructuredMesh(shared_ptr<Mesh> & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array<Point<3>>& bbbpts, const Array<string>& bbbnames, const Array<Point<3>>& hppoints, const Array<float>& hpptsfac, const Array<string>& hpbnd, const Array<float>& hpbndfac)
+  int SurfaceGeometry :: GenerateStructuredMesh(shared_ptr<Mesh> & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array<Point<3>>& bbbpts, const Array<string>& bbbnames, const Array<Point<3>>& hppoints, const Array<float>& hpptsfac, const Array<string>& hpbnd, const Array<float>& hpbndfac, Array<double> layer_thickness[4], bool layer_quad)
   {
     mesh->SetDimension(3);
 
@@ -231,40 +258,86 @@ namespace netgen
     found = false;
     Array<PointIndex> indbbbpts(bbbpts.Size());
 
+    int numx = nx;
+    int numy = ny;
 
-    Array<PointIndex> pids;
+    size_t total_layer_el[4] = {layer_thickness[0].Size(), layer_thickness[1].Size(), layer_thickness[2].Size(), layer_thickness[3].Size()};
-    Array<PointGeomInfo> pgis;
+    
-    for(int i=0; i <= ny; i++)
+    double interior_x = 1.0;
-      for(int j=0; j <= nx; j++)
+    double interior_y = 1.0;
+    for(double scale : layer_thickness[0])
+      interior_x -= scale;
+    for(double scale : layer_thickness[2])
+      interior_x -= scale;
+    for(double scale : layer_thickness[1])
+      interior_y -= scale;
+    for(double scale : layer_thickness[3])
+      interior_y -= scale;
+
+    auto AddPoint = [&] (double offsetx, double offsety, Array<PointIndex> & pids, Array<PointGeomInfo> & pgis)
     {
       PointGeomInfo pgi;
       pgi.trignum = -1;
-          pgi.u = double(j)/nx;
+      pgi.u = offsetx;
-          pgi.v = double(i)/ny;
+      pgi.v = offsety;
+
       
       Point<3> pnt = Point<3>(func(Point<2>(pgi.u,pgi.v)));
       pids.Append(mesh->AddPoint(pnt));
       pgis.Append(pgi);
       
-          for (int k = 0; k < bbbpts.Size(); k++)
+      CheckForBBBPnt(bbbpts, pnt, found, indbbbpts, pids);
+      CheckForSingularity(hppoints, pnt, hpptsfac, mesh, pids);
+    };
+
+    auto InternalLoop = [&] (double offsety, Array<PointIndex> & pids, Array<PointGeomInfo> & pgis)
     {
-              auto diff = pnt - bbbpts[k];
+      int j = 0;
-              if(diff.Length2() < 1e-14)
+      double offsetx = 0.0;
+
+      for(int l=0; l < layer_thickness[0].Size(); l++,j++)
         {
-                  found[k] = true;
+	  AddPoint(offsetx+layer_thickness[0][l]*double(j-l), offsety, pids, pgis);
-                  indbbbpts[k] = pids[pids.Size()-1];
+          offsetx += layer_thickness[0][l];
-                }
         }
       
-          for (int k = 0; k < hppoints.Size(); k++)
+      
+      for(;j <= nx-total_layer_el[2]; j++)
+        AddPoint(offsetx + interior_x*double(j-total_layer_el[0])/(nx-total_layer_el[0]-total_layer_el[2]), offsety, pids, pgis);
+      offsetx += interior_x;
+            
+      int startj = j;
+      for(int l=0; l < layer_thickness[2].Size(); l++, j++)
         {
-              auto diff = pnt - hppoints[k];
+	  AddPoint(offsetx+layer_thickness[2][layer_thickness[2].Size()-1-l]*double(j-startj-l+1), offsety, pids, pgis);
-              if(diff.Length2() < 1e-14)
+
+          offsetx += layer_thickness[2][layer_thickness[2].Size()-1-l];
+        }
+    };
+    
+    Array<PointIndex> pids;
+    Array<PointGeomInfo> pgis;
+    
+    int i = 0;
+    double offsety = 0.0;
+
+    for(int k=0; k < layer_thickness[1].Size(); k++,i++)
       {
-                  (*mesh)[pids[pids.Size()-1]].Singularity(hpptsfac[k]);
+	InternalLoop(offsety, pids, pgis);
+	offsety += layer_thickness[1][k];
       }
 
+    for(; i <= ny-total_layer_el[3]; i++)
+      {
+        InternalLoop(offsety, pids, pgis);        
+	offsety +=  interior_y/(ny-total_layer_el[1]-total_layer_el[3]);
       }
+    offsety -=  interior_y/(ny-total_layer_el[1]-total_layer_el[3]);
+
+    for(int k=0; k < layer_thickness[3].Size(); k++,i++)
+      {
+	offsety += layer_thickness[3][layer_thickness[3].Size()-1-k];
+	InternalLoop(offsety, pids, pgis); 
       }
 
     for (bool f : found)
@@ -279,14 +352,14 @@ namespace netgen
     mesh->AddFaceDescriptor(fd);
 
 
-    for(int i=0; i < ny; i++)
+    for(int i=0; i < numy; i++)
       {
-        for(int j=0; j < nx; j++)
+        for(int j=0; j < numx; j++)
           {
-            int base = i * (nx+1) + j;
+            int base = i * (numx+1) + j;
-            if (quads)
+            if (quads || (layer_quad && i < total_layer_el[1]) || (layer_quad && i > numy-1-total_layer_el[3]) || (layer_quad && j < total_layer_el[0]) || (layer_quad && j > numx-1-total_layer_el[2]) )
               {
-                int pnum[4] = {base,base+1,base+nx+2,base+nx+1};
+                int pnum[4] = {base,base+1,base+numx+2,base+numx+1};
                 Element2d el = Element2d(QUAD);
                 for (int i = 0; i < 4; i++)
                   {
@@ -305,19 +378,19 @@ namespace netgen
                   {
                     pnum1[0] = base;
                     pnum1[1] = base+1;
-                    pnum1[2] = base+nx+2;
+                    pnum1[2] = base+numx+2;
                     pnum2[0] = base;
-                    pnum2[1] = base+nx+2;
+                    pnum2[1] = base+numx+2;
-                    pnum2[2] = base+nx+1;
+                    pnum2[2] = base+numx+1;
                   }
                 else
                   {
                     pnum1[0] = base;
                     pnum1[1] = base+1;
-                    pnum1[2] = base+nx+1;
+                    pnum1[2] = base+numx+1;
                     pnum2[0] = base+1;
-                    pnum2[1] = base+nx+2;
+                    pnum2[1] = base+numx+2;
-                    pnum2[2] = base+nx+1;
+                    pnum2[2] = base+numx+1;
                   }
 
                 Element2d el = Element2d(TRIG);
@@ -357,7 +430,7 @@ namespace netgen
       }
     // needed for codim2 in 3d
     seg.edgenr = 1;
-    for(int i=0; i < nx; i++)
+    for(int i=0; i < numx; i++)
       {
         seg[0] = pids[i];
         seg[1] = pids[i+1];
@@ -388,18 +461,18 @@ namespace netgen
           }
       }
 
-    for(int i=0; i<ny; i++)
+    for(int i=0; i<numy; i++)
       {
-        seg[0] = pids[i*(nx+1)+nx];
+        seg[0] = pids[i*(numx+1)+numx];
-        seg[1] = pids[(i+1)*(nx+1)+nx];
+        seg[1] = pids[(i+1)*(numx+1)+numx];
 
-        seg.geominfo[0] = pgis[i*(nx+1)+nx];
+        seg.geominfo[0] = pgis[i*(numx+1)+numx];
-        seg.geominfo[1] = pgis[(i+1)*(nx+1)+nx];
+        seg.geominfo[1] = pgis[(i+1)*(numx+1)+numx];
-        seg.epgeominfo[0].u = pgis[i*(nx+1)+nx].u;
+        seg.epgeominfo[0].u = pgis[i*(numx+1)+numx].u;
-        seg.epgeominfo[0].v = pgis[i*(nx+1)+nx].v;
+        seg.epgeominfo[0].v = pgis[i*(numx+1)+numx].v;
         seg.epgeominfo[0].edgenr = seg.edgenr;
-        seg.epgeominfo[1].u = pgis[(i+1)*(nx+1)+nx].u;
+        seg.epgeominfo[1].u = pgis[(i+1)*(numx+1)+numx].u;
-        seg.epgeominfo[1].v = pgis[(i+1)*(nx+1)+nx].v;
+        seg.epgeominfo[1].v = pgis[(i+1)*(numx+1)+numx].v;
         seg.epgeominfo[1].edgenr = seg.edgenr;
 
         mesh->AddSegment(seg);
@@ -419,18 +492,18 @@ namespace netgen
           }
       }
 
-    for(int i=0; i<nx; i++)
+    for(int i=0; i<numx; i++)
       {
-        seg[0] = pids[ny*(nx+1)+i+1];
+        seg[0] = pids[numy*(numx+1)+i+1];
-        seg[1] = pids[ny*(nx+1)+i];
+        seg[1] = pids[numy*(numx+1)+i];
 
-        seg.geominfo[0] = pgis[ny*(nx+1)+i+1];
+        seg.geominfo[0] = pgis[numy*(numx+1)+i+1];
-        seg.geominfo[1] = pgis[ny*(nx+1)+i];
+        seg.geominfo[1] = pgis[numy*(numx+1)+i];
-        seg.epgeominfo[0].u = pgis[ny*(nx+1)+i+1].u;
+        seg.epgeominfo[0].u = pgis[numy*(numx+1)+i+1].u;
-        seg.epgeominfo[0].v = pgis[ny*(nx+1)+i+1].v;
+        seg.epgeominfo[0].v = pgis[numy*(numx+1)+i+1].v;
         seg.epgeominfo[0].edgenr = seg.edgenr;
-        seg.epgeominfo[1].u = pgis[ny*(nx+1)+i].u;
+        seg.epgeominfo[1].u = pgis[numy*(numx+1)+i].u;
-        seg.epgeominfo[1].v = pgis[ny*(nx+1)+i].v;
+        seg.epgeominfo[1].v = pgis[numy*(numx+1)+i].v;
         seg.epgeominfo[1].edgenr = seg.edgenr;
         
         mesh->AddSegment(seg);
@@ -450,18 +523,18 @@ namespace netgen
       }
 
 
-    for(int i=0; i<ny; i++)
+    for(int i=0; i<numy; i++)
       {
-        seg[0] = pids[(i+1)*(nx+1)];
+        seg[0] = pids[(i+1)*(numx+1)];
-        seg[1] = pids[i*(nx+1)];
+        seg[1] = pids[i*(numx+1)];
 
-        seg.geominfo[0] = pgis[(i+1)*(nx+1)];
+        seg.geominfo[0] = pgis[(i+1)*(numx+1)];
-        seg.geominfo[1] = pgis[i*(nx+1)];
+        seg.geominfo[1] = pgis[i*(numx+1)];
-        seg.epgeominfo[0].u = pgis[(i+1)*(nx+1)].u;
+        seg.epgeominfo[0].u = pgis[(i+1)*(numx+1)].u;
-        seg.epgeominfo[0].v = pgis[(i+1)*(nx+1)].v;
+        seg.epgeominfo[0].v = pgis[(i+1)*(numx+1)].v;
         seg.epgeominfo[0].edgenr = seg.edgenr;
-        seg.epgeominfo[1].u = pgis[i*(nx+1)].u;
+        seg.epgeominfo[1].u = pgis[i*(numx+1)].u;
-        seg.epgeominfo[1].v = pgis[i*(nx+1)].v;
+        seg.epgeominfo[1].v = pgis[i*(numx+1)].v;
         seg.epgeominfo[1].edgenr = seg.edgenr;
 
         mesh->AddSegment(seg);

libsrc/meshing/surfacegeom.hpp

@@ -62,7 +62,7 @@ namespace netgen
                       const PointGeomInfo & gi2,
                       Point<3> & newp, PointGeomInfo & newgi) const override;
 
-    int GenerateStructuredMesh(shared_ptr<Mesh> & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array<Point<3>>& bbbpts, const Array<string>& bbbnames, const Array<Point<3>>& hppoints, const Array<float>& hpptsfac, const Array<string>& hpbnd, const Array<float>& hpbndfac);
+    int GenerateStructuredMesh(shared_ptr<Mesh> & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array<Point<3>>& bbbpts, const Array<string>& bbbnames, const Array<Point<3>>& hppoints, const Array<float>& hpptsfac, const Array<string>& hpbnd, const Array<float>& hpbndfac, Array<double> layer_thickness[4], bool layer_quad);
 
   };