Merge branch 'master' into refactor_geo_meshing

libsrc/meshing/improve3.cpp

@@ -1698,7 +1698,7 @@ void MeshOptimize3d :: SwapImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal,
 }
 
 
-bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
+double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
         const NgBitArray * working_elements,
         Table<ElementIndex, PointIndex> & elementsonnode,
         INDEX_3_HASHTABLE<int> & faces,
@@ -1714,10 +1714,10 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
   Element el1b(TET), el2b(TET), el3b(TET), el4b(TET);
   ArrayMem<ElementIndex, 20> hasbothpoints;
 
-  bool do_swap = false;
+  double d_badness = 0.0;
   if (pi2 < pi1) Swap (pi1, pi2);
 
-  if (mesh.BoundaryEdge (pi1, pi2)) return false;
+  if (mesh.BoundaryEdge (pi1, pi2)) return 0.0;
 
 
   hasbothpoints.SetSize (0);
@@ -1726,7 +1726,7 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
       bool has1 = 0, has2 = 0;
       const Element & elem = mesh[elnr];
 
-      if (elem.IsDeleted()) return false;
+      if (elem.IsDeleted()) return 0.0;
 
       for (int l = 0; l < elem.GetNP(); l++)
         {
@@ -1749,27 +1749,27 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
   for (ElementIndex ei : hasbothpoints)
     {
       if (mesh[ei].GetType () != TET)
-          return false;
+          return 0.0;
 
       if (mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex())
-          return false;
+          return 0.0;
 
 
       if ((mesh.ElementType(ei)) == FIXEDELEMENT)
-          return false;
+          return 0.0;
 
       if(working_elements &&
               ei < working_elements->Size() &&
               !working_elements->Test(ei))
-          return false;
+          return 0.0;
 
       if (mesh[ei].IsDeleted())
-          return false;
+          return 0.0;
 
       if ((goal == OPT_LEGAL) &&
               mesh.LegalTet (mesh[ei]) &&
               CalcBad (mesh.Points(), mesh[ei], 0) < 1e3)
-          return false;
+          return 0.0;
     }
 
   int nsuround = hasbothpoints.Size();
@@ -1883,8 +1883,9 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
         {
           //		  (*mycout) << "3->2 " << flush;
           //		  (*testout) << "3->2 conversion" << endl;
-          do_swap = true;
+          d_badness = bad2-bad1;
-          if(check_only) return do_swap;
+          if(check_only)
+              return d_badness;
 
 
           /*
@@ -2079,12 +2080,9 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
           swap3 = !swap2 && (bad3 < bad1);
         }
 
-
+      d_badness = swap2 ? bad2-bad1 : bad3-bad1;
-      if (swap2 || swap3)
+      if(check_only)
-        {
+          return d_badness;
-          do_swap = true;
-          if(check_only) return do_swap;
-        }
 
       if (swap2)
         {
@@ -2279,8 +2277,9 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
       if (bestl != -1)
         {
           // (*mycout) << nsuround << "->" << 2 * (nsuround-2) << " " << flush;
-          do_swap = true;
+          d_badness = badopt-bad1;
-          if(check_only) return do_swap;
+          if(check_only)
+              return d_badness;
 
           for (int k = bestl+1; k <= nsuround + bestl - 2; k++)
             {
@@ -2323,7 +2322,7 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
             }
         }
     }
-  return do_swap;
+  return d_badness;
 }
 
 void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
@@ -2373,7 +2372,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
   Array<std::tuple<PointIndex,PointIndex>> edges;
   BuildEdgeList(mesh, elementsonnode, edges);
 
-  Array<int> candidate_edges(edges.Size());
+  Array<std::tuple<double, int>> candidate_edges(edges.Size());
   std::atomic<int> improvement_counter(0);
 
   tloop.Start();
@@ -2386,18 +2385,22 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
         break;
 
       auto [pi0, pi1] = edges[i];
-      if(SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, true))
+      double d_badness = SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, true);
-        candidate_edges[improvement_counter++] = i;
+      if(d_badness<0.0)
+      {
+        int index = improvement_counter++;
+        candidate_edges[index] = make_tuple(d_badness, i);
+      }
     }
   });
 
   auto edges_with_improvement = candidate_edges.Part(0, improvement_counter.load());
   QuickSort(edges_with_improvement);
 
-  for(auto ei : edges_with_improvement)
+  for(auto [d_badness, ei] : edges_with_improvement)
   {
       auto [pi0,pi1] = edges[ei];
-        if(SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, false))
+      if(SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, false) < 0.0)
           cnt++;
   }
 
@@ -3277,7 +3280,7 @@ void MeshOptimize3d :: SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal,
   2 -> 3 conversion
 */
 
-bool MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementIndex eli1, int face,
+double MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementIndex eli1, int face,
   Table<ElementIndex, PointIndex> & elementsonnode,
   TABLE<SurfaceElementIndex, PointIndex::BASE> & belementsonnode, bool check_only )
 {
@@ -3285,9 +3288,10 @@ bool MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementInd
   Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET);
   int j = face;
   double bad1, bad2;
+  double d_badness = 0.0;
 
   Element & elem = mesh[eli1];
-  if (elem.IsDeleted()) return false;
+  if (elem.IsDeleted()) return 0.0;
 
   int mattyp = elem.GetIndex();
 
@@ -3333,7 +3337,7 @@ bool MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementInd
       }
   }
 
