Merge branch 'fix_conform_segments' into 'master'

Fixes to conform to free segments

See merge request ngsolve/netgen!700

libsrc/meshing/improve3.cpp

@@ -2194,7 +2194,7 @@ void MeshOptimize3d :: SwapImproveSurface (
 
 double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
                                         Table<ElementIndex, PointIndex> & elementsonnode,
-                                        DynamicTable<SurfaceElementIndex, PointIndex> & belementsonnode, bool check_only )
+                                        DynamicTable<SurfaceElementIndex, PointIndex> & belementsonnode, bool conform_segments, bool check_only )
 {
   PointIndex pi1, pi2, pi3, pi4, pi5;
   Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET);
@@ -2228,6 +2228,8 @@ double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
   }
 
 
+  if(!conform_segments)
+    {
       bool bface = 0;
       for (int k = 0; k < belementsonnode[pi1].Size(); k++)
       {
@@ -2250,6 +2252,7 @@ double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
       }
 
       if (bface) return 0.0;
+    }
 
 
   FlatArray<ElementIndex> row = elementsonnode[pi1];
@@ -2367,11 +2370,11 @@ double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
   2 -> 3 conversion
 */
 
-void MeshOptimize3d :: SwapImprove2 ()
+void MeshOptimize3d :: SwapImprove2 (bool conform_segments)
 {
   static Timer t("MeshOptimize3d::SwapImprove2"); RegionTimer reg(t);
 
-  if (goal == OPT_CONFORM) return;
+  if (!conform_segments && goal == OPT_CONFORM) return;
 
   mesh.BuildBoundaryEdges(false);
 
@@ -2433,7 +2436,7 @@ void MeshOptimize3d :: SwapImprove2 ()
 
             for (int j = 0; j < 4; j++)
               {
-                double d_badness = SwapImprove2( eli1, j, elementsonnode, belementsonnode, true);
+                double d_badness = SwapImprove2( eli1, j, elementsonnode, belementsonnode, conform_segments, true);
                 if(d_badness<0.0)
                     my_faces_with_improvement.Append( std::make_tuple(d_badness, eli1, j) );
               }
@@ -2449,7 +2452,7 @@ void MeshOptimize3d :: SwapImprove2 ()
     {
       if(mesh[eli].IsDeleted())
           continue;
-      if(SwapImprove2( eli, j, elementsonnode, belementsonnode, false) < 0.0)
+      if(SwapImprove2( eli, j, elementsonnode, belementsonnode, conform_segments, false) < 0.0)
           cnt++;
     }
 

libsrc/meshing/improve3.hpp

@@ -45,8 +45,8 @@ public:
   void SwapImprove (const TBitArray<ElementIndex> * working_elements = NULL);
   void SwapImproveSurface (const TBitArray<ElementIndex> * working_elements = NULL,
 			   const NgArray< idmap_type* > * idmaps = NULL);
-  void SwapImprove2 ();
+  void SwapImprove2 (bool conform_segments = false);
-  double SwapImprove2 (ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, DynamicTable<SurfaceElementIndex, PointIndex> & belementsonnode, bool check_only=false );
+  double SwapImprove2 (ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, DynamicTable<SurfaceElementIndex, PointIndex> & belementsonnode, bool conform_segments, bool check_only=false );
 
   void ImproveMesh() { mesh.ImproveMesh(mp, goal); }
 

libsrc/meshing/meshfunc.cpp

@@ -783,7 +783,7 @@ namespace netgen
         free_segs.Append(segi);
 
     auto get_nonconforming = [&] (const auto & p2el) {
-      Array<size_t> nonconforming;
+      Array<SegmentIndex> nonconforming;
 
       for (auto segi : free_segs) {
         auto seg = mesh[segi];
@@ -805,16 +805,6 @@ namespace netgen
     auto split_segment = [&] (SegmentIndex segi, const auto & p2el) {
       auto seg = mesh[segi];
       auto p_new = Center(mesh[seg[0]], mesh[seg[1]]);
-      auto pi_new = mesh.AddPoint(p_new);
-      auto seg_new0 = seg;
-      auto seg_new1 = seg;
-      seg_new0[1] = pi_new;
-      seg_new1[0] = pi_new;
-
-      mesh[segi][0] = PointIndex::INVALID;
-      mesh.AddSegment(seg_new0);
-      mesh.AddSegment(seg_new1);
-
       double lam[3];
       ElementIndex ei_start = mesh.GetElementOfPoint(p_new, lam, false, domain);
 
@@ -824,6 +814,9 @@ namespace netgen
       }
       ei_start -= 1;
 
+      if(mesh[ei_start].IsDeleted())
+        return;
+
       double max_inside = -1.;
       ElementIndex ei_max_inside = -1;
 
@@ -832,6 +825,9 @@ namespace netgen
       for(auto pi : mesh[ei_start].PNums()) {
         for(auto ei1 : p2el[pi]) {
           double lam[3];
+
+          if(mesh[ei1].IsDeleted())
+            return;
           if(!mesh.PointContainedIn3DElement(p_new, lam, ei1+1))
             continue;
 
@@ -843,18 +839,34 @@ namespace netgen
         }
       }
 
+      if(max_inside < 1e-4) {
+        PrintMessage(3, "Could not find volume element with new point inside");
+        return;
+      }
+
       // split tet into 4 new tests, with new point inside
       auto el = mesh[ei_max_inside];
       if(el.GetNP() != 4) {
-        PrintMessage(1, "Only tet elements are supported to split around free segments");
+        PrintMessage(3, "Only tet elements are supported to split around free segments");
         return;
       }
 
       if(el.IsDeleted()) {
-        PrintMessage(1,"Element to split is already deleted");
+        PrintMessage(3,"Element to split is already deleted");
         return;
       }
 
+      auto pi_new = mesh.AddPoint(p_new);
+      auto seg_new0 = seg;
+      auto seg_new1 = seg;
+      seg_new0[1] = pi_new;
+      seg_new1[0] = pi_new;
+
+      mesh[segi][0] = PointIndex::INVALID;
+      mesh.AddSegment(seg_new0);
+      mesh.AddSegment(seg_new1);
+
+
       int pmap[4][4] = {
         {0,1,2,4},
         {1,3,2,4},
@@ -897,7 +909,7 @@ namespace netgen
 
       for ([[maybe_unused]] auto i : Range(3)) {
         optmesh.ImproveMesh();
-        optmesh.SwapImprove2 ();
+        optmesh.SwapImprove2(true);
         optmesh.ImproveMesh();
         optmesh.SwapImprove();
         optmesh.ImproveMesh();
@@ -910,7 +922,7 @@ namespace netgen
     auto bad_segs = get_nonconforming(p2el);
 
     if(bad_segs.Size() > 0) {
-      auto bad_seg = mesh[free_segs[bad_segs[0]]];
+      auto bad_seg = mesh[bad_segs[0]];
       if(debugparam.write_mesh_on_error)
         mesh.Save("free_segment_not_conformed_dom_"+ToString(domain)+"_seg_"+ToString(bad_seg[0])+"_"+ToString(bad_seg[1])+".vol.gz");
       throw Exception("Segment not resolved in volume mesh in domain " + ToString(domain)+ ", seg: " + ToString(bad_seg));