Merge branch 'close_edges' into 'master'

Close edges - Generate quads in narrow regions using geo._SetDomainTensorMeshing()

See merge request jschoeberl/netgen!321

libsrc/geom2d/genmesh2d.cpp

@@ -527,11 +527,16 @@ namespace netgen
 
 
           for (PointIndex pix = nextpi[c1], ix = 0; pix != c2; pix = nextpi[pix], ix++)
+          {
+            Point<3> px = (*mesh)[pix];
             for (PointIndex piy = nextpi[c2], iy = 0; piy != c3; piy = nextpi[piy], iy++)
               {
-                Point<3> p = (*mesh)[pix] + ( (*mesh)[piy] - (*mesh)[c2] );
-                pts[(nex+1)*(iy+1) + ix+1] = mesh -> AddPoint (p , 1, FIXEDPOINT);
+                double lam = Dist((*mesh)[piy],(*mesh)[c2]) / Dist((*mesh)[c3],(*mesh)[c2]);
+                auto pix1 = pts[(nex+1)*ney+ix+1];
+                auto pnew = px + lam*((*mesh)[pix1]-px);
+                pts[(nex+1)*(iy+1) + ix+1] = mesh -> AddPoint (pnew, 1, FIXEDPOINT);
               }
+          }
 
           for (int i = 0; i < ney; i++)
             for (int j = 0; j < nex; j++)

libsrc/geom2d/geometry2d.hpp

@@ -133,7 +133,7 @@ namespace netgen
     NgArray<char*> materials;
     NgArray<double> maxh;
     NgArray<bool> quadmeshing;
-    NgArray<bool> tensormeshing;
+    Array<bool> tensormeshing;
     NgArray<int> layer;
     NgArray<string*> bcnames;
     double elto0 = 1.0;
@@ -216,9 +216,20 @@ namespace netgen
     }
     bool GetDomainTensorMeshing ( int domnr ) 
     { 
-      if ( tensormeshing.Size() ) return tensormeshing[domnr-1]; 
+      if ( tensormeshing.Size()>=domnr ) return tensormeshing[domnr-1];
       else return false;
     }
+    void SetDomainTensorMeshing ( int domnr, bool tm )
+    {
+      if ( tensormeshing.Size()<domnr )
+      {
+        auto oldsize = tensormeshing.Size();
+        tensormeshing.SetSize(domnr);
+        for(auto i : IntRange(oldsize, domnr-1))
+          tensormeshing[i] = false;
+      }
+      tensormeshing[domnr-1] = tm;
+    }
     int GetDomainLayer ( int domnr ) 
     { 
       if ( layer.Size() ) return layer[domnr-1]; 

libsrc/geom2d/python_geom2d.cpp

@@ -392,6 +392,7 @@ DLL_HEADER void ExportGeom2d(py::module &m)
                 }, py::arg("mp") = nullptr,
       py::call_guard<py::gil_scoped_release>(),
       meshingparameter_description.c_str())
+    .def("_SetDomainTensorMeshing", &SplineGeometry2d::SetDomainTensorMeshing)
     ;
   
 }

libsrc/meshing/meshfunc2d.cpp

@@ -20,6 +20,19 @@ namespace netgen
       }
     mesh.Compress();
 
+    bool optimize_swap_separate_faces = false;
+    if(!mp.quad)
+    {
+      bool mixed = false;
+      ParallelFor( Range(mesh.GetNSE()), [&] (auto i) NETGEN_LAMBDA_INLINE
+          {
+            if (mesh[SurfaceElementIndex(i)].GetNP() != 3)
+                mixed = true;
+          });
+      if(mixed)
+        optimize_swap_separate_faces = true;
+    }
+
     const char * optstr = mp.optimize2d.c_str();
     int optsteps = mp.optsteps2d;
 
@@ -31,16 +44,42 @@ namespace netgen
 	    {
 	    case 's': 
 	      {  // topological swap
-		MeshOptimize2d meshopt(mesh);
+                MeshOptimize2d meshopt(mesh);
                 meshopt.SetMetricWeight (mp.elsizeweight);
-		meshopt.EdgeSwapping (0);
+
+                if(optimize_swap_separate_faces)
+                {
+                  for(auto i : Range(1, mesh.GetNFD()+1))
+                  {
+                    meshopt.SetFaceIndex(i);
+                    meshopt.EdgeSwapping (0);
+                  }
+                }
+                else
+                {
+                  meshopt.SetFaceIndex(0);
+                  meshopt.EdgeSwapping (0);
+                }
 		break;
 	      }
 	    case 'S': 
 	      {  // metric swap
 		MeshOptimize2d meshopt(mesh);
                 meshopt.SetMetricWeight (mp.elsizeweight);
-		meshopt.EdgeSwapping (1);
+
+                if(optimize_swap_separate_faces)
+                {
+                  for(auto i : Range(1, mesh.GetNFD()+1))
+                  {
+                    meshopt.SetFaceIndex(i);
+                    meshopt.EdgeSwapping (1);
+                  }
+                }
+                else
+                {
+                  meshopt.SetFaceIndex(0);
+                  meshopt.EdgeSwapping (1);
+                }
 		break;
 	      }
 	    case 'm': 

tests/pytest/test_splinegeo_tensordomainmeshing.py

@@ -0,0 +1,22 @@
+from netgen.geom2d import *
+
+def test_tensordomainmeshing():
+    geo = SplineGeometry()
+    w = 10
+    h = 0.01
+
+    p = [ (0, 0), (w, 0), (w, h), (0, h) ]
+    p = [geo.AppendPoint(*px) for px in p]
+
+    l0 = geo.Append ( ["line", p[0], p[1]], leftdomain=1, rightdomain=0 )
+    l1 = geo.Append ( ["line", p[1], p[2]], leftdomain=1, rightdomain=0)
+    geo.Append ( ["line", p[3], p[2]], leftdomain=0, rightdomain=1, copy=l0 )
+    geo.Append ( ["line", p[0], p[3]], leftdomain=0, rightdomain=1, copy=l1 )
+
+    geo._SetDomainTensorMeshing(1, True)
+
+    mesh = geo.GenerateMesh(maxh=1)
+
+    for el in mesh.Elements2D():
+        print(el.vertices)
+        assert len(el.vertices) == 4