Merge branch 'free_edges' into 'master'

OCC - support free-floating edges in solids

See merge request ngsolve/netgen!687

libsrc/meshing/basegeom.hpp

@@ -86,6 +86,9 @@ namespace netgen
   protected:
       GeometryVertex *start, *end;
   public:
+    // Neighboring domains in 2d
+    // In 3d unused, EXCEPT for free floating edges in a domain,
+    // then both are pointing to the containing domain
     int domin=-1, domout=-1;
 
     GeometryEdge( GeometryVertex &start_, GeometryVertex &end_ )
@@ -187,7 +190,10 @@ namespace netgen
   };
 
   class DLL_HEADER GeometrySolid : public GeometryShape
-  { };
+  {
+  public:
+    Array<GeometryEdge*> free_edges; // edges with no adjacent face
+  };
 
   class DLL_HEADER NetgenGeometry
   {

libsrc/meshing/delaunay.cpp

@@ -564,7 +564,7 @@ namespace netgen
 
 
 
-  void Delaunay1 (Mesh & mesh, const MeshingParameters & mp, AdFront3 * adfront,
+  void Delaunay1 (Mesh & mesh, int domainnr, const MeshingParameters & mp, AdFront3 * adfront,
 		  NgArray<DelaunayTet> & tempels,
 		  int oldnp, DelaunayTet & startel, Point3d & pmin, Point3d & pmax)
   {
@@ -623,6 +623,13 @@ namespace netgen
     for (PointIndex pi : mesh.LockedPoints())
       usep[pi] = true;
     
+    // mark points of free edge segments (no adjacent face)
+    for (auto & seg : mesh.LineSegments())
+      if(seg.domin == domainnr && seg.domout == domainnr)
+      {
+        usep[seg[0]] = true;
+        usep[seg[1]] = true;
+      }
 
     NgArray<int> freelist;
 
@@ -1554,7 +1561,7 @@ namespace netgen
 
     int np = mesh.GetNP();
 
-    Delaunay1 (mesh, mp, adfront, tempels, oldnp, startel, pmin, pmax);
+    Delaunay1 (mesh, domainnr, mp, adfront, tempels, oldnp, startel, pmin, pmax);
 
     {
       // improve delaunay - mesh by swapping !!!!

libsrc/meshing/improve3.cpp

@@ -14,6 +14,11 @@
 namespace netgen
 {
 
+static constexpr int tetedges[6][2] =
+  { { 0, 1 }, { 0, 2 }, { 0, 3 },
+    { 1, 2 }, { 1, 3 }, { 2, 3 } };
+
+
 static constexpr int IMPROVEMENT_CONFORMING_EDGE = -1e6;
 
 static inline bool NotTooBad(double bad1, double bad2)
@@ -258,6 +263,26 @@ double MeshOptimize3d :: CombineImproveEdge (
   for (auto ei : has_both_points)
       badness_old += mesh[ei].GetBadness();
 
+  if (goal == OPT_CONFORM && p0.Type() <= EDGEPOINT) {
+    // check if the optimization improves conformity with free segments
+    std::set<PointIndex> edges_before, edges_after;
+
+    for (auto ei : has_one_point) {
+        const auto el = mesh[ei];
+        for(auto i : Range(6)) {
+          auto e0 = el[tetedges[i][0]];
+          auto e1 = el[tetedges[i][1]];
+          if(e0 == pi0 || e1 == pi0) edges_before.insert(e0 == pi0 ? e1 : e0);
+          if(e0 == pi1 || e1 == pi1) edges_after.insert(e0 == pi1 ? e1 : e0);
+        }
+    }
+
+    for(auto new_edge : edges_after) {
+      if (edges_before.count(new_edge) == 0 && mesh[new_edge].Type() <= EDGEPOINT && mesh.BoundaryEdge (new_edge, pi0))
+        badness_old += GetLegalPenalty();
+    }
+  }
+
   MeshPoint pnew = p0;
   if (p0.Type() == INNERPOINT)
       pnew = Center (p0, p1);
@@ -1548,10 +1573,6 @@ void MeshOptimize3d :: SwapImproveSurface (
 	  // loop over edges
 
 	  
-	  static const int tetedges[6][2] =
-	    { { 0, 1 }, { 0, 2 }, { 0, 3 },
-	      { 1, 2 }, { 1, 3 }, { 2, 3 } };
-
 	  pi1 = elemi[tetedges[j][0]];
 	  pi2 = elemi[tetedges[j][1]];
 
@@ -2406,6 +2427,8 @@ double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
                   !mesh.LegalTet(elem2))
                   bad1 += GetLegalPenalty();
 
+              if(mesh.BoundaryEdge (pi4, pi5))
+                  bad1 += GetLegalPenalty();
 
               el31.PNum(1) = pi1;
               el31.PNum(2) = pi2;
@@ -2583,10 +2606,6 @@ double MeshOptimize3d :: SplitImprove2Element (
     return false;
 
   // search for very flat tets, with two disjoint edges nearly crossing, like a rectangle with diagonals
-  static constexpr int tetedges[6][2] =
-  { { 0, 1 }, { 0, 2 }, { 0, 3 },
-    { 1, 2 }, { 1, 3 }, { 2, 3 } };
-
   int minedge = -1;
   double mindist = 1e99;
   double minlam0, minlam1;

libsrc/meshing/meshfunc.cpp

@@ -80,6 +80,16 @@ namespace netgen
               m.AddFaceDescriptor( mesh.GetFaceDescriptor(i) );
       }
 
+      // mark interior edge points
+      for(const auto& seg : mesh.LineSegments())
+        {
+          if(seg.domin > 0 && seg.domin == seg.domout)
+            {
+              ipmap[seg.domin-1][seg[0]] = 1;
+              ipmap[seg.domin-1][seg[1]] = 1;
+            }
+        }
+
       // mark used points for each domain, add surface elements (with wrong point numbers) to domain mesh
       for(const auto & sel : mesh.SurfaceElements())
       {
@@ -522,6 +532,7 @@ namespace netgen
         }
      }
      RemoveIllegalElements (mesh, domain);
+     ConformToFreeSegments (mesh, domain);
   }
 
   void MergeMeshes( Mesh & mesh, Array<MeshingData> & md )
@@ -615,6 +626,7 @@ namespace netgen
        {
      ParallelFor( md.Range(), [&](int i)
        {
+         try {
            if (mp.checkoverlappingboundary)
              if (md[i].mesh->CheckOverlappingBoundary())
              {
@@ -627,6 +639,11 @@ namespace netgen
               FillCloseSurface( md[i] );
            CloseOpenQuads( md[i] );
            MeshDomain(md[i]);
+         }
+         catch (const Exception & e) {
+           cerr << "Meshing of domain " << i+1 << " failed with error: " << e.what() << endl;
+           throw e;
+         }
        }, md.Size());
        }
      catch(...)
@@ -734,6 +751,65 @@ namespace netgen
   }
 
 
+  void ConformToFreeSegments (Mesh & mesh, int domain)
+  {
+    auto geo = mesh.GetGeometry();
+    if(!geo) return;
+    auto n_solids = geo->GetNSolids();
+    if(!n_solids) return;
+    if(geo->GetSolid(domain-1).free_edges.Size() == 0)
+      return;
+
+    Segment bad_seg;
+    Array<Segment> free_segs;
+    for (auto seg : mesh.LineSegments())
+      if(seg.domin == domain && seg.domout == domain)
+        free_segs.Append(seg);
+
+    auto num_nonconforming = [&] () {
+      size_t count = 0;
+
+      auto p2el = mesh.CreatePoint2ElementTable();
+
+      for (auto seg : free_segs) {
+
+        auto has_p0 = p2el[seg[0]];
+        bool has_both = false;
+
+        for(auto ei : has_p0) {
+          if(mesh[ei].PNums().Contains(seg[1]))
+            has_both = true;
+        }
+
+        if(!has_both) {
+          bad_seg = seg;
+          count++;
+        }
+      }
+      return count;
+    };
+
+    for (auto i : Range(5)) {
+      auto num_bad_segs = num_nonconforming();
+      PrintMessage(1, "Non-conforming free segments in domain ", domain, ": ", num_bad_segs);
+
+      if(num_bad_segs == 0)
+        return;
+
+      MeshingParameters dummymp;
+      MeshOptimize3d optmesh(mesh, dummymp, OPT_CONFORM);
+
+      for (auto i : Range(3)) {
+        optmesh.SwapImprove2 ();
+        optmesh.SwapImprove();
+        optmesh.CombineImprove();
+      }
+    }
+
+    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));
+  }
 
 
   void RemoveIllegalElements (Mesh & mesh3d, int domain)

libsrc/meshing/meshfunc.hpp

@@ -31,6 +31,7 @@ DLL_HEADER MESHING3_RESULT OptimizeVolume (const MeshingParameters & mp, Mesh& m
 //			       const CSGeometry * geometry = NULL);
 
 DLL_HEADER void RemoveIllegalElements (Mesh & mesh3d, int domain = 0);
+DLL_HEADER void ConformToFreeSegments (Mesh & mesh3d, int domain);
 
 
 enum MESHING_STEP { 

libsrc/occ/occgeom.cpp

@@ -1132,6 +1132,21 @@ namespace netgen
          }
       }
       */
+
+      std::map<int, ArrayMem<int, 10>> free_edges_in_solid;
+      for(auto i1 : Range(1, somap.Extent()+1))
+      {
+          auto s = somap(i1);
+          for (auto edge : MyExplorer(s, TopAbs_EDGE, TopAbs_WIRE))
+            if (!emap.Contains(edge))
+              {
+                free_edges_in_solid[i1].Append(emap.Add (edge));
+                for (auto vertex : MyExplorer(edge, TopAbs_VERTEX))
+                  if (!vmap.Contains(vertex))
+                    vmap.Add (vertex);
+              }
+      }
+
       for (auto edge : MyExplorer(shape, TopAbs_EDGE, TopAbs_WIRE))
         if (!emap.Contains(edge))
           {
@@ -1256,6 +1271,16 @@ namespace netgen
               if(face.Shape().Orientation() == TopAbs_INTERNAL)
                   face.domout = k;
           }
+
+          if(free_edges_in_solid.count(i1))
+            for(auto ei : free_edges_in_solid[i1])
+            {
+              auto & edge = GetEdge(emap(ei));
+              edge.properties.maxh = min(edge.properties.maxh, occ_solid->properties.maxh);
+              edge.domin = k;
+              edge.domout = k;
+              occ_solid->free_edges.Append(&GetEdge(emap(ei)));
+            }
           solids.Append(std::move(occ_solid));
       }