Extra optimization steps for bad elements

libsrc/meshing/improve3.cpp

@@ -136,29 +136,36 @@ static double SplitElementBadness (const Mesh::T_POINTS & points, const MeshingP
 }
 
 
-tuple<double, double> MeshOptimize3d :: UpdateBadness()
+tuple<double, double, int> MeshOptimize3d :: UpdateBadness()
 {
   static Timer tbad("UpdateBadness");
   RegionTimer reg(tbad);
 
   double totalbad = 0.0;
   double maxbad = 0.0;
+  atomic<int> bad_elements = 0;
+
   ParallelForRange(Range(mesh.GetNE()), [&] (auto myrange) {
     double totalbad_local = 0.0;
     double maxbad_local = 0.0;
+    int bad_elements_local = 0;
     for (ElementIndex ei : myrange)
     {
       auto & el = mesh[ei];
       if(mp.only3D_domain_nr && mp.only3D_domain_nr != el.GetIndex()) continue;
       if(!el.BadnessValid())
         el.SetBadness(CalcBad(mesh.Points(), el, 0));
-      totalbad_local += el.GetBadness();
+      double bad = el.GetBadness();
-      maxbad_local = max(maxbad_local, static_cast<double>(el.GetBadness()));
+      totalbad_local += bad;
+      maxbad_local = max(maxbad_local, bad);
+      if(bad > min_badness)
+        bad_elements_local++;
     }
     AtomicAdd(totalbad, totalbad_local);
     AtomicMax(maxbad, maxbad_local);
+    bad_elements += bad_elements_local;
   });
-  return {totalbad, maxbad};
+  return {totalbad, maxbad, bad_elements};
 }
 
 bool MeshOptimize3d :: HasBadElement(FlatArray<ElementIndex> els)
@@ -177,6 +184,13 @@ bool MeshOptimize3d :: HasIllegalElement(FlatArray<ElementIndex> els)
   return false;
 }
 
+bool MeshOptimize3d :: NeedsOptimization(FlatArray<ElementIndex> els)
+{
+  if(goal == OPT_LEGAL) return HasIllegalElement(els);
+  if(goal == OPT_QUALITY) return HasBadElement(els);
+  return true;
+}
+
 
 /*
   Combine two points to one.
@@ -440,7 +454,7 @@ double MeshOptimize3d :: SplitImproveEdge (Table<ElementIndex,PointIndex> & elem
           if(mp.only3D_domain_nr != mesh[ei].GetIndex())
               return 0.0;
 
-  if ((goal == OPT_LEGAL) && !HasIllegalElement(hasbothpoints))
+  if (!NeedsOptimization(hasbothpoints))
     return 0.0;
 
   double bad1 = 0.0;
@@ -733,7 +747,7 @@ double MeshOptimize3d :: SwapImproveEdge (
           return 0.0;
     }
 
-  if ((goal == OPT_LEGAL) && !HasIllegalElement(hasbothpoints))
+  if(!NeedsOptimization(hasbothpoints))
     return 0.0;
 
   int nsuround = hasbothpoints.Size();
@@ -1354,6 +1368,8 @@ void MeshOptimize3d :: SwapImprove (const NgBitArray * working_elements)
   Array<std::tuple<double, int>> candidate_edges(edges.Size());
   std::atomic<int> improvement_counter(0);
 
+  UpdateBadness();
+
   tloop.Start();
 
   auto num_elements_before = mesh.VolumeElements().Range().Next();
@@ -2366,6 +2382,10 @@ double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
           if (elem2.GetType() != TET)
               continue;
 
+          ArrayMem<ElementIndex, 2> elis = {eli1, eli2};
+          if(!NeedsOptimization(elis))
+            continue;
+
           int comnodes=0;
           for (int l = 1; l <= 4; l++)
               if (elem2.PNum(l) == pi1 || elem2.PNum(l) == pi2 ||
@@ -2380,8 +2400,7 @@ double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
 
           if (comnodes == 3)
           {
-              bad1 = CalcBad (mesh.Points(), elem, 0) +
+              bad1 = elem.GetBadness() + elem2.GetBadness();
-                CalcBad (mesh.Points(), elem2, 0);
 
               if (!mesh.LegalTet(elem) ||
                   !mesh.LegalTet(elem2))
@@ -2494,6 +2513,8 @@ void MeshOptimize3d :: SwapImprove2 ()
   Array<std::tuple<double, ElementIndex, int>> faces_with_improvement;
   Array<Array<std::tuple<double, ElementIndex, int>>> faces_with_improvement_threadlocal(num_threads);
 
+  UpdateBadness();
+
   ParallelForRange( Range(ne), [&]( auto myrange )
       {
         int tid = ngcore::TaskManager::GetThreadId();

libsrc/meshing/improve3.hpp

@@ -17,9 +17,9 @@ class MeshOptimize3d
   OPTIMIZEGOAL goal = OPT_QUALITY;
   double min_badness = 0;
 
-  tuple<double, double> UpdateBadness();
   bool HasBadElement(FlatArray<ElementIndex> els);
   bool HasIllegalElement(FlatArray<ElementIndex> els);
+  bool NeedsOptimization(FlatArray<ElementIndex> els);
 
 public:
 
@@ -29,6 +29,8 @@ public:
   void SetGoal(OPTIMIZEGOAL agoal) { goal = agoal; }
   void SetMinBadness(double badness) { min_badness = badness; }
 
+  tuple<double, double, int> UpdateBadness();
+
   double CombineImproveEdge (
             Table<ElementIndex, PointIndex> & elements_of_point,
             PointIndex pi0, PointIndex pi1,

libsrc/meshing/meshfunc.cpp

@@ -654,13 +654,28 @@ namespace netgen
     */
 
     mesh3d.CalcSurfacesOfNode();