-  if (bface) return false;
+  if (bface) return 0.0;
 
 
   FlatArray<ElementIndex> row = elementsonnode[pi1];
@@ -3403,14 +3407,16 @@ bool MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementInd
                   bad2 += 1e4;
 
 
-              bool do_swap = (bad2 < bad1);
+              d_badness = bad2 - bad1;
 
               if ( ((bad2 < 1e6) || (bad2 < 10 * bad1)) &&
                   mesh.BoundaryEdge (pi4, pi5))
-                  do_swap = 1;
+                  d_badness = -1e4;
 
+              if(check_only)
+                  return d_badness;
 
-              if (!check_only && do_swap)
+              if (d_badness<0.0)
               {
                   el31.flags.illegal_valid = 0;
                   el32.flags.illegal_valid = 0;
@@ -3422,12 +3428,11 @@ bool MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementInd
                   mesh.AddVolumeElement (el32);
                   mesh.AddVolumeElement (el33);
               }
-              return do_swap;
+              return d_badness;
           }
       }
   }
-
+  return d_badness;
-  return false;
 }
 
 /*
@@ -3570,8 +3575,8 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
 
   int num_threads = ngcore::TaskManager::GetNumThreads();
 
-  Array<std::tuple<ElementIndex, int>> faces_with_improvement;
+  Array<std::tuple<double, ElementIndex, int>> faces_with_improvement;
-  Array<Array<std::tuple<ElementIndex, int>>> faces_with_improvement_threadlocal(num_threads);
+  Array<Array<std::tuple<double, ElementIndex, int>>> faces_with_improvement_threadlocal(num_threads);
 
   ParallelForRange( Range(ne), [&]( auto myrange )
       {
@@ -3598,8 +3603,9 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
 
             for (int j = 0; j < 4; j++)
               {
-                if(SwapImprove2( mesh, goal, eli1, j, elementsonnode, belementsonnode, true))
+                double d_badness = SwapImprove2( mesh, goal, eli1, j, elementsonnode, belementsonnode, true);
-                  my_faces_with_improvement.Append( std::make_tuple(eli1, j) );
+                if(d_badness<0.0)
+                    my_faces_with_improvement.Append( std::make_tuple(d_badness, eli1, j) );
               }
           }
       });
@@ -3609,8 +3615,9 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
 
   QuickSort(faces_with_improvement);
 
-  for (auto [eli,j] : faces_with_improvement)
+  for (auto [dummy, eli,j] : faces_with_improvement)
-    cnt += SwapImprove2( mesh, goal, eli, j, elementsonnode, belementsonnode, false);
+      if(SwapImprove2( mesh, goal, eli, j, elementsonnode, belementsonnode, false) < 0.0)
+          cnt++;
 
   PrintMessage (5, cnt, " swaps performed");
 

libsrc/meshing/improve3.hpp

@@ -25,7 +25,7 @@ public:
 
   void SplitImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
 
-  bool SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, const NgBitArray * working_elements, Table<ElementIndex,PointIndex> & elementsonnode, INDEX_3_HASHTABLE<int> & faces, PointIndex pi1, PointIndex pi2, bool check_only=false);
+  double SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, const NgBitArray * working_elements, Table<ElementIndex,PointIndex> & elementsonnode, INDEX_3_HASHTABLE<int> & faces, PointIndex pi1, PointIndex pi2, bool check_only=false);
   void SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,
 		    const NgBitArray * working_elements = NULL);
   void SwapImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,
@@ -35,7 +35,7 @@ public:
 			   const NgArray< NgArray<int,PointIndex::BASE>* > * idmaps = NULL);
   void SwapImprove2Sequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
   void SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
-  bool SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, TABLE<SurfaceElementIndex, PointIndex::BASE> & belementsonnode, bool check_only=false );
+  double SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, TABLE<SurfaceElementIndex, PointIndex::BASE> & belementsonnode, bool check_only=false );
 
   double 
   CalcBad (const Mesh::T_POINTS & points, const Element & elem, double h)

tests/pytest/compare_results.py

@@ -0,0 +1,91 @@
+import json
+import sys
+import subprocess
+import statistics
+
+def readData(a, files):
+    amin=[]
+    amax=[]
+    amin1=[]
+    amax1=[]
+    bad=[]
+    ne1d=[]
+    ne2d=[]
+    ne3d=[]
+    for f in files:
+        for t in a[f]:
+            if t['ne1d']>0:
+                ne1d.append(t['ne1d'])
+            if t['ne2d']>0:
+                ne2d.append(t['ne2d'])
+            if t['ne3d']>0:
+                ne3d.append(t['ne3d'])
+            if t['total_badness']>0.0:
+                bad.append(t['total_badness'])
+                if 'angles_tet' in t:
+                    amin.append(t['angles_tet'][0])
+                    amax.append(t['angles_tet'][1])
+                if 'angles_trig' in t:
+                    amin1.append(t['angles_trig'][0])
+                    amax1.append(t['angles_trig'][1])
+    return {
+            "min tet angle":amin,
+            "max tet angle" : amax,
+            "min trig angle":amin1,
+            "max trig angle" : amax1,
+            "badness" : bad,
+            "#edges" : ne1d,
+            "#trigs" : ne2d,
+            "#tets" : ne3d,
+            }
+
+import matplotlib.pyplot as plt
+
+ref = 'master'
+if len(sys.argv)>1:
+    ref = sys.argv[1]
+
+res = subprocess.run(['git','show','{}:./results.json'.format(ref)], capture_output=True)
+s = json.loads(res.stdout.decode())
+
+if len(sys.argv) > 2:
+    ref2 = sys.argv[2]
+    res = subprocess.run(['git','show','{}:./results.json'.format(ref2)], capture_output=True)
+    s2 = res.stdout.decode()
+else:
+    ref2 = 'current'
+    s2 = open('results.json','r').read()
+s2 = json.loads(s2)
+
+filenames = [f for f in s if f in s2]
+data = readData(s, filenames)
+data2 = readData(s2, filenames)
+
+n = len(data)+1
+fig,ax = plt.subplots(figsize=(10,7))
+for i,d in enumerate(['min trig angle','min tet angle','max trig angle','max tet angle']):
+    ax = plt.subplot(2,5,i+1)
+    plt.title(d)
+    ax.set_xticks([1,2])
+    if len(data[d])==0 or len(data2[d])==0:
+        continue
+
+    plt.violinplot([data[d],data2[d]], showmedians=True)
+    med = statistics.median(data[d])
+    plt.hlines(med, 1,2, linestyle='dotted')
+    if d=='badness':
+        ax.set_yscale('log')
+    ax.set_xticklabels([ref, ref2])
+
+
+for i,d in enumerate(['badness','#edges','#trigs','#tets']):
+    ax = plt.subplot(2,5,6+i)
+    plt.title('difference '+d+' (in %)')
+#     plt.violinplot([(y-x)/x for x,y in zip(data[d],data2[d])], showmedians=True)
+    plt.boxplot([(y-x)/x for x,y in zip(data[d],data2[d])])
+    plt.hlines(0.0, 0.5,1.5, linestyle='dotted')
+
+
+# plt.savefig('comparison.png', dpi=100)
+plt.show()
+ 

tests/pytest/test_tutorials.py

@@ -14,13 +14,23 @@ except ImportError:
 
 SetMessageImportance(0)
 
+def round(x, digits=11):
+    try:
+        return float(("{:."+str(digits)+"g}").format(x))
+    except: #list
+        return [float(("{:."+str(digits)+"g}").format(y)) for y in x]
+
+
 def getData(mesh, mp):
     out = {}
     out['ne1d'] = len(mesh.Elements1D())
     out['ne2d'] = len(mesh.Elements2D())
     out['ne3d'] = len(mesh.Elements3D())
     # round badness to avoid fluctuations in last digits
-    out["total_badness"] = float("{:.11g}".format(mesh.CalcTotalBadness(mp)))
+    out["total_badness"] = round(mesh.CalcTotalBadness(mp))
+    angles = mesh.CalcMinMaxAngle()
+    out["angles_trig"] = round(angles["trig"], 5)
+    out["angles_tet"] = round(angles["tet"], 5)
     out["quality_histogram"] = str(list(mesh.GetQualityHistogram()))
     return out
 
@@ -31,6 +41,8 @@ def checkData(mesh, mp, ref):
     assert ref['ne3d'] == data['ne3d']
     assert json.loads(ref['quality_histogram']) == pytest.approx(json.loads(data['quality_histogram']), abs=1, rel=0.4)
     assert ref['total_badness'] == pytest.approx(data['total_badness'], rel=1e-5)
+    assert ref['angles_trig'] == pytest.approx(data['angles_trig'], rel=1e-5)
+    assert ref['angles_tet'] == pytest.approx(data['angles_tet'], rel=1e-5)
 
 # get tutorials
 def getFiles(fileEnding):