some more bugfixing

libsrc/meshing/boundarylayer.cpp

@@ -530,8 +530,12 @@ struct GrowthVectorLimiter {
           if(!pts.Contains(other))
             pts.Append(other);
       }
-      if(pts.Size() != 2)
+      if(pts.Size() != 2) {
+        cout << "getEdgeTangent pi = " << pi << ", edgenr = " << edgenr << endl;
+        for(auto segi : topo.GetVertexSegments(pi))
+          cout  << mesh[segi] << endl;
         throw Exception("Something went wrong in getEdgeTangent!");
+      }
       tangent = mesh[pts[1]] - mesh[pts[0]];
       return tangent.Normalize();
   }
@@ -1380,12 +1384,16 @@ struct GrowthVectorLimiter {
         new_max_edge_nr = seg.edgenr;
 
     auto getGW = [&] (PointIndex pi) -> Vec<3>& {
+        static Vec<3> zero(0.,0.,0.);
+        if(growth_vector_map.count(pi))
           return *get<0>(growth_vector_map[pi]);
+        zero = {0.,0.,0.};
+        return zero;
     };
   // cout << "edge range " << max_edge_nr << ", " << new_max_edge_nr << endl;
 
     // interpolate tangential component of growth vector along edge
-    for(auto edgenr : Range(max_edge_nr, new_max_edge_nr))
+    for(auto edgenr : Range(max_edge_nr+1, new_max_edge_nr))
       {
     // cout << "SEARCH EDGE " << edgenr +1 << endl;
         // if(!is_edge_moved[edgenr+1]) continue;
@@ -1492,22 +1500,25 @@ struct GrowthVectorLimiter {
           //     gt2 = 0.;
           // }
 
+        // cout << "edgenr " << edgenr << endl;
+        // cout << "points " << endl << points << endl;
+
         double len = 0.;
-        for(size_t i = 1; i < points.Size()-1; i++)
+        for(auto i : IntRange(1, points.Size()-1))
           {
             auto pi = points[i];
             len += (mesh[pi] - mesh[points[i-1]]).Length();
             auto t = getEdgeTangent(pi, edgenr);
             auto lam = len/edge_len;
             auto interpol = (1-lam) * (gt1 * t) * t + lam * (gt2 * t) * t;
-      if(pi==89) {
+      // if(pi==89) {
-        cout << "points " << points << endl;
+      //   cout << "points " << points << endl;
-        cout << "INTERPOL" << len << ',' << t << ',' << lam << ',' << interpol << endl;
+      //   cout << "INTERPOL" << len << ',' << t << ',' << lam << ',' << interpol << endl;
-        cout << gt1 << endl;
+      //   cout << gt1 << endl;
-        cout << gt2 << endl;
+      //   cout << gt2 << endl;
-        cout << getGW(pi) << endl;
+      //   cout << getGW(pi) << endl;
 
-      }
+      // }
             getGW(pi) += interpol;
           }
       }
@@ -1537,6 +1548,11 @@ struct GrowthVectorLimiter {
           growth_vector_map[pi_new] = { &growth_vector, params.heights[i] };
           if(special_boundary_points.count(pi) == 0)
             mesh.GetIdentifications().Add(pi_last, pi_new, identnr);
+          else
+          {
+            cout << "add locked point " << pi_new << endl;
+            mesh.AddLockedPoint(pi_new);
+          }
           pi_last = pi_new;
         }
     };

libsrc/meshing/meshclass.cpp

@@ -4081,6 +4081,7 @@ namespace netgen
           }
         else
           {
+            cout << "unused point " << pi << endl;
             op2np[pi].Invalidate(); 
           }
       

libsrc/meshing/meshfunc.cpp

@@ -50,6 +50,7 @@ namespace netgen
           ret[0].mp = mp;
           return ret;
       }
+  cout << "divide mesh" << endl;
       ret.SetSize(num_domains);
 
       Array<Array<PointIndex, PointIndex>> ipmap;
@@ -71,6 +72,7 @@ namespace netgen
 
           m.SetLocalH(mesh.GetLocalH());
 
+    cout << "set imap size " << i << " to " << num_points << endl;
           ipmap[i].SetSize(num_points);
           ipmap[i] = PointIndex::INVALID;
           m.SetDimension( mesh.GetDimension() );
@@ -155,6 +157,8 @@ namespace netgen
           auto & imap = ipmap[i];
           auto nmax = identifications.GetMaxNr ();
           auto & m_ident = m.GetIdentifications();
+          cout << "imap " << imap << endl;
+          cout << imap.Size() << endl;
 
           for (auto & sel : m.SurfaceElements())
             for(auto & pi : sel.PNums())
@@ -171,6 +175,7 @@ namespace netgen
 
               for(auto pair : pairs)
               {
+        // cout << "get pair " << pair[0] << ',' << pair[1] << endl;
                   auto pi0 = imap[pair[0]];
                   auto pi1 = imap[pair[1]];
                   if(!pi0.IsValid() || !pi1.IsValid())
@@ -289,6 +294,7 @@ namespace netgen
     for (int qstep = 0; qstep <= 3; qstep++)
      {
        if (qstep == 0 && !mp.try_hexes) continue;
+       if (qstep == 1) continue;
        
        if (mesh.HasOpenQuads())
          {
@@ -301,6 +307,7 @@ namespace netgen
                rulep = hexrules;
                break;
              case 1:
+               cout << "Apply prismrules " << endl;
                rulep = prismrules2;
                break;
              case 2: // connect pyramid to triangle
@@ -428,7 +435,11 @@ namespace netgen
          if (cntsteps > mp.maxoutersteps) 
          {
            if(debugparam.write_mesh_on_error)
+           {
              md.mesh->Save("meshing_error_domain_"+ToString(md.domain)+".vol.gz");
+             if(mesh.GetNOpenElements())
+               GetOpenElements(*md.mesh, md.domain)->Save("meshing_error_rest_" + ToString(md.domain)+".vol.gz");
+           }
            throw NgException ("Stop meshing since too many attempts in domain " + ToString(md.domain));
          }
 
@@ -532,6 +543,7 @@ namespace netgen
          md[0].mesh.release();
          return;
      }
+  cout << "divide mesh" << endl;
 
      mesh.VolumeElements().DeleteAll();
      for(auto & m_ : md)
@@ -584,7 +596,9 @@ namespace netgen
   {
     static Timer t("MeshVolume"); RegionTimer reg(t);
 
+  cout << "initial compress " << endl;
      mesh3d.Compress();
+  cout << "initial compress done" << endl;
 
 
 

libsrc/visualization/mvdraw.hpp

@@ -242,6 +242,9 @@ namespace netgen
     ngcore::INT<2> Project(Point<3> p);
   };
 
+  void DrawElement(const Mesh & mesh, SurfaceElementIndex sei);
+  void DrawElement(const Mesh & mesh, ElementIndex sei);
+
   NGGUI_API extern VisualSceneMesh vsmesh;