Parallel SplitImprove2, update test results

Due to prallelization, the order of splits is changed (sort by
    improvement of badness, like in other optimization passes)

libsrc/meshing/improve3.cpp

@@ -95,6 +95,51 @@ static ArrayMem<Element, 3> SplitElement (Element old, PointIndex pi0, PointInde
   return new_elements;
 };
 
+static double SplitElementBadness (const Mesh::T_POINTS & points, const MeshingParameters & mp, Element old, PointIndex pi0, PointIndex pi1, MeshPoint & pnew)
+{
+  double badness = 0;
+  auto np = old.GetNP();
+  PointIndex dummy{-1};
+  if(np == 4)
+  {
+    // Split tet into two tets
+    badness += CalcBadReplacePoints ( points, mp, old, 0, pi0, dummy, pnew );
+    badness += CalcBadReplacePoints ( points, mp, old, 0, pi1, dummy, pnew );
+  }
+  else if (np == 5)
+  {
+    // split pyramid into pyramid and two tets
+    auto pibase = (pi0==old[4]) ? pi1 : pi0;
+    auto pitop = (pi0==old[4]) ? pi0 : pi1;
+
+    badness += CalcBadReplacePoints ( points, mp, old, 0, pitop, dummy, pnew );
+
+    Element tet = old;
+    tet.SetType(TET);
+
+    size_t pibase_index=0;
+    for(auto i : Range(4))
+      if(old[i]==pibase)
+        pibase_index = i;
+
+    MeshPoint p[4];
+    p[0] = points[old[(pibase_index+1)%4]];
+    p[1] = points[old[(pibase_index+2)%4]];
+    p[2] = pnew;
+    p[3] = points[pitop];
+    badness += CalcTetBadness (p[0], p[1], p[2], p[3], 0, mp);
+
+    p[0] = points[old[(pibase_index+2)%4]];
+    p[1] = points[old[(pibase_index+3)%4]];
+    p[2] = pnew;
+    p[3] = points[pitop];
+    badness += CalcTetBadness (p[0], p[1], p[2], p[3], 0, mp);
+  }
+
+  return badness;
+};
+
+
 /*
   Combine two points to one.
   Set new point into the center, if both are
@@ -3970,13 +4015,148 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
   (*testout) << "swapimprove2 done" << "\n";
 }
 
+double MeshOptimize3d :: SplitImprove2Element (Mesh & mesh,
+                            ElementIndex ei,
+                            const Table<ElementIndex, PointIndex> & elements_of_point,
+                            const Array<double> & el_badness,
+                            bool check_only)
+{
+  auto & el = mesh[ei];
+  if(el.GetType() != TET)
+    return false;
+
+  // Optimize only bad elements
+  if(el_badness[ei] < 100)
+    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;
+
+  for (int i : Range(3))
+  {
+    auto pi0 = el[tetedges[i][0]];
+    auto pi1 = el[tetedges[i][1]];
+    auto pi2 = el[tetedges[5-i][0]];
+    auto pi3 = el[tetedges[5-i][1]];
+
+    double lam0, lam1;
+    double dist = MinDistLL2(mesh[pi0], mesh[pi1], mesh[pi2], mesh[pi3], lam0, lam1 );
+    if(dist<mindist)
+    {
+      mindist = dist;
+      minedge = i;
+      minlam0 = lam0;
+      minlam1 = lam1;
+    }
+  }
+
+  if(minedge==-1)
+    return false;
+
+  auto pi0 = el[tetedges[minedge][0]];
+  auto pi1 = el[tetedges[minedge][1]];
+  auto pi2 = el[tetedges[5-minedge][0]];
+  auto pi3 = el[tetedges[5-minedge][1]];
+
+  // we cannot split edges on the boundary
+  if(mesh.BoundaryEdge (pi0,pi1) || mesh.BoundaryEdge(pi2, pi3))
+    return false;
+
+  ArrayMem<ElementIndex, 50> has_both_points0;
+  ArrayMem<ElementIndex, 50> has_both_points1;
+
+  Point3d p[4] = { mesh[el[0]], mesh[el[1]], mesh[el[2]], mesh[el[3]] };
+  auto center = Center(p[0]+minlam0*(p[1]-p[0]), p[2]+minlam1*(p[3]-p[2]));
+  MeshPoint pnew;
+
+  pnew(0) = center.X();
+  pnew(1) = center.Y();
+  pnew(2) = center.Z();
+
+  // find all tets with edge (pi0,pi1) or (pi2,pi3)
+  for (auto ei0 : elements_of_point[pi0] )
+  {
+    Element & elem = mesh[ei0];
+    if (elem.IsDeleted()) return false;
+    if (ei0 == ei) continue;
+
+    if (elem[0] == pi1 || elem[1] == pi1 || elem[2] == pi1 || elem[3] == pi1 || (elem.GetNP()==5 && elem[4]==pi1) )
+      if(!has_both_points0.Contains(ei0))
+        has_both_points0.Append (ei0);
+  }
+
+  for (auto ei1 : elements_of_point[pi2] )
+  {
+    Element & elem = mesh[ei1];
+    if (elem.IsDeleted()) return false;
+    if (ei1 == ei) continue;
+
+    if (elem[0] == pi3 || elem[1] == pi3 || elem[2] == pi3 || elem[3] == pi3 || (elem.GetNP()==5 && elem[4]==pi3))
+      if(!has_both_points1.Contains(ei1))
+        has_both_points1.Append (ei1);
+  }
+
+  double badness_before = el_badness[ei];
+  double badness_after = 0.0;
+
+  for (auto ei0 : has_both_points0)
+  {
+    badness_before += el_badness[ei0];
+    badness_after += SplitElementBadness (mesh.Points(), mp, mesh[ei0], pi0, pi1, pnew);
+  }
+  for (auto ei1 : has_both_points1)
+  {
+    badness_before += el_badness[ei1];
+    badness_after += SplitElementBadness (mesh.Points(), mp, mesh[ei1], pi2, pi3, pnew);
+  }
+
+  if(check_only)
+    return badness_after-badness_before;
+
+  if(badness_after<badness_before)
+  {
+    PointIndex pinew = mesh.AddPoint (center);
+    el.flags.illegal_valid = 0;
+    el.Delete();
+
+    for (auto ei1 : has_both_points0)
+    {
+      auto new_els = SplitElement(mesh[ei1], pi0, pi1, pinew);
+      for(const auto & el : new_els)
+        mesh.AddVolumeElement(el);
+      mesh[ei1].Delete();
+    }
+    for (auto ei1 : has_both_points1)
+    {
+      auto new_els = SplitElement(mesh[ei1], pi2, pi3, pinew);
+      for(const auto & el : new_els)
+        mesh.AddVolumeElement(el);
+      mesh[ei1].Delete();
+    }
+  }
+  return badness_after-badness_before;
+}
+
 // Split two opposite edges of very flat tet and let all 4 new segments have one common vertex
 // Imagine a square with 2 diagonals -> new point where diagonals cross, remove the flat tet
-void MeshOptimize3d :: SplitImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
+void MeshOptimize3d :: SplitImprove2 (Mesh & mesh)
 {
-  double bad1 = mesh.CalcTotalBad (mp);
+  static Timer t("MeshOptimize3d::SplitImprove2"); RegionTimer reg(t);
+  static Timer tsearch("Search");
+  static Timer topt("Optimize");
+
   int ne = mesh.GetNE();
   auto elements_of_point = mesh.CreatePoint2ElementTable();
+  int ntasks = 4*ngcore::TaskManager::GetNumThreads();
+
+  const char * savetask = multithread.task;
+  multithread.task = "Optimize Volume: Split Improve 2";
 
   Array<double> el_badness (ne);
 
@@ -3985,141 +4165,51 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
     for (ElementIndex ei : myrange)
     {
       if(mp.only3D_domain_nr && mp.only3D_domain_nr != mesh[ei].GetIndex())
-	continue;
+        continue;
       el_badness[ei] = CalcBad (mesh.Points(), mesh[ei], 0);
     }
   });
 
   mesh.BoundaryEdge (1,2); // ensure the boundary-elements table is built
 
-  // TODO: Parallelization
+  Array<std::tuple<double, ElementIndex>> split_candidates(ne);
-  for (ElementIndex ei : Range(ne) )
+  std::atomic<int> improvement_counter(0);
+
+  tsearch.Start();
+  ParallelForRange(Range(ne), [&] (auto myrange)
   {
-    auto & el = mesh[ei];
+    for(ElementIndex ei : myrange)
-    if(el.GetType() != TET)
-      continue;
-
-    // Optimize only bad elements
-    if(el_badness[ei] < 100)
-      continue;
-
-    // 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;
-
-    for (int i : Range(3))
     {
-      auto pi0 = el[tetedges[i][0]];
+      if(mp.only3D_domain_nr && mp.only3D_domain_nr != mesh[ei].GetIndex())
-      auto pi1 = el[tetedges[i][1]];
+        continue;
-      auto pi2 = el[tetedges[5-i][0]];
+      double d_badness = SplitImprove2Element(mesh, ei, elements_of_point, el_badness, true);
-      auto pi3 = el[tetedges[5-i][1]];
+      if(d_badness<0.0)
-
-      double lam0, lam1;
-      double dist = MinDistLL2(mesh[pi0], mesh[pi1], mesh[pi2], mesh[pi3], lam0, lam1 );
-      if(dist<mindist)
       {
-        mindist = dist;
+        int index = improvement_counter++;
-        minedge = i;
+        split_candidates[index] = make_tuple(d_badness, ei);
-        minlam0 = lam0;
-        minlam1 = lam1;
       }
     }
+  }, ntasks);
+  tsearch.Stop();
 
-    if(minedge==-1)
+  auto elements_with_improvement = split_candidates.Part(0, improvement_counter.load());
-      continue;
+  QuickSort(elements_with_improvement);
 
-    auto pi0 = el[tetedges[minedge][0]];
+  size_t cnt = 0;
-    auto pi1 = el[tetedges[minedge][1]];
+  topt.Start();
-    auto pi2 = el[tetedges[5-minedge][0]];
+  for(auto [d_badness, ei] : elements_with_improvement)
-    auto pi3 = el[tetedges[5-minedge][1]];
+  {
-
+    if( SplitImprove2Element(mesh, ei, elements_of_point, el_badness, false) < 0.0)
-    // we cannot split edges on the boundary
+      cnt++;
-    if(mesh.BoundaryEdge (pi0,pi1) || mesh.BoundaryEdge(pi2, pi3))
-      continue;
-
-    ArrayMem<ElementIndex, 50> has_both_points0;
-    ArrayMem<ElementIndex, 50> has_both_points1;
-
-    Point3d p[4] = { mesh[el[0]], mesh[el[1]], mesh[el[2]], mesh[el[3]] };
-    auto center = Center(p[0]+minlam0*(p[1]-p[0]), p[2]+minlam1*(p[3]-p[2]));
-    MeshPoint pnew;
-
-    pnew(0) = center.X();
-    pnew(1) = center.Y();
-    pnew(2) = center.Z();
-
-    bool valid = true;
-    // find all tets with edge (pi0,pi1) or (pi2,pi3)
-    for (auto ei0 : elements_of_point[pi0] )
-    {
-      Element & elem = mesh[ei0];
-      if (elem.IsDeleted()) valid = false;
-      if (ei0 == ei) continue;
-
-      if (elem[0] == pi1 || elem[1] == pi1 || elem[2] == pi1 || elem[3] == pi1 || (elem.GetNP()==5 && elem[4]==pi1) )
-        if(!has_both_points0.Contains(ei0))
-          has_both_points0.Append (ei0);
-    }
-
-    for (auto ei1 : elements_of_point[pi2] )
-    {
-      Element & elem = mesh[ei1];
-      if (elem.IsDeleted()) valid = false;
-      if (ei1 == ei) continue;
-
-      if (elem[0] == pi3 || elem[1] == pi3 || elem[2] == pi3 || elem[3] == pi3 || (elem.GetNP()==5 && elem[4]==pi3))
-        if(!has_both_points1.Contains(ei1))
-          has_both_points1.Append (ei1);
-    }
-    if(!valid) continue;
-    double badness_before = el_badness[ei];
-
-    double badness_after = 0.0;
-    PointIndex dummy{-1};
-
-    PointIndex pinew = mesh.AddPoint (center);
-    for (auto ei0 : has_both_points0)
-    {
-      badness_before += el_badness[ei0];
-      auto new_els = SplitElement(mesh[ei0], pi0, pi1, pinew);
-      for(const auto & el : new_els)
-        badness_after += CalcBad (mesh.Points(), el, 0);
-    }
-    for (auto ei1 : has_both_points1)
-    {
-      badness_before += el_badness[ei1];
-      auto new_els = SplitElement(mesh[ei1], pi2, pi3, pinew);
-      for(const auto & el : new_els)
-        badness_after += CalcBad (mesh.Points(), el, 0);
-    }
-    if(badness_after<badness_before)
-    {
-      el.flags.illegal_valid = 0;
-      el.Delete();
-
-      for (auto ei1 : has_both_points0)
-      {
-        auto new_els = SplitElement(mesh[ei1], pi0, pi1, pinew);
-        for(const auto & el : new_els)
-          mesh.AddVolumeElement(el);
-        mesh[ei1].Delete();
-      }
-      for (auto ei1 : has_both_points1)
-      {
-        auto new_els = SplitElement(mesh[ei1], pi2, pi3, pinew);
-        for(const auto & el : new_els)
-          mesh.AddVolumeElement(el);
-        mesh[ei1].Delete();
-      }
-    }
   }
-  mesh.Compress();
+  topt.Stop();
-  bad1 = mesh.CalcTotalBad (mp);
+
+  PrintMessage (5, cnt, " elements split");
+  (*testout) << "SplitImprove2 done" << "\n";
+
+  if(cnt>0)
+    mesh.Compress();
+  multithread.task = savetask;
 }
 
 

libsrc/meshing/improve3.hpp

@@ -39,7 +39,9 @@ public:
   void SplitImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
   void SplitImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
   double SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table<ElementIndex,PointIndex> & elementsonnode, Array<double> &elerrs, NgArray<INDEX_3> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only=false);
-  void SplitImprove2 (Mesh & mesh, OPTIMIZEGOAL goal);
+
+  void SplitImprove2 (Mesh & mesh);
+  double SplitImprove2Element (Mesh & mesh, ElementIndex ei, const Table<ElementIndex, PointIndex> & elements_of_point, const Array<double> & elerrs, bool check_only);
   
 
   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);

libsrc/meshing/meshfunc.cpp

@@ -684,7 +684,7 @@ namespace netgen
 	      {
 	      case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break;
 	      case 'd': optmesh.SplitImprove(mesh3d); break;
-	      case 'D': optmesh.SplitImprove2(mesh3d, OPT_QUALITY); break;
+	      case 'D': optmesh.SplitImprove2(mesh3d); break;
 	      case 's': optmesh.SwapImprove(mesh3d); break;
                 // case 'u': optmesh.SwapImproveSurface(mesh3d); break;
 	      case 't': optmesh.SwapImprove2(mesh3d); break;

tests/pytest/results.json

@@ -473,7 +473,7 @@
             "ne2d": 726,
             "ne3d": 2167,
             "quality_histogram": "[0, 4, 17, 35, 75, 117, 114, 112, 77, 51, 58, 86, 115, 177, 248, 293, 239, 204, 118, 27]",
-            "total_badness": 4176.9281305
+            "total_badness": 4176.9278168
         },
         {
             "angles_tet": [
@@ -1176,8 +1176,8 @@
             "ne1d": 432,
             "ne2d": 9544,
             "ne3d": 69846,
-            "quality_histogram": "[0, 0, 0, 0, 0, 0, 1, 7, 30, 71, 221, 652, 1560, 3667, 7029, 11233, 14630, 15524, 11810, 3411]",
+            "quality_histogram": "[0, 0, 0, 0, 0, 0, 1, 7, 30, 71, 221, 652, 1560, 3667, 7028, 11234, 14630, 15524, 11810, 3411]",
-            "total_badness": 85625.273734
+            "total_badness": 85625.275421
         }
     ],
     "ellipticcyl.geo": [
@@ -1364,6 +1364,53 @@
             "total_badness": 639.78974791
         }
     ],
+    "frame.step": [
+        {
+            "angles_tet": [
+                2.9095,
+                171.1
+            ],
+            "angles_trig": [
+                2.6092,
+                172.05
+            ],
+            "ne1d": 10108,
+            "ne2d": 30160,
+            "ne3d": 153012,
+            "quality_histogram": "[0, 3, 1, 3, 6, 20, 57, 149, 536, 1257, 2919, 5836, 10439, 16388, 21788, 26067, 26565, 22916, 14350, 3712]",
+            "total_badness": 202656.25887
+        },
+        {
+            "angles_tet": [
+                2.296,
+                175.61
+            ],
+            "angles_trig": [
+                1.8443,
+                175.57
+            ],
+            "ne1d": 5988,
+            "ne2d": 11102,
+            "ne3d": 29343,
+            "quality_histogram": "[3, 4, 5, 8, 14, 42, 121, 248, 691, 1040, 1542, 2504, 3118, 3920, 4331, 4281, 3366, 2421, 1367, 317]",
+            "total_badness": 43497.876838
+        },
+        {
+            "angles_tet": [
+                2.5792,
+                174.11
+            ],
+            "angles_trig": [
+                2.2053,
+                174.13
+            ],
+            "ne1d": 9622,
+            "ne2d": 23964,
+            "ne3d": 80995,
+            "quality_histogram": "[2, 14, 4, 20, 18, 40, 94, 225, 485, 1115, 2415, 4537, 7493, 10248, 12753, 13190, 12020, 9207, 5660, 1455]",
+            "total_badness": 111934.52308
+        }
+    ],
     "hinge.stl": [
         {
             "angles_tet": [
@@ -1652,9 +1699,9 @@
             ],
             "ne1d": 5886,
             "ne2d": 48052,
-            "ne3d": 178832,
+            "ne3d": 178844,
-            "quality_histogram": "[0, 0, 0, 0, 0, 5, 11, 72, 289, 828, 2308, 6210, 10937, 18854, 27390, 30600, 31363, 26926, 18496, 4543]",
+            "quality_histogram": "[0, 0, 0, 0, 0, 5, 11, 71, 288, 829, 2312, 6198, 10948, 18856, 27414, 30640, 31333, 26927, 18474, 4538]",
-            "total_badness": 233895.88826
+            "total_badness": 233920.54177
         },
         {
             "angles_tet": [
@@ -1722,8 +1769,8 @@
         },
         {
             "angles_tet": [
-                7.9601,
+                9.3063,
-                167.83
+                165.3
             ],
             "angles_trig": [
                 7.9174,
@@ -1731,9 +1778,9 @@
             ],
             "ne1d": 106,
             "ne2d": 610,
-            "ne3d": 1658,
+            "ne3d": 1659,
-            "quality_histogram": "[0, 1, 14, 48, 79, 156, 193, 158, 156, 122, 138, 142, 136, 108, 66, 40, 31, 41, 25, 4]",
+            "quality_histogram": "[0, 1, 13, 49, 81, 155, 190, 149, 159, 132, 135, 145, 137, 107, 66, 37, 33, 40, 26, 4]",
-            "total_badness": 4098.9846426
+            "total_badness": 4094.4605262
         },
         {
             "angles_tet": [
@@ -1783,7 +1830,7 @@
         {
             "angles_tet": [
                 18.203,
-                145.26
+                145.38
             ],
             "angles_trig": [
                 17.821,
@@ -1791,9 +1838,9 @@
             ],
             "ne1d": 418,
             "ne2d": 5968,
-            "ne3d": 101104,
+            "ne3d": 101113,
-            "quality_histogram": "[0, 0, 0, 0, 0, 0, 4, 7, 52, 103, 348, 992, 2548, 5550, 10212, 16062, 20713, 22250, 16917, 5346]",
+            "quality_histogram": "[0, 0, 0, 0, 0, 0, 4, 7, 51, 102, 349, 993, 2550, 5563, 10196, 16090, 20698, 22258, 16911, 5341]",
-            "total_badness": 124141.43144
+            "total_badness": 124155.81178
         }
     ],
     "ortho.geo": [
@@ -2144,9 +2191,9 @@
             ],
             "ne1d": 2992,
             "ne2d": 23322,
-            "ne3d": 281901,
+            "ne3d": 281896,
-            "quality_histogram": "[4, 10, 12, 10, 9, 21, 29, 61, 95, 247, 747, 2145, 5539, 13549, 27678, 44557, 59948, 64037, 48291, 14912]",
+            "quality_histogram": "[4, 10, 12, 10, 9, 21, 29, 61, 95, 246, 747, 2146, 5540, 13546, 27675, 44558, 59947, 64037, 48291, 14912]",
-            "total_badness": 344516.23097
+            "total_badness": 344508.9779
         }
     ],
     "revolution.geo": [
@@ -2241,6 +2288,53 @@
             "total_badness": 244286.77424
         }
     ],
+    "screw.step": [
+        {
+            "angles_tet": [
+                15.139,
+                148.36
+            ],
+            "angles_trig": [
+                17.363,
+                140.59
+            ],
+            "ne1d": 400,
+            "ne2d": 1436,
+            "ne3d": 2342,
+            "quality_histogram": "[0, 0, 0, 0, 0, 0, 16, 61, 92, 165, 209, 246, 307, 280, 262, 272, 182, 148, 84, 18]",
+            "total_badness": 3718.1755695
+        },
+        {
+            "angles_tet": [
+                21.55,
+                146.38
+            ],
+            "angles_trig": [
+                17.221,
+                126.09
+            ],
+            "ne1d": 528,
+            "ne2d": 2792,
+            "ne3d": 8129,
+            "quality_histogram": "[0, 0, 0, 0, 0, 1, 4, 8, 35, 96, 188, 298, 537, 817, 1057, 1323, 1446, 1284, 798, 237]",
+            "total_badness": 10753.086327
+        },
+        {
+            "angles_tet": [
+                20.515,
+                144.03
+            ],
+            "angles_trig": [
+                24.891,
+                120.48
+            ],
+            "ne1d": 666,
+            "ne2d": 4922,
+            "ne3d": 31540,
+            "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 1, 5, 31, 92, 290, 707, 1762, 3221, 4997, 6712, 6966, 5146, 1610]",
+            "total_badness": 38689.280913
+        }
+    ],
     "sculpture.geo": [
         {
             "angles_tet": [
@@ -2914,7 +3008,7 @@
             "ne2d": 692,
             "ne3d": 2737,
             "quality_histogram": "[17, 200, 365, 335, 363, 301, 234, 187, 154, 143, 106, 84, 56, 48, 38, 45, 27, 19, 12, 3]",
-            "total_badness": 13234.68014
+            "total_badness": 13234.755766
         },
         {
             "angles_tet": [
@@ -2985,13 +3079,13 @@
             ],
             "angles_trig": [
                 14.916,
-                130.79
+                132.02
             ],
             "ne1d": 690,
             "ne2d": 1684,
-            "ne3d": 5186,
+            "ne3d": 5177,
-            "quality_histogram": "[0, 0, 1, 0, 0, 9, 26, 34, 106, 192, 291, 357, 463, 570, 662, 695, 628, 542, 471, 139]",
+            "quality_histogram": "[0, 0, 1, 0, 1, 8, 27, 37, 108, 191, 285, 369, 461, 565, 670, 690, 621, 536, 462, 145]",
-            "total_badness": 7461.9914431
+            "total_badness": 7461.1502455
         },
         {
             "angles_tet": [
@@ -3019,9 +3113,9 @@
             ],
             "ne1d": 512,
             "ne2d": 874,
-            "ne3d": 2382,
+            "ne3d": 2381,
-            "quality_histogram": "[0, 0, 0, 3, 9, 13, 41, 68, 124, 140, 196, 215, 305, 390, 343, 237, 128, 98, 47, 25]",
+            "quality_histogram": "[0, 0, 0, 3, 9, 13, 41, 68, 124, 140, 196, 214, 302, 390, 345, 237, 128, 98, 47, 26]",
-            "total_badness": 3929.8055104
+            "total_badness": 3927.0434195
         },
         {
             "angles_tet": [
@@ -3030,13 +3124,13 @@
             ],
             "angles_trig": [
                 14.916,
-                130.79
+                132.02
             ],
             "ne1d": 690,
             "ne2d": 1684,
-            "ne3d": 5110,
+            "ne3d": 5095,
-            "quality_histogram": "[0, 0, 1, 0, 0, 4, 18, 33, 104, 176, 266, 350, 443, 555, 682, 703, 613, 551, 468, 143]",
+            "quality_histogram": "[0, 0, 1, 0, 0, 4, 19, 34, 101, 181, 263, 354, 439, 548, 689, 696, 611, 547, 467, 141]",
-            "total_badness": 7297.2366495
+            "total_badness": 7282.7477612
         },
         {
             "angles_tet": [
@@ -3049,9 +3143,9 @@
             ],
             "ne1d": 1050,
             "ne2d": 3812,
-            "ne3d": 18010,
+            "ne3d": 18003,
-            "quality_histogram": "[0, 0, 0, 0, 0, 0, 3, 15, 34, 61, 181, 572, 1425, 2192, 2302, 2698, 2704, 2739, 2390, 694]",
+            "quality_histogram": "[0, 0, 0, 0, 0, 0, 3, 15, 34, 61, 181, 573, 1428, 2189, 2298, 2700, 2707, 2735, 2389, 690]",
-            "total_badness": 23478.48161
+            "total_badness": 23471.146878
         },
         {
             "angles_tet": [
@@ -3064,26 +3158,26 @@
             ],
             "ne1d": 1722,
             "ne2d": 10042,
-            "ne3d": 84793,
+            "ne3d": 84812,
-            "quality_histogram": "[0, 0, 0, 0, 0, 2, 50, 1428, 718, 373, 701, 1171, 2465, 5462, 8878, 13214, 16426, 16935, 12864, 4106]",
+            "quality_histogram": "[0, 0, 0, 0, 0, 2, 49, 1423, 720, 374, 704, 1174, 2454, 5477, 8890, 13211, 16429, 16935, 12870, 4100]",
-            "total_badness": 108481.65242
+            "total_badness": 108503.84867
         }
     ],
     "twobricks.geo": [
         {
             "angles_tet": [
-                25.655,
+                29.453,
-                137.37
+                134.56
             ],
             "angles_trig": [
-                23.972,
+                26.574,
-                121.23
+                91.538
             ],
             "ne1d": 72,
             "ne2d": 50,
-            "ne3d": 43,
+            "ne3d": 36,
-            "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 5, 6, 3, 16, 3, 6, 1, 1, 0, 0]",
+            "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 8, 2, 18, 2, 4, 0, 0, 0, 0]",
-            "total_badness": 67.089294311
+            "total_badness": 55.618194358
         },
         {
             "angles_tet": [
@@ -3117,18 +3211,18 @@
         },
         {
             "angles_tet": [
-                25.655,
+                29.453,
-                137.37
+                134.56
             ],
             "angles_trig": [
-                23.971,
+                26.574,
-                121.23
+                91.538
             ],
             "ne1d": 72,
             "ne2d": 50,
-            "ne3d": 43,
+            "ne3d": 36,
-            "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 5, 6, 3, 16, 3, 6, 1, 1, 0, 0]",
+            "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 8, 2, 18, 2, 4, 0, 0, 0, 0]",
-            "total_badness": 67.089294246
+            "total_badness": 55.618194358
         },
         {
             "angles_tet": [
@@ -3164,18 +3258,18 @@
     "twocubes.geo": [
         {
             "angles_tet": [
-                25.655,
+                29.453,
-                137.37
+                134.56
             ],
             "angles_trig": [
-                23.972,
+                26.574,
-                121.23
+                91.538
             ],
             "ne1d": 72,
             "ne2d": 50,
-            "ne3d": 43,
+            "ne3d": 36,
-            "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 5, 6, 3, 16, 3, 6, 1, 1, 0, 0]",
+            "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 8, 2, 18, 2, 4, 0, 0, 0, 0]",
-            "total_badness": 67.089294311
+            "total_badness": 55.618194358
         },
         {
             "angles_tet": [
@@ -3209,18 +3303,18 @@
         },
         {
             "angles_tet": [
-                25.655,
+                29.453,
-                137.37
+                134.56
             ],
             "angles_trig": [
-                23.971,
+                26.574,
-                121.23
+                91.538
             ],
             "ne1d": 72,
             "ne2d": 50,
-            "ne3d": 43,
+            "ne3d": 36,
-            "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 5, 6, 3, 16, 3, 6, 1, 1, 0, 0]",
+            "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 8, 2, 18, 2, 4, 0, 0, 0, 0]",
-            "total_badness": 67.089294246
+            "total_badness": 55.618194358
         },
         {
             "angles_tet": [