+
+    MeshOptimize3d optmesh(mesh3d, mp);
+
+    // optimize only bad elements first
+    optmesh.SetMinBadness(1000.);
+    for(auto i : Range(mp.optsteps3d))
+      {
+        auto [total_badness, max_badness, bad_els] = optmesh.UpdateBadness();
+        if(bad_els==0) break;
+        optmesh.SplitImprove();
+        optmesh.SwapImprove();
+        optmesh.SwapImprove2();
+      }
+
+    // Now optimize all elements
+    optmesh.SetMinBadness(0);
+
     for (auto i : Range(mp.optsteps3d))
       {
 	if (multithread.terminate)
 	  break;
 
-	MeshOptimize3d optmesh(mesh3d, mp);
-
 	// teterrpow = mp.opterrpow;
 	// for (size_t j = 1; j <= strlen(mp.optimize3d); j++)
         for (auto j : Range(mp.optimize3d.size()))

tests/pytest/compare_results.py

@@ -88,29 +88,34 @@ for bad1,bad2, f1, f2 in zip(data['#trigs'], data2['#trigs'], data['file'], data
     if bad2>0 and bad2<0.8*bad1:
         print(f"{GREEN}ntrigs {f1} got better: {bad1} -> {bad2}".ljust(w) + diff + RESET)
 
-n = len(data)+1
+n = len(data) + 1
-fig,ax = plt.subplots(figsize=(10,7))
+fig, ax = plt.subplots(figsize=(15, 7))  # Adjust figsize as needed
-for i,d in enumerate(['min trig angle','min tet angle','max trig angle','max tet angle']):
+plt.xticks([])
-    ax = plt.subplot(2,5,i+1)
+plt.yticks([])
+ax.yaxis.grid(False)
+ax.xaxis.grid(False)
+for i, d in enumerate(['min trig angle', 'min tet angle', 'max trig angle', 'max tet angle']):
+    plt.subplot(2, 4, i + 1)  # Remove ax = 
     plt.title(d)
-    ax.set_xticks([1,2])
+    # plt.xticks([1, 2])
-    if len(data[d])==0 or len(data2[d])==0:
+    if len(data[d]) == 0 or len(data2[d]) == 0:
         continue
 
-    plt.violinplot([data[d],data2[d]], showmedians=True)
+    plt.violinplot([data[d], data2[d]], showmedians=True)
     med = statistics.median(data[d])
-    plt.hlines(med, 1,2, linestyle='dotted')
+    plt.hlines(med, 1, 2, linestyle='dotted')
-    if d=='badness':
+    if d == 'badness':
-        ax.set_yscale('log')
+        plt.yscale('log')
-    ax.set_xticklabels([ref, ref2])
+    plt.xticks([1, 2], [ref, ref2])
 
+for i, d in enumerate(['badness', '#edges', '#trigs', '#tets']):
+    plt.xticks([])
+    plt.subplot(2, 4, 5 + i)
+    plt.title('difference ' + d + ' (in %)')
+    plt.boxplot([100.0 * (y - x) / x for x, y in zip(data[d], data2[d])])
+    plt.hlines(0.0, 0.5, 1.5, linestyle='dotted')
 
-for i,d in enumerate(['badness','#edges','#trigs','#tets']):
+plt.tight_layout()  # Adjust layout
-    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([100.0*(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)

tests/pytest/test_nonnative_master

@@ -0,0 +1,93 @@
+tstart (must be vertically): (0, 1.82574)
+tend (must be vertically): (0, -1.82574)
+sin (must not be 0) = 0.447214sin (must not be 0) = 0.447214tstart (must be vertically): (0, 1.82574)
+tend (must be vertically): (0, -1.82574)
+sin (must not be 0) = 0.447214sin (must not be 0) = 0.447214tstart (must be vertically): (0, 1.82574)
+tend (must be vertically): (0, -1.82574)
+sin (must not be 0) = 0.447214sin (must not be 0) = 0.447214tstart (must be vertically): (0, 1.82574)
+tend (must be vertically): (0, -1.82574)
+sin (must not be 0) = 0.447214sin (must not be 0) = 0.447214tstart (must be vertically): (0, 1.82574)
+tend (must be vertically): (0, -1.82574)
+sin (must not be 0) = 0.447214sin (must not be 0) = 0.447214tstart (must be vertically): (0, 1.82574)
+tend (must be vertically): (0, -1.82574)
+generate boundarycondition.geo
+needed 0.12342000007629395 seconds
+generate boxcyl.geo
+needed 0.801245927810669 seconds
+generate circle_on_cube.geo
+needed 0.48620080947875977 seconds
+generate cone.geo
+needed 1.186929702758789 seconds
+generate cube.geo
+needed 0.09043073654174805 seconds
+generate cubeandring.geo
+needed 2.1420021057128906 seconds
+generate cubeandspheres.geo
+needed 0.26427721977233887 seconds
+generate cubemcyl.geo
+needed 18.373886108398438 seconds
+generate cubemsphere.geo
+needed 3.6954052448272705 seconds
+generate cylinder.geo
+needed 0.44164204597473145 seconds
+generate cylsphere.geo
+needed 0.8774328231811523 seconds
+generate ellipsoid.geo
+needed 1.4510962963104248 seconds
+generate ellipticcone.geo
+needed 3.0906074047088623 seconds
+generate ellipticcyl.geo
+needed 2.0780415534973145 seconds
+generate extrusion.geo
+needed 0.40680599212646484 seconds
+generate fichera.geo
+needed 0.1265270709991455 seconds
+generate hinge.stl
+needed 5.519087553024292 seconds
+generate lense.in2d
+needed 0.05641365051269531 seconds
+generate lshape3d.geo
+needed 0.09937620162963867 seconds
+generate manyholes.geo
+needed 9.43863034248352 seconds
+generate manyholes2.geo
+needed 7.40019965171814 seconds
+generate matrix.geo
+needed 3.734792709350586 seconds
+generate ortho.geo
+needed 0.09516167640686035 seconds
+generate part1.stl
+needed 2.6940107345581055 seconds
+generate period.geo
+needed 2.709449291229248 seconds
+generate revolution.geo
+needed 7.7064368724823 seconds
+generate sculpture.geo
+needed 0.6283819675445557 seconds
+generate shaft.geo
+needed 2.921243190765381 seconds
+generate shell.geo
+generate sphere.geo
+needed 0.18424725532531738 seconds
+generate sphereincube.geo
+needed 0.6060984134674072 seconds
+generate square.in2d
+needed 0.021883487701416016 seconds
+generate squarecircle.in2d
+needed 0.04081606864929199 seconds
+generate squarehole.in2d
+needed 0.03681302070617676 seconds
+generate torus.geo
+needed 6.590093612670898 seconds
+generate trafo.geo
+needed 3.2712368965148926 seconds
+generate twobricks.geo
+needed 0.13849091529846191 seconds
+generate twocubes.geo
+needed 0.13361692428588867 seconds
+generate twocyl.geo
+needed 0.8036918640136719 seconds
+generate plane.stl
+needed 15.712460041046143 seconds
+done
+sin (must not be 0) = 0.447214sin (must not be 0) = 0.447214