return invalid surfaceindex (index is 0 based)

Don't check for intersecting mapped trigs if only volume elements are
added (and no mapped surface elements).

In case of only inserted volume elements, we don't care about the limit
at special points.

cmake/SuperBuild.cmake

@@ -57,7 +57,7 @@ set_vars(SUBPROJECT_CMAKE_ARGS CMAKE_C_COMPILER)
 set_vars(SUBPROJECT_CMAKE_ARGS CMAKE_CXX_COMPILER)
 set_vars(SUBPROJECT_CMAKE_ARGS CMAKE_BUILD_TYPE)
 
-set(SUBPROJECT_CMAKE_ARGS "${SUBPROJECT_CMAKE_ARGS};-DCMAKE_POSITION_INDEPENDENT_CODE=ON" CACHE INTERNAL "")
+set(SUBPROJECT_CMAKE_ARGS "${SUBPROJECT_CMAKE_ARGS};-DCMAKE_POSITION_INDEPENDENT_CODE=ON;-DCMAKE_POLICY_VERSION_MINIMUM=3.5" CACHE INTERNAL "")
 
 if(USE_CCACHE)
   find_program(CCACHE_FOUND NAMES ccache ccache.bat)

cmake/external_projects/tcltk.cmake

@@ -109,6 +109,7 @@ if(APPLE)
            -DCMAKE_INSTALL_PREFIX=${CMAKE_INSTALL_PREFIX}/Contents/MacOS
 	   -DTCL_INCLUDE_PATH=${CMAKE_INSTALL_PREFIX}/Contents/Frameworks/Tcl.framework/Headers
 	   -DTK_INCLUDE_PATH=${CMAKE_INSTALL_PREFIX}/Contents/Frameworks/Tk.framework/Headers
+     -DCMAKE_POLICY_VERSION_MINIMUM=3.5
     ${SUBPROJECT_ARGS}
   )
 

libsrc/csg/genmesh.cpp

@@ -718,7 +718,7 @@ namespace netgen
 
 	    mesh -> LoadLocalMeshSize (mparam.meshsizefilename);
             for (auto mspnt : mparam.meshsize_points)
-              mesh -> RestrictLocalH (mspnt.pnt, mspnt.h);
+              mesh -> RestrictLocalH (mspnt.pnt, mspnt.h, mspnt.layer);
 	  }
 
 	spoints.SetSize(0);

libsrc/include/nginterface_v2_impl.hpp

@@ -72,13 +72,17 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<0> (size_t nr) const
   ret.facets.base = POINTINDEX_BASE;
   ret.facets.ptr = (int*)&el.pnum;
 
+  /*
   if (mesh->GetDimension() == 1)
     ret.mat = *(mesh->GetBCNamePtr(el.index-1));
   else if (mesh->GetDimension() == 2)
     ret.mat = *(mesh->GetCD2NamePtr(el.index-1));
   else
     ret.mat = *(mesh->GetCD3NamePtr(el.index-1));
-  
+  */
+  ret.mat = mesh->GetRegionName(0, el.index);
+    
+  ret.is_curved = false;
   return ret;
 }
 
@@ -96,6 +100,8 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (size_t nr) const
     ret.index = el.edgenr;
   else
     ret.index = el.si;
+
+  /*
   if (mesh->GetDimension() == 2)
     ret.mat = *(mesh->GetBCNamePtr(el.si-1));
   else
@@ -105,6 +111,8 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (size_t nr) const
       else
         ret.mat = *(mesh->GetMaterialPtr(el.si));
     }
+  */
+  ret.mat = mesh->GetRegionName(1, ret.index);
 
   ret.points.num = el.GetNP();
   ret.points.ptr = (int*)&(el[0]);

libsrc/interface/nginterface.cpp

@@ -651,14 +651,14 @@ int Ng_FindElementOfPoint (double * p, double * lami, int build_searchtree,
     {
       Point3d p3d(p[0], p[1], p[2]);
       ind = 
-	mesh->GetElementOfPoint(p3d, lami, dummy, build_searchtree != 0);
+	mesh->GetElementOfPoint(p3d, lami, dummy, build_searchtree != 0) + 1;
     }
   else
     {
       double lam3[3];
       Point3d p2d(p[0], p[1], 0);
       ind = 
-	mesh->GetElementOfPoint(p2d, lam3, dummy, build_searchtree != 0);
+	mesh->GetSurfaceElementOfPoint(p2d, lam3, dummy, build_searchtree != 0) + 1;
 
       if (ind > 0)
 	{
@@ -697,7 +697,7 @@ int Ng_FindSurfaceElementOfPoint (double * p, double * lami, int build_searchtre
     {
       Point3d p3d(p[0], p[1], p[2]);
       ind = 
-	mesh->GetSurfaceElementOfPoint(p3d, lami, dummy, build_searchtree != 0);
+	mesh->GetSurfaceElementOfPoint(p3d, lami, dummy, build_searchtree != 0) + 1;
     }
   else
     {

libsrc/interface/nginterface_v2.cpp

@@ -1011,68 +1011,28 @@ namespace netgen
    int * const indices, int numind) const
 
   {
-    switch (mesh->GetDimension())
+    Point<3> p(hp[0], 0., 0.);
+    if(mesh->GetDimension() > 1)
+      p[1] = hp[1];
+    if(mesh->GetDimension() == 3)
+      p[2] = hp[2];
+
+    for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
       {
-      case 1:
-        {
-          Point<3> p(hp[0], 0,0);
-          for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
-            {
-              auto & seg = (*mesh)[si];
-              Point<3> p1 = (*mesh)[seg[0]];
-              Point<3> p2 = (*mesh)[seg[1]];
-              double lam = (p(0)-p1(0)) / (p2(0)-p1(0));
-              if (lam >= -1e-10 && lam <= 1+1e-10)
-                {
-                  lami[0] = 1-lam;
-                  return si;
-                }
-            }
-        }
-        break;
-      case 2:
-        {
-          Point<3> p(hp[0], hp[1],0);
-          try
-            {
-              auto ind = mesh->GetSurfaceElementOfPoint(p, lami, nullptr,
-                                                        build_searchtree);
-              return ind - 1;
-            }
-          catch(const NgException & e) // quads not implemented curved yet
-            {
-          for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
-            {
-              auto & seg = (*mesh)[si];
-              Point<3> p1 = (*mesh)[seg[0]];
-              Point<3> p2 = (*mesh)[seg[1]];
-              double lam;
-              double r;
-              if (fabs(p2[0]-p1[0]) >= fabs(p2[1]-p1[1]))
-                {
-                  lam = (p[0]-p1[0])/(p2[0]-p1[0]);
-                  r = p[1] - p1[1] - lam*(p2[1]-p1[1]);
-                }
-              else
-                {
-                  lam = (p[1]-p1[1])/(p2[1]-p1[1]);
-                  r = p[0] - p1[0] - lam*(p2[0]-p1[0]);
-                }
-              if ( lam >= -1e-10 && lam <= 1+1e-10 && fabs(r) <= 1e-10 )
-                {
-                  lami[0] = 1-lam;
-                  return si;
-                }
-            }
-            }
-        }
-        break;
-      case 3:
-      default:
-        throw Exception("FindElementOfPoint<1> only implemented for mesh-dimension 1 and 2!");
-        break;
+        auto & seg = (*mesh)[si];
+        Point<3> p1 = (*mesh)[seg[0]];
+        Point<3> p2 = (*mesh)[seg[1]];
+        Vec<3> v1 = p2-p1;
+        Vec<3> v2 = p-p1;
+        double lam = v1*v2 / v1.Length2();
+        double lam2 = (v2 - lam * v1).Length() / v1.Length();
+
+        if (lam >= -1e-10 && lam <= 1+1e-10 && lam2 < 1e-10)
+          {
+            lami[0] = 1-lam;
+            return si;
+          }
       }
- 
     return -1;
   }
 
@@ -1083,37 +1043,26 @@ namespace netgen
    int * const indices, int numind) const
 
   {
-    NgArray<int> dummy(numind);
-    for (int i = 0; i < numind; i++) dummy[i] = indices[i]+1;
-    
+    Point<3> pp(p[0], p[1], 0.);
+    if(mesh->GetDimension() == 3)
+      pp[2] = p[2];
+    FlatArray<int> ind(numind, indices);
     double lam3[3];
-    int ind;
-
-    if (mesh->GetDimension() == 2)
+    auto elnr = mesh->GetSurfaceElementOfPoint(pp, lam3, ind, build_searchtree);
+    if(elnr.IsValid())
       {
-        Point<3> p2d(p[0], p[1], 0);
-        ind = mesh->GetElementOfPoint(p2d, lam3, &dummy, build_searchtree);
-      }
-    else
-      {
-        Point3d p3d(p[0], p[1], p[2]);
-        ind = mesh->GetSurfaceElementOfPoint(p3d, lam3, &dummy, build_searchtree);
-      }
-    
-    if (ind > 0)
-      {
-        if(mesh->SurfaceElement(ind).GetType()==QUAD || mesh->SurfaceElement(ind).GetType()==TRIG6)
+        if((*mesh)[elnr].GetType() == QUAD || (*mesh)[elnr].GetType() == TRIG6)
           {
             lami[0] = lam3[0];
             lami[1] = lam3[1];
           }
-        else 
+        else
           {
             lami[0] = 1-lam3[0]-lam3[1];
             lami[1] = lam3[0];
           }
       }
-    return ind-1;
+    return elnr;
   }
 
 
@@ -1124,13 +1073,9 @@ namespace netgen
    int * const indices, int numind) const
 
   {
-    NgArray<int> dummy(numind);
-    for (int i = 0; i < numind; i++) dummy[i] = indices[i]+1;
-    
-    Point<3> p3d(p[0], p[1], p[2]);
-    int ind = 
-      mesh->GetElementOfPoint(p3d, lami, &dummy, build_searchtree);
-    return ind-1;
+    Point<3> pp(p[0], p[1], p[2]);
+    FlatArray<int> ind(numind, indices);
+    return mesh->GetElementOfPoint(pp, lami, ind, build_searchtree);
   }
 
   void Ngx_Mesh :: Curve (int order)

libsrc/interface/writeOpenFOAM15x.cpp

@@ -607,7 +607,7 @@ namespace netgen
       const_cast<Mesh&> (mesh).Compress();
       const_cast<Mesh&> (mesh).CalcSurfacesOfNode();
       const_cast<Mesh&> (mesh).RebuildSurfaceElementLists();
-      const_cast<Mesh&> (mesh).BuildElementSearchTree();
+      const_cast<Mesh&> (mesh).BuildElementSearchTree(3);
 
 
       int np = mesh.GetNP();

libsrc/interface/writefluent.cpp

@@ -66,7 +66,7 @@ void WriteFluentFormat (const Mesh & mesh,
   int i2, j2;
   NgArray<INDEX_3> surfaceelp;
   NgArray<int> surfaceeli;
-  NgArray<int> locels;
+  Array<ElementIndex> locels;
 
   //no cells=no tets
   //no faces=2*tets
@@ -79,7 +79,7 @@ void WriteFluentFormat (const Mesh & mesh,
   snprintf (str, size(str), "(13 (4 1 %x 2 3)(",noverbface); //hexadecimal!!!
   outfile << str << endl;
 
-  const_cast<Mesh&> (mesh).BuildElementSearchTree();
+  const_cast<Mesh&> (mesh).BuildElementSearchTree(3);
 
   for (i = 1; i <= ne; i++)
     {
@@ -117,12 +117,9 @@ void WriteFluentFormat (const Mesh & mesh,
 	  int eli2 = 0;
 	  int stopsig = 0;
 	      
-	  for (i2 = 1; i2 <= nel; i2++)
+	  for (auto locind : locels)
 	    {
-	      locind = locels.Get(i2);
-	      //cout << "  locind=" << locind << endl;
-
-	      Element el2 = mesh.VolumeElement(locind);
+	      Element el2 = mesh[locind];
 	      //if (inverttets)
 	      //  el2.Invert();
 
@@ -130,7 +127,7 @@ void WriteFluentFormat (const Mesh & mesh,
 		{
 		  el2.GetFace(j2, face2);
 
-		  if (face2.HasFace(face)) {eli2 = locind; stopsig = 1; break;}
+		  if (face2.HasFace(face)) {eli2 = locind+1; stopsig = 1; break;}
 		}
 	      if (stopsig) break;
 	    }

libsrc/meshing/basegeom.cpp

@@ -474,7 +474,7 @@ namespace netgen
       }
 
     for(const auto& mspnt : mparam.meshsize_points)
-      mesh.RestrictLocalH(mspnt.pnt, mspnt.h);
+      mesh.RestrictLocalH(mspnt.pnt, mspnt.h, mspnt.layer);
 
     mesh.LoadLocalMeshSize(mparam.meshsizefilename);
   }
@@ -1220,6 +1220,7 @@ namespace netgen
       {
         PrintMessage(3, "Optimization step ", i);
         meshopt.SetFaceIndex(k+1);
+        meshopt.SetMetricWeight (mparam.elsizeweight);
         int innerstep = 0;
         for(auto optstep : mparam.optimize2d)
           {
@@ -1314,6 +1315,13 @@ namespace netgen
     if(multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHEDGES)
       return 0;
 
+    if(dimension == 1)
+      {
+        FinalizeMesh(*mesh);
+        mesh->SetDimension(1);
+        return 0;
+      }
+
     if (mparam.perfstepsstart <= MESHCONST_MESHSURFACE)
       {
         MeshSurface(*mesh, mparam);

libsrc/meshing/bisect.cpp

@@ -4061,6 +4061,8 @@ namespace netgen
           {
             PointIndices<2> oi2(identmap[i2[0]], 
                                 identmap[i2[1]]);
+            if((!oi2[0].IsValid()) || (!oi2[1].IsValid()))
+              continue;
             oi2.Sort();
             if (cutedges.Used (oi2))
               {

libsrc/meshing/boundarylayer.cpp

@@ -14,7 +14,7 @@ namespace netgen
 
 struct SpecialPointException : public Exception
 {
-  SpecialPointException()
+  SpecialPointException ()
     : Exception("") {}
 };
 
@@ -101,14 +101,14 @@ Vec<3> CalcGrowthVector (FlatArray<Vec<3>> ns)
   return gw;
 }
 
-SpecialBoundaryPoint ::GrowthGroup ::GrowthGroup(FlatArray<int> faces_,
-                                                 FlatArray<Vec<3>> normals)
+SpecialBoundaryPoint ::GrowthGroup ::GrowthGroup (FlatArray<int> faces_,
+                                                  FlatArray<Vec<3>> normals)
 {
   faces = faces_;
   growth_vector = CalcGrowthVector(normals);
 }
 
-SpecialBoundaryPoint ::SpecialBoundaryPoint(
+SpecialBoundaryPoint ::SpecialBoundaryPoint (
   const std::map<int, Vec<3>>& normals)
 {
   // find opposing face normals
@@ -146,7 +146,7 @@ SpecialBoundaryPoint ::SpecialBoundaryPoint(
   growth_groups.Append(GrowthGroup(g2_faces, normals2));
 }
 
-Vec<3> BoundaryLayerTool ::getEdgeTangent(PointIndex pi, int edgenr, FlatArray<Segment*> segs)
+Vec<3> BoundaryLayerTool ::getEdgeTangent (PointIndex pi, int edgenr, FlatArray<Segment*> segs)
 {
   Vec<3> tangent = 0.0;
   ArrayMem<PointIndex, 2> pts;
@@ -171,7 +171,7 @@ Vec<3> BoundaryLayerTool ::getEdgeTangent(PointIndex pi, int edgenr, FlatArray<S
   return tangent.Normalize();
 }
 
-void BoundaryLayerTool ::LimitGrowthVectorLengths()
+void BoundaryLayerTool ::LimitGrowthVectorLengths ()
 {
   static Timer tall("BoundaryLayerTool::LimitGrowthVectorLengths");
   RegionTimer rtall(tall);
@@ -220,7 +220,7 @@ bool HaveSingleSegments (const Mesh& mesh)
 
 // duplicates segments (and sets seg.si accordingly) to have a unified data
 // structure for all geometry types
-void BuildSegments (Mesh& mesh, bool have_single_segments, Array<Segment> & segments, Array<Segment> & free_segments)
+void BuildSegments (Mesh& mesh, bool have_single_segments, Array<Segment>& segments, Array<Segment>& free_segments)
 {
   // auto& topo = mesh.GetTopology();
 
@@ -234,11 +234,11 @@ void BuildSegments (Mesh& mesh, bool have_single_segments, Array<Segment> & segm
           free_segments.Append(seg);
           continue;
         }
-      if(!have_single_segments)
-      {
-        segments.Append(seg);
-        continue;
-      }
+      if (!have_single_segments)
+        {
+          segments.Append(seg);
+          continue;
+        }
       mesh.GetTopology().GetSegmentSurfaceElements(segi + 1, surf_els);
       for (auto seli : surf_els)
         {
@@ -284,8 +284,8 @@ void MergeAndAddSegments (Mesh& mesh, FlatArray<Segment> segments, FlatArray<Seg
     addSegment(seg);
 }
 
-BoundaryLayerTool::BoundaryLayerTool(Mesh& mesh_,
-                                     const BoundaryLayerParameters& params_)
+BoundaryLayerTool::BoundaryLayerTool (Mesh& mesh_,
+                                      const BoundaryLayerParameters& params_)
   : mesh(mesh_), topo(mesh_.GetTopology()), params(params_)
 {
   static Timer timer("BoundaryLayerTool::ctor");
@@ -318,7 +318,7 @@ BoundaryLayerTool::BoundaryLayerTool(Mesh& mesh_,
     si_map[i] = i;
 }
 
-void BoundaryLayerTool ::CreateNewFaceDescriptors()
+void BoundaryLayerTool ::CreateNewFaceDescriptors ()
 {
   surfacefacs.SetSize(nfd_old + 1);
   surfacefacs = 0.0;
@@ -360,7 +360,7 @@ void BoundaryLayerTool ::CreateNewFaceDescriptors()
       throw Exception("Surface " + to_string(si) + " is not a boundary of the domain to be grown into!");
 }
 
-void BoundaryLayerTool ::CreateFaceDescriptorsSides()
+void BoundaryLayerTool ::CreateFaceDescriptorsSides ()
 {
   if (insert_only_volume_elements)
     return;
@@ -400,7 +400,7 @@ void BoundaryLayerTool ::CreateFaceDescriptorsSides()
     }
 }
 
-void BoundaryLayerTool ::CalculateGrowthVectors()
+void BoundaryLayerTool ::CalculateGrowthVectors ()
 {
   growthvectors.SetSize(np);
   growthvectors = 0.;
@@ -457,7 +457,7 @@ void BoundaryLayerTool ::CalculateGrowthVectors()
 }
 
 Array<Array<pair<SegmentIndex, int>>, SegmentIndex>
-BoundaryLayerTool ::BuildSegMap()
+BoundaryLayerTool ::BuildSegMap ()
 {
   // Bit array to keep track of segments already processed
   BitArray segs_done(nseg + 1);
@@ -536,7 +536,7 @@ BoundaryLayerTool ::BuildSegMap()
   return segmap;
 }
 
-BitArray BoundaryLayerTool ::ProjectGrowthVectorsOnSurface()
+BitArray BoundaryLayerTool ::ProjectGrowthVectorsOnSurface ()
 {
   BitArray in_surface_direction(nfd_old + 1);
   in_surface_direction.Clear();
@@ -596,7 +596,7 @@ BitArray BoundaryLayerTool ::ProjectGrowthVectorsOnSurface()
   return in_surface_direction;
 }
 
-void BoundaryLayerTool ::InsertNewElements(
+void BoundaryLayerTool ::InsertNewElements (
   FlatArray<Array<pair<SegmentIndex, int>>, SegmentIndex> segmap,
   const BitArray& in_surface_direction)
 {
@@ -879,61 +879,61 @@ void BoundaryLayerTool ::InsertNewElements(
   auto p2el = mesh.CreatePoint2ElementTable();
   for (SurfaceElementIndex si = 0; si < nse; si++)
     {
-      // copy because surfaceels array will be resized!
       const auto sel = mesh[si];
-      if (moved_surfaces.Test(sel.GetIndex()))
+      const auto iface = sel.GetIndex();
+
+      if (moved_surfaces.Test(iface))
         {
-          Array<PointIndex> points(sel.PNums());
-          if (surfacefacs[sel.GetIndex()] > 0)
+          const auto np = sel.GetNP();
+          ArrayMem<PointIndex, 4> points(sel.PNums());
+          if (surfacefacs[iface] > 0)
             Swap(points[0], points[2]);
           ArrayMem<int, 4> groups(points.Size());
           for (auto i : Range(points))
-            groups[i] = getClosestGroup(sel[i], si);
+            groups[i] = getClosestGroup(points[i], si);
           bool add_volume_element = true;
-          for (auto pi : sel.PNums())
+          for (auto pi : points)
             if (numGroups(pi) > 1)
               add_volume_element = false;
+
+          Element el(2 * np);
+          el.PNums().Range(np, 2 * np) = points;
+          auto new_index = new_mat_nrs[iface];
+          if (new_index == -1)
+            throw Exception("Boundary " + ToString(iface) + " with name " + mesh.GetBCName(iface - 1) + " extruded, but no new material specified for it!");
+          el.SetIndex(new_index);
+
           for (auto j : Range(par_heights))
             {
-              auto eltype = points.Size() == 3 ? PRISM : HEX;
-              Element el(eltype);
-              for (auto i : Range(points))
-                el[i] = points[i];
-              for (auto i : Range(points))
-                points[i] = newPoint(sel.PNums()[i], j, groups[i]);
-              if (surfacefacs[sel.GetIndex()] > 0)
-                Swap(points[0], points[2]);
-              for (auto i : Range(points))
-                el[sel.PNums().Size() + i] = points[i];
-              auto new_index = new_mat_nrs[sel.GetIndex()];
-              if (new_index == -1)
-                throw Exception("Boundary " + ToString(sel.GetIndex()) + " with name " + mesh.GetBCName(sel.GetIndex() - 1) + " extruded, but no new material specified for it!");
-              el.SetIndex(new_mat_nrs[sel.GetIndex()]);
+              el.PNums().Range(0, np) = el.PNums().Range(np, 2 * np);
+              for (auto i : Range(np))
+                el[np + i] = newPoint(points[i], j, groups[i]);
               if (add_volume_element)
                 mesh.AddVolumeElement(el);
               else
                 {
                   // Let the volume mesher fill the hole with pyramids/tets
-                  // To insert pyramids, we need close surface identifications on open
-                  // quads
-                  for (auto i : Range(points))
-                    if (numGroups(sel[i]) == 1)
+                  // To insert pyramids, we need close surface identifications on open quads
+                  for (auto i : Range(np))
+                    if (numGroups(el[i]) == 1)
                       {
                         auto pi0 = el[i];
-                        auto pi1 = el[i + points.Size()];
+                        auto pi1 = el[np + i];
                         auto nr = identifications.Get(pi0, pi1);
                         if (nr == 0)
-                          identifications.Add(el[i], el[i + points.Size()], identnr);
+                          identifications.Add(pi0, pi1, identnr);
                       }
                 }
             }
           Element2d newel = sel;
-          for (auto i : Range(points))
-            newel[i] = newPoint(sel[i], -1, groups[i]);
-          newel.SetIndex(si_map[sel.GetIndex()]);
+          for (auto i : Range(np))
+            newel[i] = newPoint(points[i], -1, groups[i]);
+          if (surfacefacs[iface] > 0)
+            Swap(newel[0], newel[2]); // swap back
+          newel.SetIndex(si_map[iface]);
           new_sels.Append(newel);
         }
-      if (is_boundary_moved.Test(sel.GetIndex()))
+      if (is_boundary_moved.Test(iface))
         {
           auto& sel = mesh[si];
           for (auto& p : sel.PNums())
@@ -1065,7 +1065,7 @@ void BoundaryLayerTool ::InsertNewElements(
     }
 }
 
-void BoundaryLayerTool ::SetDomInOut()
+void BoundaryLayerTool ::SetDomInOut ()
 {
   if (insert_only_volume_elements)
     return;
@@ -1081,7 +1081,7 @@ void BoundaryLayerTool ::SetDomInOut()
       }
 }
 
-void BoundaryLayerTool ::SetDomInOutSides()
+void BoundaryLayerTool ::SetDomInOutSides ()
 {
   // Set the domin/domout entries for face descriptors on the "side" of new boundary layers
   if (insert_only_volume_elements)
@@ -1155,7 +1155,7 @@ void BoundaryLayerTool ::SetDomInOutSides()
     }
 }
 
-void BoundaryLayerTool ::AddSegments()
+void BoundaryLayerTool ::AddSegments ()
 {
   auto& new_segs =
     insert_only_volume_elements ? new_segments_on_moved_bnd : new_segments;
@@ -1199,14 +1199,14 @@ void BoundaryLayerTool ::AddSegments()
     mesh.AddSegment(seg);
 }
 
-void BoundaryLayerTool ::AddSurfaceElements()
+void BoundaryLayerTool ::AddSurfaceElements ()
 {
   for (auto& sel :
        insert_only_volume_elements ? new_sels_on_moved_bnd : new_sels)
     mesh.AddSurfaceElement(sel);
 }
 
-void BoundaryLayerTool ::ProcessParameters()
+void BoundaryLayerTool ::ProcessParameters ()
 {
   if (int* bc = get_if<int>(&params.boundary); bc)
     {
@@ -1374,7 +1374,7 @@ void BoundaryLayerTool ::ProcessParameters()
     domains.Invert();
 }
 
-void BoundaryLayerTool ::Perform()
+void BoundaryLayerTool ::Perform ()
 {
   if (domains.NumSet() == 0)
     return;

libsrc/meshing/boundarylayer.hpp

@@ -22,15 +22,15 @@ struct SpecialBoundaryPoint
     Vec<3> growth_vector;
     Array<PointIndex> new_points;
 
-    GrowthGroup(FlatArray<int> faces_, FlatArray<Vec<3>> normals);
-    GrowthGroup(const GrowthGroup&) = default;
-    GrowthGroup() = default;
+    GrowthGroup (FlatArray<int> faces_, FlatArray<Vec<3>> normals);
+    GrowthGroup (const GrowthGroup&) = default;
+    GrowthGroup () = default;
   };
   Array<GrowthGroup> growth_groups;
   Vec<3> separating_direction;
 
-  SpecialBoundaryPoint(const std::map<int, Vec<3>>& normals);
-  SpecialBoundaryPoint() = default;
+  SpecialBoundaryPoint (const std::map<int, Vec<3>>& normals);
+  SpecialBoundaryPoint () = default;
 };
 
 DLL_HEADER void GenerateBoundaryLayer (Mesh& mesh,
@@ -41,7 +41,7 @@ DLL_HEADER int /* new_domain_number */ GenerateBoundaryLayer2 (Mesh& mesh, int d
 class BoundaryLayerTool
 {
 public:
-  BoundaryLayerTool(Mesh& mesh_, const BoundaryLayerParameters& params_);
+  BoundaryLayerTool (Mesh& mesh_, const BoundaryLayerParameters& params_);
   void ProcessParameters ();
   void Perform ();
 

libsrc/meshing/boundarylayer2d.cpp

@@ -243,10 +243,6 @@ namespace netgen
 
      Array<SegmentIndex> segments;
 
-    // surface index map
-    Array<int> si_map(mesh.GetNFD()+2);
-    si_map = -1;
-
     // int fd_old = mesh.GetNFD();
 
     int max_edge_nr = -1;
@@ -254,8 +250,8 @@ namespace netgen
 
     for(const auto& seg : line_segments)
     {
-      if(seg.epgeominfo[0].edgenr > max_edge_nr)
-        max_edge_nr = seg.epgeominfo[0].edgenr;
+      if(seg.edgenr > max_edge_nr)
+        max_edge_nr = seg.edgenr;
       if(seg.domin > max_domain)
          max_domain = seg.domin;
       if(seg.domout > max_domain)
@@ -263,6 +259,7 @@ namespace netgen
     }
 
     int new_domain = max_domain+1;
+    int new_edge_nr = max_edge_nr+1;
 
     BitArray active_boundaries(max_edge_nr+1);
     BitArray active_segments(nseg);
@@ -288,7 +285,10 @@ namespace netgen
         // int new_fd_index =
         mesh.AddFaceDescriptor(new_fd);
         if(should_make_new_domain)
-           mesh.SetBCName(new_domain-1, "mapped_" + mesh.GetBCName(domain-1));
+          {
+           mesh.SetMaterial(new_domain, "layer_" + mesh.GetMaterial(domain));
+           mesh.SetBCName(new_edge_nr - 1, "moved");
+          }
     }
 
     for(auto segi : Range(line_segments))
@@ -618,8 +618,6 @@ namespace netgen
            }
         }
 
-     map<pair<PointIndex, PointIndex>, int> seg2edge;
-
      // insert new elements ( and move old ones )
      for(auto si : moved_segs)
      {
@@ -629,8 +627,6 @@ namespace netgen
         auto & pm0 = mapto[seg[0]];
         auto & pm1 = mapto[seg[1]];
 
-        // auto newindex = si_map[domain];
-
         Segment s = seg;
         s.geominfo[0] = {};
         s.geominfo[1] = {};
@@ -638,10 +634,10 @@ namespace netgen
         s[1] = pm1.Last();
         s[2] = PointIndex::INVALID;
         auto pair = s[0] < s[1] ? make_pair(s[0], s[1]) : make_pair(s[1], s[0]);
-        if(seg2edge.find(pair) == seg2edge.end())
-           seg2edge[pair] = ++max_edge_nr;
-        s.edgenr = seg2edge[pair];
-        s.si = seg.si;
+        s.edgenr = new_edge_nr;
+        s.epgeominfo[0].edgenr = -1;
+        s.epgeominfo[1].edgenr = -1;
+        s.si = s.edgenr;
         mesh.AddSegment(s);
 
         for ( auto i : Range(thicknesses))

libsrc/meshing/boundarylayer_interpolate.cpp

@@ -70,7 +70,7 @@ BuildNeighbors (FlatArray<PointIndex> points, const Mesh& mesh)
   return neighbors;
 }
 
-void BoundaryLayerTool ::InterpolateGrowthVectors()
+void BoundaryLayerTool ::InterpolateGrowthVectors ()
 {
   point_types.SetSize(mesh.GetNP());
   for (auto p : mesh.Points().Range())
@@ -139,7 +139,7 @@ void BoundaryLayerTool ::InterpolateGrowthVectors()
                   faces.Append(sei);
             }
 
-          if (faces.Size() == 2)
+          if (faces.Size() == 2 && mesh[faces[0]].GetIndex() != mesh[faces[1]].GetIndex())
             {
               auto n0 = getNormal(mesh[faces[0]]);
               auto n1 = getNormal(mesh[faces[1]]);
@@ -299,7 +299,7 @@ void BoundaryLayerTool ::InterpolateGrowthVectors()
     }
 }
 
-void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors()
+void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors ()
 {
   static Timer tall("InterpolateSurfaceGrowthVectors");
   RegionTimer rtall(tall);
@@ -399,7 +399,7 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors()
     growthvectors[pi] += corrections[pi];
 }
 
-void BoundaryLayerTool ::FixSurfaceElements()
+void BoundaryLayerTool ::FixSurfaceElements ()
 {
   static Timer tall("FixSurfaceElements");
   RegionTimer rtall(tall);

libsrc/meshing/boundarylayer_limiter.hpp

@@ -29,7 +29,7 @@ struct GrowthVectorLimiter
   Array<PointIndex, PointIndex> map_from;
   Table<SurfaceElementIndex, PointIndex> p2sel;
 
-  GrowthVectorLimiter(BoundaryLayerTool& tool_)
+  GrowthVectorLimiter (BoundaryLayerTool& tool_)
     : tool(tool_), params(tool_.params), mesh(tool_.mesh), height(tool_.total_height), growthvectors(tool_.growthvectors), map_from(mesh.Points().Size())
   {
     changed_domains = tool.domains;
@@ -287,7 +287,7 @@ struct GrowthVectorLimiter
     if (factor == 0.0)
       return;
     // for (PointIndex pi : IntRange(tool.np, mesh.GetNP()))
-    for (PointIndex pi : mesh.Points().Range().Modify(tool.np,0))
+    for (PointIndex pi : mesh.Points().Range().Modify(tool.np, 0))
       {
         // auto pi_from = map_from[pi];
         std::set<PointIndex> pis;
@@ -464,8 +464,7 @@ struct GrowthVectorLimiter
   }
 
   template <typename TFunc>
-  void FindTreeIntersections (double trig_shift, double seg_shift, TFunc f,
-                              TBitArray<PointIndex> * relevant_points = nullptr)
+  void FindTreeIntersections (double trig_shift, double seg_shift, TFunc f, TBitArray<PointIndex>* relevant_points = nullptr)
   {
     static Timer t("GrowthVectorLimiter::FindTreeIntersections");
     RegionTimer rt(t);
@@ -506,6 +505,25 @@ struct GrowthVectorLimiter
     RegionTimer reg(t);
     // check if surface trigs are intersecting each other
     bool changed = true;
+    std::set<PointIndex> special_points;
+
+    if (tool.insert_only_volume_elements)
+      for (auto [pi, special_point] : tool.special_boundary_points)
+        {
+          special_points.insert(pi);
+          for (auto& group : special_point.growth_groups)
+            special_points.insert(group.new_points.Last());
+        }
+
+    auto skip_trig = [&] (const Element2d& tri) {
+      if (!tool.insert_only_volume_elements)
+        return false;
+      for (auto pi : tri.PNums())
+        if (special_points.find(pi) != special_points.end())
+          return true;
+      return false;
+    };
+
     while (changed)
       {
         changed = false;
@@ -517,6 +535,9 @@ struct GrowthVectorLimiter
           {
             const Element2d& tri = Get(sei);
 
+            if (skip_trig(tri))
+              continue;
+
             Box<3> box(Box<3>::EMPTY_BOX);
             for (PointIndex pi : tri.PNums())
               box.Add(GetPoint(pi, 1.0, true));
@@ -529,6 +550,9 @@ struct GrowthVectorLimiter
           {
             const Element2d& tri = Get(sei);
 
+            if (skip_trig(tri))
+              continue;
+
             Box<3> box(Box<3>::EMPTY_BOX);
             for (PointIndex pi : tri.PNums())
               box.Add(GetPoint(pi, 1.0, true));
@@ -629,7 +653,7 @@ struct GrowthVectorLimiter
     size_t limit_counter = 1;
 
     TBitArray<PointIndex> relevant_points, relevant_points_next;
-    relevant_points.SetSize(mesh.Points().Size() + 1);   
+    relevant_points.SetSize(mesh.Points().Size() + 1);
     relevant_points_next.SetSize(mesh.Points().Size() + 1);
     relevant_points.Set();
 
@@ -685,8 +709,9 @@ struct GrowthVectorLimiter
     for (auto pi : Range(growthvectors))
       check_point(pi);
 
-    for (auto& [special_pi, special_point] : tool.special_boundary_points)
-      check_point(special_pi);
+    if (!tool.insert_only_volume_elements)
+      for (auto& [special_pi, special_point] : tool.special_boundary_points)
+        check_point(special_pi);
   }
 
   void Perform ()

libsrc/meshing/fieldlines.cpp

@@ -9,6 +9,25 @@
 
 namespace netgen
 {
+  inline int GetVolElement(const Mesh& mesh, const Point<3>& p,
+                           double* lami)
+  {
+    if(mesh.GetDimension() == 3)
+      {
+        auto ei = mesh.GetElementOfPoint(p, lami, true);
+        if(!ei.IsValid())
+          return -1;
+        return ei;
+      }
+    else
+      {
+        auto ei = mesh.GetSurfaceElementOfPoint(p, lami, true);
+        if(!ei.IsValid())
+          return -1;
+        return ei;
+      }
+  }
+
   RKStepper :: ~RKStepper() 
   {
     delete a;
@@ -190,9 +209,9 @@ namespace netgen
 	
 	for(int i=0; i<potential_startpoints.Size(); i++)
 	  {
-	    int elnr = mesh.GetElementOfPoint(potential_startpoints[i],lami,true) - 1;
-	    if(elnr == -1)
-	      continue;
+	    int elnr = GetVolElement(mesh, potential_startpoints[i], lami);
+            if (elnr == -1)
+              continue;
 
 	    mesh.SetPointSearchStartElement(elnr);
 	    
@@ -289,8 +308,7 @@ namespace netgen
     dirstart.SetSize(0);
     dirstart.Append(0);
 
-
-    int startelnr = mesh.GetElementOfPoint(startpoint,startlami,true) - 1;
+    int startelnr = GetVolElement(mesh, startpoint,startlami);
     (*testout) << "p = " << startpoint << "; elnr = " << startelnr << endl;
     if (startelnr == -1)
       return;
@@ -342,7 +360,7 @@ namespace netgen
             Point<3> newp;
 	    while(!stepper.GetNextData(newp,dummyt,h) && elnr != -1)
 	      {
-		elnr = mesh.GetElementOfPoint(newp,lami,true) - 1;
+		elnr = GetVolElement(mesh, newp, lami);
 		if(elnr != -1)
 		  {
 		    mesh.SetPointSearchStartElement(elnr);

libsrc/meshing/improve3.cpp

@@ -2194,7 +2194,7 @@ void MeshOptimize3d :: SwapImproveSurface (
 
 double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
                                         Table<ElementIndex, PointIndex> & elementsonnode,
-                                        DynamicTable<SurfaceElementIndex, PointIndex> & belementsonnode, bool check_only )
+                                        DynamicTable<SurfaceElementIndex, PointIndex> & belementsonnode, bool conform_segments, bool check_only )
 {
   PointIndex pi1, pi2, pi3, pi4, pi5;
   Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET);
@@ -2228,28 +2228,31 @@ double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
   }
 
 
-  bool bface = 0;
-  for (int k = 0; k < belementsonnode[pi1].Size(); k++)
-  {
-      const Element2d & bel =
-        mesh[belementsonnode[pi1][k]];
+  if(!conform_segments)
+    {
+      bool bface = 0;
+      for (int k = 0; k < belementsonnode[pi1].Size(); k++)
+      {
+          const Element2d & bel =
+            mesh[belementsonnode[pi1][k]];
 
-      bool bface1 = 1;
-      for (int l = 0; l < 3; l++)
-          if (bel[l] != pi1 && bel[l] != pi2 && bel[l] != pi3)
+          bool bface1 = 1;
+          for (int l = 0; l < 3; l++)
+              if (bel[l] != pi1 && bel[l] != pi2 && bel[l] != pi3)
+              {
+                  bface1 = 0;
+                  break;
+              }
+
+          if (bface1)
           {
-              bface1 = 0;
+              bface = 1;
               break;
           }
-
-      if (bface1)
-      {
-          bface = 1;
-          break;
       }
-  }
 
-  if (bface) return 0.0;
+      if (bface) return 0.0;
+    }
 
 
   FlatArray<ElementIndex> row = elementsonnode[pi1];
@@ -2367,11 +2370,11 @@ double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
   2 -> 3 conversion
 */
 
-void MeshOptimize3d :: SwapImprove2 ()
+void MeshOptimize3d :: SwapImprove2 (bool conform_segments)
 {
   static Timer t("MeshOptimize3d::SwapImprove2"); RegionTimer reg(t);
 
-  if (goal == OPT_CONFORM) return;
+  if (!conform_segments && goal == OPT_CONFORM) return;
 
   mesh.BuildBoundaryEdges(false);
 
@@ -2433,7 +2436,7 @@ void MeshOptimize3d :: SwapImprove2 ()
 
             for (int j = 0; j < 4; j++)
               {
-                double d_badness = SwapImprove2( eli1, j, elementsonnode, belementsonnode, true);
+                double d_badness = SwapImprove2( eli1, j, elementsonnode, belementsonnode, conform_segments, true);
                 if(d_badness<0.0)
                     my_faces_with_improvement.Append( std::make_tuple(d_badness, eli1, j) );
               }
@@ -2449,7 +2452,7 @@ void MeshOptimize3d :: SwapImprove2 ()
     {
       if(mesh[eli].IsDeleted())
           continue;
-      if(SwapImprove2( eli, j, elementsonnode, belementsonnode, false) < 0.0)
+      if(SwapImprove2( eli, j, elementsonnode, belementsonnode, conform_segments, false) < 0.0)
           cnt++;
     }
 

libsrc/meshing/improve3.hpp

@@ -45,8 +45,8 @@ public:
   void SwapImprove (const TBitArray<ElementIndex> * working_elements = NULL);
   void SwapImproveSurface (const TBitArray<ElementIndex> * working_elements = NULL,
 			   const NgArray< idmap_type* > * idmaps = NULL);
-  void SwapImprove2 ();
-  double SwapImprove2 (ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, DynamicTable<SurfaceElementIndex, PointIndex> & belementsonnode, bool check_only=false );
+  void SwapImprove2 (bool conform_segments = false);
+  double SwapImprove2 (ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, DynamicTable<SurfaceElementIndex, PointIndex> & belementsonnode, bool conform_segments, bool check_only=false );
 
   void ImproveMesh() { mesh.ImproveMesh(mp, goal); }
 

libsrc/meshing/localh.cpp

@@ -196,54 +196,41 @@ namespace netgen
             fabs (p(1) - root->xmid[1]) > root->h2)
           return;
         
-        if (GetH(p) <= 1.2 * h) return;
-        
-        GradingBox * box = root;
-        GradingBox * nbox = root;
-        GradingBox * ngb;
-        int childnr;
-        double x1[3], x2[3];
-        
-        while (nbox)
-          {
-            box = nbox;
-            childnr = 0;
-            if (p(0) > box->xmid[0]) childnr += 1;
-            if (p(1) > box->xmid[1]) childnr += 2;
-            nbox = box->childs[childnr];
-          };
+        GradingBox * box = Find(p);
+        if (box->HOpt() <= 1.2 * h) return;
         
         while (2 * box->h2 > h)
           {
-            childnr = 0;
-            if (p(0) > box->xmid[0]) childnr += 1;
-            if (p(1) > box->xmid[1]) childnr += 2;
-            
+            int childnr = 0;
+            double x1[3], x2[3];
+        
             double h2 = box->h2;
-            if (childnr & 1)
+            if (p(0) > box->xmid[0])
               {
+                childnr += 1;
                 x1[0] = box->xmid[0];
-                x2[0] = x1[0]+h2;   // box->x2[0];
+                x2[0] = x1[0]+h2;
               }
             else
               {
                 x2[0] = box->xmid[0];
-                x1[0] = x2[0]-h2;   // box->x1[0];
+                x1[0] = x2[0]-h2;
               }
             
-            if (childnr & 2)
+            if (p(1) > box->xmid[1])
               {
+                childnr += 2;
                 x1[1] = box->xmid[1];
-                x2[1] = x1[1]+h2;   // box->x2[1];
+                x2[1] = x1[1]+h2;
               }
             else
               {
                 x2[1] = box->xmid[1];
-                x1[1] = x2[1]-h2;   // box->x1[1];
+                x1[1] = x2[1]-h2;
               }
             x1[2] = x2[2] = 0;
 
-            ngb = new GradingBox (x1, x2);
+            auto ngb = new GradingBox (x1, x2);
             box->childs[childnr] = ngb;
             ngb->father = box;
             
@@ -276,67 +263,52 @@ namespace netgen
             fabs (p(2) - root->xmid[2]) > root->h2)
           return;
         
-        if (GetH(p) <= 1.2 * h) return;
-        
-        GradingBox * box = root;
-        GradingBox * nbox = root;
-        GradingBox * ngb;
-        int childnr;
-        double x1[3], x2[3];
-        
-        while (nbox)
-          {
-            box = nbox;
-            childnr = 0;
-            if (p(0) > box->xmid[0]) childnr += 1;
-            if (p(1) > box->xmid[1]) childnr += 2;
-            if (p(2) > box->xmid[2]) childnr += 4;
-            nbox = box->childs[childnr];
-          };
-        
+        GradingBox * box = Find(p);
+        if (box->HOpt() <= 1.2 * h) return;
         
         while (2 * box->h2 > h)
           {
-            childnr = 0;
-            if (p(0) > box->xmid[0]) childnr += 1;
-            if (p(1) > box->xmid[1]) childnr += 2;
-            if (p(2) > box->xmid[2]) childnr += 4;
-            
+            double x1[3], x2[3];
+            int childnr = 0;
             double h2 = box->h2;
-            if (childnr & 1)
+
+            if (p(0) > box->xmid[0])
               {
+                childnr += 1;
                 x1[0] = box->xmid[0];
-                x2[0] = x1[0]+h2;   // box->x2[0];
+                x2[0] = x1[0]+h2;
               }
             else
               {
                 x2[0] = box->xmid[0];
-                x1[0] = x2[0]-h2;   // box->x1[0];
+                x1[0] = x2[0]-h2;
               }
             
-            if (childnr & 2)
+            if (p(1) > box->xmid[1])
               {
+                childnr += 2;
                 x1[1] = box->xmid[1];
-                x2[1] = x1[1]+h2;   // box->x2[1];
+                x2[1] = x1[1]+h2;
               }
             else
               {
                 x2[1] = box->xmid[1];
-                x1[1] = x2[1]-h2;   // box->x1[1];
+                x1[1] = x2[1]-h2;
               }
             
-            if (childnr & 4)
+            if (p(2) > box->xmid[2])
               {
+                childnr += 4;
                 x1[2] = box->xmid[2];
-                x2[2] = x1[2]+h2;  // box->x2[2];
+                x2[2] = x1[2]+h2;
               }
             else
               {
                 x2[2] = box->xmid[2];
-                x1[2] = x2[2]-h2;  // box->x1[2];
+                x1[2] = x2[2]-h2;
               }
             
-            ngb = new GradingBox (x1, x2);
+            auto ngb = new GradingBox (x1, x2);
             box->childs[childnr] = ngb;
             ngb->father = box;
             
@@ -366,37 +338,7 @@ namespace netgen
 
   double LocalH :: GetH (Point<3> x) const
   {
-    const GradingBox * box = root;
-    if (dimension == 2)
-      {
-        while (1)
-          {
-            int childnr = 0;
-            if (x(0) > box->xmid[0]) childnr += 1;
-            if (x(1) > box->xmid[1]) childnr += 2;
-            
-            if (box->childs[childnr])
-              box = box->childs[childnr];
-            else
-              return box->hopt;
-          }
-      }
-    else
-      {
-        while (1)
-          {
-            int childnr = 0;
-            if (x(0) > box->xmid[0]) childnr += 1;
-            if (x(1) > box->xmid[1]) childnr += 2;
-            if (x(2) > box->xmid[2]) childnr += 4;
-            
-            if (box->childs[childnr])
-              box = box->childs[childnr];
-            else
-              return box->hopt;
-          }
-      }
-      
+    return Find(x)->HOpt();
   }
 
 
@@ -488,6 +430,41 @@ namespace netgen
   }
 
 
+  GradingBox * LocalH :: Find (Point<3> p) const
+  {
+    GradingBox * box = root;
+    if (dimension == 2)
+      {
+        while (1)
+          {
+            int childnr = 0;
+            if (p(0) > box->xmid[0]) childnr += 1;
+            if (p(1) > box->xmid[1]) childnr += 2;
+
+            if (box->childs[childnr])
+              box = box->childs[childnr];
+            else
+              return box;
+          }
+      }
+    else
+      {
+        while (1)
+          {
+            int childnr = 0;
+            if (p(0) > box->xmid[0]) childnr += 1;
+            if (p(1) > box->xmid[1]) childnr += 2;
+            if (p(2) > box->xmid[2]) childnr += 4;
+
+            if (box->childs[childnr])
+              box = box->childs[childnr];
+            else
+              return box;
+          }
+      }
+    return nullptr;
+  }
+
 
   void LocalH :: FindInnerBoxes (const AdFront3 & adfront,
 				 int (*testinner)(const Point3d & p1))

libsrc/meshing/localh.hpp

@@ -54,6 +54,7 @@ namespace netgen
 
     Point<3> PMid() const { return Point<3> (xmid[0], xmid[1], xmid[2]); }
     double H2() const { return h2; }
+    double HOpt() const { return hopt; }
 
     bool HasChilds() const
     {
@@ -119,6 +120,8 @@ namespace netgen
 
     void CutBoundary (const Box<3> & box)
     { CutBoundaryRec (box.PMin(), box.PMax(), root); }
+
+    GradingBox * Find(Point<3> p) const;
   
     /// find inner boxes
     void FindInnerBoxes (const class AdFront3 & adfront,

libsrc/meshing/meshclass.cpp

@@ -12,15 +12,15 @@
 
 namespace netgen
 {
-  int Find3dElement (const Mesh& mesh,
-                     const netgen::Point<3> & p,
-		     double * lami,
-		     const NgArray<int> * const indices,
-		     BoxTree<3> * searchtree,
-		     const bool allowindex = true)
+  ElementIndex Find3dElement (const Mesh& mesh,
+                              const netgen::Point<3> & p,
+                              double * lami,
+                              optional<FlatArray<int>> indices,
+                              BoxTree<3, ElementIndex> * searchtree,
+                              const bool allowindex)
   {
     int ne = 0;
-    NgArray<int> locels;
+    Array<ElementIndex> locels;
     if (searchtree)
       {
         searchtree->GetIntersecting (p, p, locels);
@@ -29,92 +29,90 @@ namespace netgen
     else
       ne = mesh.GetNE();
 
-    for (int i = 1; i <= ne; i++)
+    for (auto i : Range(ne))
       {
-        int ii;
+        ElementIndex ei;
 
         if (searchtree)
-          ii = locels.Get(i);
+          ei = locels[i];
         else
-          ii = i;
+          ei = i;
 
-        if(indices != NULL && indices->Size() > 0)
+        if(indices && indices->Size() > 0)
           {
-            bool contained = indices->Contains(mesh.VolumeElement(ii).GetIndex());
+            bool contained = indices->Contains(mesh[ei].GetIndex());
             if((allowindex && !contained) || (!allowindex && contained)) continue;
           }
 
-        if(mesh.PointContainedIn3DElement(p,lami,ii))
-            return ii;
+        if(mesh.PointContainedIn3DElement(p,lami,ei))
+          return ei;
       }
 
     // Not found, try uncurved variant:
-    for (int i = 1; i <= ne; i++)
+    for (auto i : Range(ne))
       {
-        int ii;
+        ElementIndex ei;
 
         if (searchtree)
-          ii = locels.Get(i);
+          ei = locels[i];
         else
-          ii = i;
+          ei = i;
 
-        if(indices != NULL && indices->Size() > 0)
+        if(indices && indices->Size() > 0)
           {
-            bool contained = indices->Contains(mesh.VolumeElement(ii).GetIndex());
+            bool contained = indices->Contains(mesh[ei].GetIndex());
             if((allowindex && !contained) || (!allowindex && contained)) continue;
           }
 
 
-        if(mesh.PointContainedIn3DElementOld(p,lami,ii))
+        if(mesh.PointContainedIn3DElementOld(p,lami,ei))
           {
             (*testout) << "WARNING: found element of point " << p <<" only for uncurved mesh" << endl;
-            return ii;
+            return ei;
           }
       }
-    return 0;
+    return ElementIndex::INVALID;
   }
 
-  int Find2dElement (const Mesh& mesh,
-                     const netgen::Point<3> & p,
-		     double * lami,
-		     const NgArray<int> * const indices,
-		     BoxTree<3> * searchtree,
-		     const bool allowindex = true)
+  SurfaceElementIndex
+  Find2dElement (const Mesh& mesh,
+                 const netgen::Point<3> & p,
+                 double * lami,
+                 std::optional<FlatArray<int>> indices,
+                 BoxTree<3, SurfaceElementIndex> * searchtree,
+                 bool allowindex)
   {
     double vlam[3];
-    int velement = 0;
+    ElementIndex velement = ElementIndex::INVALID;
 
     if(mesh.GetNE())
-      velement = Find3dElement(mesh, p,vlam,NULL,searchtree,allowindex);
+      {
+        if(searchtree)
+          const_cast<Mesh&>(mesh).BuildElementSearchTree(3);
+        velement = Find3dElement(mesh, p,vlam, nullopt,searchtree ? mesh.GetElementSearchTree() : nullptr,allowindex);
+      }
 
     //(*testout) << "p " << p << endl;
     //(*testout) << "velement " << velement << endl;
 
     // first try to find a volume element containing p and project to face
-    if(velement!=0)
+    if(velement.IsValid())
     {
       auto & topology = mesh.GetTopology();
-      // NgArray<int> faces;
-      // topology.GetElementFaces(velement,faces);
-      auto faces = Array<int> (topology.GetFaces(ElementIndex(velement-1)));
-
-      //(*testout) << "faces " << faces << endl;
-
-      for(int i=0; i<faces.Size(); i++)
-        faces[i] = topology.GetFace2SurfaceElement(faces[i])+1;
-
-      //(*testout) << "surfel " << faces << endl;
+      const auto & fnrs = topology.GetFaces(velement);
+      auto faces = ArrayMem<SurfaceElementIndex,4>();
+      for(auto face : fnrs)
+        faces.Append(topology.GetFace2SurfaceElement(face));
 
       for(int i=0; i<faces.Size(); i++)
         {
-          if(faces[i] == 0)
+          if(!faces[i].IsValid())
             continue;
           auto sel = mesh.SurfaceElement(faces[i]);
-          if(indices && indices->Size() != 0 && !indices->Contains(sel.GetIndex()))
+          if(indices && indices->Size() > 0 && !indices->Contains(sel.GetIndex()))
             continue;
 
-          auto & el = mesh.VolumeElement(velement);
-
+          auto & el = mesh[velement];
           if (el.GetType() == TET)
           {
             double lam4[4] = { vlam[0], vlam[1], vlam[2], 1.0-vlam[0]-vlam[1]-vlam[2] };
@@ -127,7 +125,7 @@ namespace netgen
                 for(auto k : Range(4))
                   if(sel[j] == el[k])
                     lami[j-1] = lam4[k]/(1.0-face_lam);
-              return faces[i];
+              return SurfaceElementIndex(faces[i]);
             }
           }
 
@@ -139,9 +137,8 @@ namespace netgen
     // Did't find any matching face of a volume element, search 2d elements directly
     int ne;
 
-    NgArray<int> locels;
-    // TODO: build search tree for surface elements
-    if (!mesh.GetNE() && searchtree)
+    Array<SurfaceElementIndex> locels;
+    if (searchtree)
       {
         searchtree->GetIntersecting (p, p, locels);
         ne = locels.Size();
@@ -149,37 +146,38 @@ namespace netgen
     else
       ne = mesh.GetNSE();
 
-    for (int i = 1; i <= ne; i++)
+    for (auto i : Range(ne))
       {
-        int ii;
+        SurfaceElementIndex ii;
 
         if (locels.Size())
-          ii = locels.Get(i);
+          ii = locels[i];
         else
           ii = i;
 
-        if(indices != NULL && indices->Size() > 0)
+        if(indices && indices->Size() > 0)
           {
-            bool contained = indices->Contains(mesh.SurfaceElement(ii).GetIndex());
+            bool contained = indices->Contains(mesh[ii].GetIndex());
             if((allowindex && !contained) || (!allowindex && contained)) continue;
           }
-
-        if(mesh.PointContainedIn2DElement(p,lami,ii)) return ii;
-
+        if(mesh.PointContainedIn2DElement(p,lami,ii))
+          return ii;
       }
-    return 0;
+    return SurfaceElementIndex::INVALID;
   }
 
-  int Find1dElement (const Mesh& mesh,
-                     const netgen::Point<3> & p,
-		     double * lami,
-		     const NgArray<int> * const indices,
-		     BoxTree<3> * searchtree,
-		     const bool allowindex = true)
+  SegmentIndex Find1dElement (const Mesh& mesh,
+                              const netgen::Point<3> & p,
+                              double * lami,
+                              std::optional<FlatArray<int>> indices,
+                              BoxTree<3> * searchtree,
+                              const bool allowindex = true)
   {
     double vlam[3];
-    int velement = Find2dElement(mesh, p, vlam, NULL, searchtree, allowindex);
-    if(velement == 0)
+    if(searchtree)
+      const_cast<Mesh&>(mesh).BuildElementSearchTree(2);
+    auto velement = Find2dElement(mesh, p, vlam, nullopt, searchtree ? mesh.GetSurfaceElementSearchTree() : nullptr, allowindex);
+    if(!velement.IsValid())
       return 0;
 
     vlam[2] = 1.-vlam[0] - vlam[1];
@@ -235,8 +233,8 @@ namespace netgen
     : topology(*this), surfarea(*this)
   {
     lochfunc = {nullptr};
-    elementsearchtree = nullptr;
-    elementsearchtreets = NextTimeStamp();
+    for(auto i : Range(4))
+      elementsearchtreets[i] = NextTimeStamp();
     majortimestamp = timestamp = NextTimeStamp();
     hglob = 1e10;
     hmin = 0;
@@ -4821,6 +4819,21 @@ namespace netgen
         for (auto & seg : LineSegments())
           seg.si = seg.edgenr;
       }
+    if (dimension == 3 && dim == 1)
+      {
+        for(auto str : materials)
+          delete str;
+        materials.SetSize(0);
+        for(auto str : bcnames)
+          delete str;
+        bcnames.SetSize(0);
+        for(auto str: cd2names)
+          materials.Append(str);
+        cd2names.SetSize(0);
+        for(auto str : cd3names)
+          bcnames.Append(str);
+        cd3names.SetSize(0);
+      }
     dimension = dim;
   }
 
@@ -5251,122 +5264,95 @@ namespace netgen
     RebuildSurfaceElementLists();
   }
 
-  void Mesh :: BuildElementSearchTree ()
+  void Mesh :: BuildElementSearchTree (int dim)
   {
-    if (elementsearchtreets == GetTimeStamp()) return;
+    if(dim < 2)
+      return;
+    if (elementsearchtreets[dim] == GetTimeStamp())
+      return;
 
     {
       std::lock_guard<std::mutex> guard(buildsearchtree_mutex);
-      if (elementsearchtreets != GetTimeStamp())
+      // check again to see if some other thread built while waiting for lock
+      if (elementsearchtreets[dim] == GetTimeStamp()) return;
+
+      elementsearchtreets[dim] = GetTimeStamp();
+
+      PrintMessage (4, "Rebuild element searchtree dim " + ToString(dim));
+          
+
+      Point3d pmin, pmax;
+      GetBox(pmin, pmax);
+      Box<3> box(pmin, pmax);
+      if (dim == 3)
+        elementsearchtree_vol = make_unique<BoxTree<3, ElementIndex>>(box);
+      else
+        elementsearchtree_surf = make_unique<BoxTree<3, SurfaceElementIndex>>(box);
+
+      if (dim == 3)
         {
-          NgLock lock(mutex);
-          lock.Lock();
-          
-          PrintMessage (4, "Rebuild element searchtree");
-          
-          elementsearchtree = nullptr;
-          
-          int ne = (dimension == 2) ? GetNSE() : GetNE();
-          if (dimension == 3 && !GetNE() && GetNSE())
-            ne = GetNSE();
-
-          if (ne) 
+          for(auto ei : volelements.Range())
             {
-              if (dimension == 2 || (dimension == 3 && !GetNE()) )
+              const auto& el = volelements[ei];
+              Box<3> box (Box<3>::EMPTY_BOX);
+              for (auto pi : el.PNums())
+                box.Add (points[pi]);
+
+              if(el.IsCurved() && curvedelems->IsElementCurved(ei))
                 {
-                  Box<3> box (Box<3>::EMPTY_BOX);
-                  for (SurfaceElementIndex sei = 0; sei < ne; sei++)
-                    // box.Add (points[surfelements[sei].PNums()]);
-                    for (auto pi : surfelements[sei].PNums())
-                      box.Add (points[pi]);
-                  
-                  box.Increase (1.01 * box.Diam());
-                  elementsearchtree = make_unique<BoxTree<3>> (box);
-                  
-                  for (SurfaceElementIndex sei = 0; sei < ne; sei++)
-                    {
-                      //  box.Set (points[surfelements[sei].PNums()]);
+                  // add edge/face midpoints to box
+                  auto eltype = el.GetType();
+                  const auto verts = topology.GetVertices(eltype);
 
-                      Box<3> box (Box<3>::EMPTY_BOX);
-                      for (auto pi : surfelements[sei].PNums())
-                        box.Add (points[pi]);
+                  const auto edges = FlatArray<const ELEMENT_EDGE>(topology.GetNEdges(eltype), topology.GetEdges0(eltype));
+                  for (const auto & edge: edges) {
+                    netgen::Point<3> lam = netgen::Point<3>(0.5* (verts[edge[0]] + verts[edge[1]]));
+                    auto p = netgen::Point<3>(0.0);
+                    curvedelems->CalcElementTransformation(lam,ei,p);
+                    box.Add(p);
+                  }
 
-                      auto & el = surfelements[sei];
-                      if(el.IsCurved() && curvedelems->IsSurfaceElementCurved(sei))
-                        {
-                          netgen::Point<2>  lami [4] = {netgen::Point<2>(0.5,0), netgen::Point<2>(0,0.5), netgen::Point<2>(0.5,0.5), netgen::Point<2>(1./3,1./3)};
-                          for (auto lam : lami)
-                            {
-                              netgen::Point<3> x;
-                              Mat<3,2> Jac;
-                          
-                              curvedelems->CalcSurfaceTransformation(lam,sei,x,Jac);
-                              box.Add (x);
-                            }
-                          box.Scale(1.2);
-                        }
-                      elementsearchtree -> Insert (box, sei+1);
+                  const auto faces = FlatArray<const ELEMENT_FACE>(topology.GetNFaces(eltype), topology.GetFaces0(eltype));
+                  for (const auto & face: faces) {
+                    netgen::Vec<3> lam = netgen::Vec<3>(verts[face[0]] + verts[face[1]] + verts[face[2]]);
+                    if(face[3] != -1) {
+                      lam += netgen::Vec<3>(verts[face[3]]);
+                      lam *= 0.25;
                     }
+                    else
+                      lam *= 1.0/3;
+                    auto p = netgen::Point<3>(0.0);
+                    curvedelems->CalcElementTransformation(netgen::Point<3>(lam),ei,p);
+                    box.Add(p);
+                  }
                 }
-              else
+              box.Scale(1.2);
+              elementsearchtree_vol -> Insert (box, ei);
+            }
+        }
+      else if (dim == 2)
+        {
+          for (auto ei : Range(surfelements))
+            {
+              const auto& el = surfelements[ei];
+              Box<3> box (Box<3>::EMPTY_BOX);
+              for (auto pi : el.PNums())
+                box.Add (points[pi]);
+
+              if(el.IsCurved() && curvedelems->IsSurfaceElementCurved(ei))
                 {
-                  Box<3> box (Box<3>::EMPTY_BOX);
-                  for (ElementIndex ei = 0; ei < ne; ei++)
-                    // box.Add (points[volelements[ei].PNums()]);
-                    for (auto pi : volelements[ei].PNums())
-                      box.Add (points[pi]);
-                  
-                  box.Increase (1.01 * box.Diam());
-                  elementsearchtree = make_unique<BoxTree<3>> (box);
-                  
-                  for (ElementIndex ei = 0; ei < ne; ei++)
+                  netgen::Point<2>  lami [4] = {netgen::Point<2>(0.5,0), netgen::Point<2>(0,0.5), netgen::Point<2>(0.5,0.5), netgen::Point<2>(1./3,1./3)};
+                  for (auto lam : lami)
                     {
-                      // box.Set (points[volelements[ei].PNums()]);
+                      netgen::Point<3> x;
+                      Mat<3,2> Jac;
 
-                      Box<3> box (Box<3>::EMPTY_BOX);
-                      for (auto pi : volelements[ei].PNums())
-                        box.Add (points[pi]);
-
-                      auto & el = volelements[ei];
-                      if(el.IsCurved() && curvedelems->IsElementCurved(ei))
-                      {
-                        // add edge/face midpoints to box
-                        auto eltype = el.GetType();
-                        const auto verts = topology.GetVertices(eltype);
-
-
-                        const auto edges = FlatArray<const ELEMENT_EDGE>(topology.GetNEdges(eltype), topology.GetEdges0(eltype));
-                        for (const auto & edge: edges) {
-                          netgen::Point<3> lam = netgen::Point<3>(0.5* (verts[edge[0]] + verts[edge[1]]));
-                          auto p = netgen::Point<3>(0.0);
-                          curvedelems->CalcElementTransformation(lam,ei,p);
-                          box.Add(p);
-                        }
-
-                        const auto faces = FlatArray<const ELEMENT_FACE>(topology.GetNFaces(eltype), topology.GetFaces0(eltype));
-                        for (const auto & face: faces) {
-                          netgen::Vec<3> lam = netgen::Vec<3>(verts[face[0]] + verts[face[1]] + verts[face[2]]);
-                          if(face[3] != -1) {
-                            lam += netgen::Vec<3>(verts[face[3]]);
-                            lam *= 0.25;
-                          }
-                          else
-                            lam *= 1.0/3;
-                          auto p = netgen::Point<3>(0.0);
-                          curvedelems->CalcElementTransformation(netgen::Point<3>(lam),ei,p);
-                          box.Add(p);
-                        }
-
-
-                        box.Scale(1.2);
-                      }
-
-
-                      elementsearchtree -> Insert (box, ei+1);
+                      curvedelems->CalcSurfaceTransformation(lam,ei,x,Jac);
+                      box.Add (x);
                     }
+                  box.Scale(1.2);
                 }
-              
-              elementsearchtreets = GetTimeStamp();
+              elementsearchtree_surf -> Insert (box, ei);
             }
         }
     }
@@ -5406,7 +5392,7 @@ namespace netgen
 
   bool Mesh :: PointContainedIn2DElement(const Point3d & p,
                                          double lami[3],
-                                         const int element,
+                                         SurfaceElementIndex ei,
                                          bool consider3D) const
   {
     Vec3d col1, col2, col3;
@@ -5417,9 +5403,9 @@ namespace netgen
 
     
     //SZ 
-    if(SurfaceElement(element).GetType()==QUAD)
+    if(surfelements[ei].GetType()==QUAD)
       {
-        const Element2d & el = SurfaceElement(element); 
+        const Element2d & el = surfelements[ei];
 
         const Point3d & p1 = Point(el.PNum(1)); 
         const Point3d & p2 = Point(el.PNum(2));
@@ -5438,7 +5424,7 @@ namespace netgen
           int i = 0;
           while(delta > 1e-16 && i < maxits)
             {
-              curvedelems->CalcSurfaceTransformation(lam,element-1,x,Jac);
+              curvedelems->CalcSurfaceTransformation(lam,ei,x,Jac);
               rhs = p - x;
               Jac.Solve(rhs,deltalam);
               lam += deltalam;
@@ -5767,7 +5753,7 @@ namespace netgen
       {
         //	  SurfaceElement(element).GetTets (loctets);
         loctrigs.SetSize(1);
-        loctrigs.Elem(1) = SurfaceElement(element);
+        loctrigs.Elem(1) = surfelements[ei];
 
 
 
@@ -5802,11 +5788,11 @@ namespace netgen
             //(*testout) << "col1 " << col1 << " col2 " << col2 << " col3 " << col3 << " rhs " << rhs << endl;
             //(*testout) << "sol " << sol << endl;
 
-            if (SurfaceElement(element).GetType() ==TRIG6 || curvedelems->IsSurfaceElementCurved(element-1))
+            if (surfelements[ei].GetType() ==TRIG6 || curvedelems->IsSurfaceElementCurved(ei))
               {
                 // netgen::Point<2> lam(1./3,1./3);
                 netgen::Point<2> lam(sol.X(), sol.Y());
-                if(SurfaceElement(element).GetType() != TRIG6)
+                if(surfelements[ei].GetType() != TRIG6)
                   {
                     lam[0] = 1-sol.X()-sol.Y();
                     lam[1] = sol.X();
@@ -5825,7 +5811,7 @@ namespace netgen
                 const int maxits = 30;
                 while(delta > 1e-16 && i<maxits)
                   {
-                    curvedelems->CalcSurfaceTransformation(lam,element-1,x,Jac);
+                    curvedelems->CalcSurfaceTransformation(lam,ei,x,Jac);
                     rhs = p-x;
                     Jac.Solve(rhs,deltalam);
                     
@@ -5844,7 +5830,7 @@ namespace netgen
                 sol.X() = lam(0);
                 sol.Y() = lam(1);
 
-                if (SurfaceElement(element).GetType() !=TRIG6 )
+                if (surfelements[ei].GetType() !=TRIG6 )
                   {
                     sol.Z() = sol.X();
                     sol.X() = sol.Y();
@@ -5876,7 +5862,7 @@ namespace netgen
 
   bool Mesh :: PointContainedIn3DElement(const Point3d & p,
                                          double lami[3],
-                                         const int element) const
+                                         ElementIndex ei) const
   {
     //bool oldresult = PointContainedIn3DElementOld(p,lami,element);
     //(*testout) << "old result: " << oldresult
@@ -5887,7 +5873,7 @@ namespace netgen
 
 
     const double eps = 1.e-4;
-    const Element & el = VolumeElement(element);
+    const Element & el = volelements[ei];
 
     netgen::Point<3> lam = 0.0;
 
@@ -5922,7 +5908,7 @@ namespace netgen
     const int maxits = 30;
     while(delta > 1e-16 && i<maxits)
       {
-        curvedelems->CalcElementTransformation(lam,element-1,x,Jac);
+        curvedelems->CalcElementTransformation(lam,ei,x,Jac);
         rhs = p-x;
         Jac.Solve(rhs,deltalam);
 
@@ -6044,86 +6030,72 @@ namespace netgen
   }
 
 
-  int Mesh :: GetElementOfPoint (const netgen::Point<3> & p,
-                                 double lami[3],
-                                 bool build_searchtree,
-                                 const int index,
-                                 const bool allowindex) const
+  ElementIndex Mesh :: GetElementOfPoint (const netgen::Point<3> & p,
+                                          double* lami,
+                                          bool build_searchtree,
+                                          int index,
+                                          bool allowindex) const
   {
     if(index != -1) 
       {
-        NgArray<int> dummy(1);
+        Array<int> dummy(1);
         dummy[0] = index;
-        return GetElementOfPoint(p,lami,&dummy,build_searchtree,allowindex);
+        return GetElementOfPoint(p,lami,dummy,build_searchtree,allowindex);
       }
     else
-      return GetElementOfPoint(p,lami,NULL,build_searchtree,allowindex);
+      return GetElementOfPoint(p,lami,nullopt,build_searchtree,allowindex);
   }
 
 
 
 
-  int Mesh :: GetElementOfPoint (const netgen::Point<3> & p,
-                                 double lami[3],
-                                 const NgArray<int> * const indices,
-                                 bool build_searchtree,
-                                 const bool allowindex) const
+  ElementIndex Mesh :: GetElementOfPoint (const netgen::Point<3> & p,
+                                          double* lami,
+                                          std::optional<FlatArray<int>> indices,
+                                          bool build_searchtree,
+                                          bool allowindex) const
   {
-    if ( (dimension == 2 && !GetNSE()) ||
-    	 (dimension == 3 && !GetNE() && !GetNSE()) )
-      return -1;
-
     if (build_searchtree)
-      const_cast<Mesh&>(*this).BuildElementSearchTree ();
-
-    if (dimension == 2 || (dimension==3 && !GetNE() && GetNSE()))
-      return Find2dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex);
-
-    return Find3dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex);
+      const_cast<Mesh&>(*this).BuildElementSearchTree (3);
+    return Find3dElement(*this, p, lami, indices, elementsearchtree_vol.get(), allowindex);
   }
 
 
 
-  int Mesh :: GetSurfaceElementOfPoint (const netgen::Point<3> & p,
-                                        double lami[3],
-                                        bool build_searchtree,
-                                        const int index,
-                                        const bool allowindex) const
+  SurfaceElementIndex Mesh ::
+  GetSurfaceElementOfPoint (const netgen::Point<3> & p,
+                            double* lami,
+                            bool build_searchtree,
+                            int index,
+                            bool allowindex) const
   {
-    if(index != -1) 
+    if(index != -1)
       {
-        NgArray<int> dummy(1);
+        Array<int> dummy(1);
         dummy[0] = index;
-        return GetSurfaceElementOfPoint(p,lami,&dummy,build_searchtree,allowindex);
+        return GetSurfaceElementOfPoint(p,lami,dummy,build_searchtree,allowindex);
       }
     else
-      return GetSurfaceElementOfPoint(p,lami,NULL,build_searchtree,allowindex);
+      return GetSurfaceElementOfPoint(p,lami,nullopt,build_searchtree,allowindex);
   }
 
-
-
-
-  int Mesh :: GetSurfaceElementOfPoint (const netgen::Point<3> & p,
-                                        double lami[3],
-                                        const NgArray<int> * const indices,
-                                        bool build_searchtree,
-                                        const bool allowindex) const
+  SurfaceElementIndex Mesh ::
+  GetSurfaceElementOfPoint (const netgen::Point<3> & p,
+                            double* lami,
+                            std::optional<FlatArray<int>> indices,
+                            bool build_searchtree,
+                            bool allowindex) const
   {
-    if (!GetNE() && build_searchtree)
-      const_cast<Mesh&>(*this).BuildElementSearchTree ();
-
-    if (dimension == 2)
-        return Find1dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex);
-    else
-        return Find2dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex);
-    return 0;
+    if (build_searchtree)
+      const_cast<Mesh&>(*this).BuildElementSearchTree(2);
+    return Find2dElement(*this, p, lami, indices, elementsearchtree_surf.get(), allowindex);
   }
 
 
   void Mesh::GetIntersectingVolEls(const Point3d& p1, const Point3d& p2, 
-                                   NgArray<int> & locels) const
+                                   Array<ElementIndex> & locels) const
   {
-    elementsearchtree->GetIntersecting (p1, p2, locels);
+    elementsearchtree_vol->GetIntersecting (p1, p2, locels);
   }
 
   void Mesh :: SplitIntoParts()
@@ -7441,9 +7413,10 @@ namespace netgen
   }
 
   string Mesh :: defaultmat = "default";
+  string_view Mesh :: defaultmat_sv = "default";  
   const string & Mesh :: GetMaterial (int domnr) const
   {
-    if (domnr <= materials.Size())
+    if (domnr <= materials.Size() && materials[domnr-1])
       return *materials[domnr-1];
     static string emptystring("default");
     return emptystring;

libsrc/meshing/meshclass.hpp

@@ -128,6 +128,13 @@ namespace netgen
     */
     NgArray<EdgeDescriptor> edgedecoding;
 
+    Array<string*> region_name_cd[4];
+    Array<string*> & materials = region_name_cd[0];
+    Array<string*> & bcnames   = region_name_cd[1];
+    Array<string*> & cd2names  = region_name_cd[2];
+    Array<string*> & cd3names  = region_name_cd[3];
+
+    /*
     /// sub-domain materials 
     Array<string*> materials;
 
@@ -139,7 +146,8 @@ namespace netgen
 
     /// labels for co dim 3 bbboundary conditions
     Array<string*> cd3names;
-
+    */
+    
     /// Periodic surface, close surface, etc. identifications
     unique_ptr<Identifications> ident;
 
@@ -148,9 +156,10 @@ namespace netgen
     int numvertices;
 
     /// geometric search tree for interval intersection search
-    unique_ptr<BoxTree<3>> elementsearchtree;
+    unique_ptr<BoxTree<3, ElementIndex>> elementsearchtree_vol;
+    unique_ptr<BoxTree<3, SurfaceElementIndex>> elementsearchtree_surf;
     /// time stamp for tree
-    mutable int elementsearchtreets;
+    mutable size_t elementsearchtreets[4];
 
     /// element -> face, element -> edge etc ...
     MeshTopology topology;
@@ -200,11 +209,11 @@ namespace netgen
 
     DLL_HEADER bool PointContainedIn2DElement(const Point3d & p,
 				   double lami[3],
-				   const int element,
+				   SurfaceElementIndex element,
 				   bool consider3D = false) const;
     DLL_HEADER bool PointContainedIn3DElement(const Point3d & p,
 				   double lami[3],
-				   const int element) const;
+				   ElementIndex element) const;
     DLL_HEADER bool PointContainedIn3DElementOld(const Point3d & p,
 				      double lami[3],
 				      const int element) const;
@@ -663,35 +672,48 @@ namespace netgen
 
 
     /// build box-search tree
-    DLL_HEADER void BuildElementSearchTree ();
+    DLL_HEADER void BuildElementSearchTree (int dim);
+    BoxTree<3, ElementIndex>* GetElementSearchTree () const
+    {
+        return elementsearchtree_vol.get();
+    }
+
+    BoxTree<3, SurfaceElementIndex>* GetSurfaceElementSearchTree () const
+    {
+      return elementsearchtree_surf.get();
+    }
 
     void SetPointSearchStartElement(const int el) const {ps_startelement = el;}
 
     /// gives element of point, barycentric coordinates
-    DLL_HEADER int GetElementOfPoint (const netgen::Point<3> & p,
-			   double * lami,
-			   bool build_searchtree = 0,
-			   const int index = -1,
-			   const bool allowindex = true) const;
-    DLL_HEADER int GetElementOfPoint (const netgen::Point<3> & p,
-			   double * lami,
-			   const NgArray<int> * const indices,
-			   bool build_searchtree = 0,
-			   const bool allowindex = true) const;
-    DLL_HEADER int GetSurfaceElementOfPoint (const netgen::Point<3> & p,
-				  double * lami,
-				  bool build_searchtree = 0,
-				  const int index = -1,
-				  const bool allowindex = true) const;
-    DLL_HEADER int GetSurfaceElementOfPoint (const netgen::Point<3> & p,
-				  double * lami,
-				  const NgArray<int> * const indices,
-				  bool build_searchtree = 0,
-				  const bool allowindex = true) const;
+    DLL_HEADER ElementIndex
+    GetElementOfPoint (const netgen::Point<3> & p,
+                       double * lami,
+                       bool build_searchtree = false,
+                       int index = -1,
+                       bool allowindex = true) const;
+    DLL_HEADER ElementIndex
+    GetElementOfPoint (const netgen::Point<3> & p,
+                       double * lami,
+                       std::optional<FlatArray<int>> indices,
+                       bool build_searchtree = 0,
+                       bool allowindex = true) const;
+    DLL_HEADER SurfaceElementIndex
+    GetSurfaceElementOfPoint (const netgen::Point<3> & p,
+                              double * lami,
+                              bool build_searchtree = false,
+                              int index = -1,
+                              bool allowindex = true) const;
+    DLL_HEADER SurfaceElementIndex
+    GetSurfaceElementOfPoint (const netgen::Point<3> & p,
+                              double * lami,
+                              std::optional<FlatArray<int>> indices,
+                              bool build_searchtree = false,
+                              bool allowindex = true) const;
 
     /// give list of vol elements which are int the box(p1,p2)
     void GetIntersectingVolEls(const Point3d& p1, const Point3d& p2, 
-			       NgArray<int> & locels) const;
+			       Array<ElementIndex> & locels) const;
 
     ///
     int AddFaceDescriptor(const FaceDescriptor& fd)
@@ -761,6 +783,15 @@ namespace netgen
     std::string_view GetRegionName(SegmentIndex ei) const { return GetRegionName((*this)[ei]); }
     std::string_view GetRegionName(SurfaceElementIndex ei) const { return GetRegionName((*this)[ei]); }
     std::string_view GetRegionName(ElementIndex ei) const { return GetRegionName((*this)[ei]); }
+
+    DLL_HEADER static string_view defaultmat_sv;
+    std::string_view GetRegionName (int dim, int domnr) // 1-based domnr
+    {
+      domnr--;
+      auto & names = region_name_cd[dimension-dim];
+      if (domnr < names.Size() && names[domnr]) return *names[domnr];
+      return defaultmat_sv;
+    }
     
     ///
     void ClearFaceDescriptors()
@@ -1036,7 +1067,6 @@ namespace netgen
   }
 
   DLL_HEADER void AddFacesBetweenDomains(Mesh & mesh);
-  
 }
 
 #endif // NETGEN_MESHCLASS_HPP

libsrc/meshing/meshfunc.cpp

@@ -783,7 +783,7 @@ namespace netgen
         free_segs.Append(segi);
 
     auto get_nonconforming = [&] (const auto & p2el) {
-      Array<size_t> nonconforming;
+      Array<SegmentIndex> nonconforming;
 
       for (auto segi : free_segs) {
         auto seg = mesh[segi];
@@ -805,34 +805,29 @@ namespace netgen
     auto split_segment = [&] (SegmentIndex segi, const auto & p2el) {
       auto seg = mesh[segi];
       auto p_new = Center(mesh[seg[0]], mesh[seg[1]]);
-      auto pi_new = mesh.AddPoint(p_new);
-      auto seg_new0 = seg;
-      auto seg_new1 = seg;
-      seg_new0[1] = pi_new;
-      seg_new1[0] = pi_new;
-
-      mesh[segi][0] = PointIndex::INVALID;
-      mesh.AddSegment(seg_new0);
-      mesh.AddSegment(seg_new1);
-
       double lam[3];
       ElementIndex ei_start = mesh.GetElementOfPoint(p_new, lam, false, domain);
 
-      if(ei_start == 0) {
+      if(!ei_start.IsValid()) {
         PrintMessage(1, "Could not find volume element with new point");
         return;
       }
-      ei_start -= 1;
+
+      if(mesh[ei_start].IsDeleted())
+        return;
 
       double max_inside = -1.;
-      ElementIndex ei_max_inside = -1;
+      ElementIndex ei_max_inside = ElementIndex::INVALID;
 
       // search for adjacent volume element, where the new point is "most inside",
       // i.e. the minimal barycentric coordinate is maximal
       for(auto pi : mesh[ei_start].PNums()) {
         for(auto ei1 : p2el[pi]) {
           double lam[3];
-          if(!mesh.PointContainedIn3DElement(p_new, lam, ei1+1))
+
+          if(mesh[ei1].IsDeleted())
+            return;
+          if(!mesh.PointContainedIn3DElement(p_new, lam, ei1))
             continue;
 
           double inside = min(min(lam[0], lam[1]), min(lam[2], 1.0-lam[0]-lam[1]));
@@ -843,18 +838,34 @@ namespace netgen
         }
       }
 
+      if(max_inside < 1e-4) {
+        PrintMessage(3, "Could not find volume element with new point inside");
+        return;
+      }
+
       // split tet into 4 new tests, with new point inside
       auto el = mesh[ei_max_inside];
       if(el.GetNP() != 4) {
-        PrintMessage(1, "Only tet elements are supported to split around free segments");
+        PrintMessage(3, "Only tet elements are supported to split around free segments");
         return;
       }
 
       if(el.IsDeleted()) {
-        PrintMessage(1,"Element to split is already deleted");
+        PrintMessage(3,"Element to split is already deleted");
         return;
       }
 
+      auto pi_new = mesh.AddPoint(p_new);
+      auto seg_new0 = seg;
+      auto seg_new1 = seg;
+      seg_new0[1] = pi_new;
+      seg_new1[0] = pi_new;
+
+      mesh[segi][0] = PointIndex::INVALID;
+      mesh.AddSegment(seg_new0);
+      mesh.AddSegment(seg_new1);
+
+
       int pmap[4][4] = {
         {0,1,2,4},
         {1,3,2,4},
@@ -897,7 +908,7 @@ namespace netgen
 
       for ([[maybe_unused]] auto i : Range(3)) {
         optmesh.ImproveMesh();
-        optmesh.SwapImprove2 ();
+        optmesh.SwapImprove2(true);
         optmesh.ImproveMesh();
         optmesh.SwapImprove();
         optmesh.ImproveMesh();
@@ -910,7 +921,7 @@ namespace netgen
     auto bad_segs = get_nonconforming(p2el);
 
     if(bad_segs.Size() > 0) {
-      auto bad_seg = mesh[free_segs[bad_segs[0]]];
+      auto bad_seg = mesh[bad_segs[0]];
       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));

libsrc/meshing/meshtype.hpp

@@ -1635,7 +1635,7 @@ namespace netgen
       Point<3> pnt;
       double h;
       int layer = 1;
-      MeshSizePoint (Point<3> _pnt, double _h) : pnt(_pnt), h(_h) { ; }
+      MeshSizePoint (Point<3> pnt_, double h_, int layer_ = 1) : pnt(pnt_), h(h_), layer(layer_) { ; }
       MeshSizePoint () = default;
       MeshSizePoint (const MeshSizePoint &) = default;
       MeshSizePoint (MeshSizePoint &&) = default;

libsrc/meshing/parallelmesh.cpp

@@ -966,7 +966,7 @@ namespace netgen
     self.ident = make_unique<Identifications> (self);
     self.topology = MeshTopology(*this);
     self.topology.Update();
-    self.BuildElementSearchTree();
+    self.BuildElementSearchTree(3);
     
     // const_cast<Mesh&>(*this).DeleteMesh();
 

libsrc/meshing/python_mesh.cpp

@@ -1300,6 +1300,10 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     .def("IdentifyPeriodicBoundaries", &Mesh::IdentifyPeriodicBoundaries,
          py::arg("identification_name"), py::arg("face1"), py::arg("mapping"),
 py::arg("point_tolerance") = -1.)
+    .def("GetCurveOrder", [] (Mesh & self)
+          {
+            return self.GetCurvedElements().GetOrder();
+          })
     .def("GetNrIdentifications", [](Mesh& self)
                                  {
                                    return self.GetIdentifications().GetMaxNr();
@@ -1480,7 +1484,8 @@ py::arg("point_tolerance") = -1.)
            }))
     */
     
-    .def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard<py::gil_scoped_release>())
+    .def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard<py::gil_scoped_release>(),
+          py::arg("dim")=3)
 
     .def ("BoundaryLayer2", GenerateBoundaryLayer2, py::arg("domain"), py::arg("thicknesses"), py::arg("make_new_domain")=true, py::arg("boundaries")=Array<int>{})
     .def ("BoundaryLayer", [](Mesh & self, variant<string, int, std::vector<int>> boundary,
@@ -1682,25 +1687,25 @@ py::arg("point_tolerance") = -1.)
                     return mp;
                   }), py::arg("mp")=nullptr, meshingparameter_description.c_str())
     .def("__str__", &ToString<MP>)
-    .def("RestrictH", [](MP & mp, double x, double y, double z, double h)
+    .def("RestrictH", [](MP & mp, double x, double y, double z, double h, int layer)
           {
-            mp.meshsize_points.Append ( MeshingParameters::MeshSizePoint(Point<3> (x,y,z), h));
-          }, py::arg("x"), py::arg("y"), py::arg("z"), py::arg("h")
+            mp.meshsize_points.Append ( MeshingParameters::MeshSizePoint(Point<3> (x,y,z), h, layer));
+          }, py::arg("x"), py::arg("y"), py::arg("z"), py::arg("h"), py::arg("layer")=1
          )
-    .def("RestrictH", [](MP & mp, const Point<3>& p, double h)
+    .def("RestrictH", [](MP & mp, const Point<3>& p, double h, int layer)
     {
-      mp.meshsize_points.Append ({p, h});
-    }, py::arg("p"), py::arg("h"))
+      mp.meshsize_points.Append ({p, h, layer});
+    }, py::arg("p"), py::arg("h"), py::arg("layer")=1)
     .def("RestrictHLine", [](MP& mp, const Point<3>& p1, const Point<3>& p2,
-                             double maxh)
+                             double maxh, int layer)
     {
       int steps = int(Dist(p1, p2) / maxh) + 2;
       auto v = p2 - p1;
       for (int i = 0; i <= steps; i++)
         {
-          mp.meshsize_points.Append({p1 + double(i)/steps * v, maxh});
+          mp.meshsize_points.Append({p1 + double(i)/steps * v, maxh, layer});
         }
-    }, py::arg("p1"), py::arg("p2"), py::arg("maxh"))
+    }, py::arg("p1"), py::arg("p2"), py::arg("maxh"), py::arg("layer")=1)
     ;
 
   m.def("SetTestoutFile", FunctionPointer ([] (const string & filename)

libsrc/occ/occgeom.cpp

@@ -2288,10 +2288,22 @@ namespace netgen
             XCAFPrs::CollectStyleSettings(label, loc, set);
             XCAFPrs_Style aStyle;
             set.FindFromKey(e.Current(), aStyle);
-
-            auto & prop = OCCGeometry::GetProperties(e.Current());
             if(aStyle.IsSetColorSurf())
-                prop.col = step_utils::ReadColor(aStyle.GetColorSurfRGBA());
+              {
+                for(TopExp_Explorer e2(e.Current(), TopAbs_FACE); e2.More(); e2.Next())
+                  {
+                    auto & prop = OCCGeometry::GetProperties(e2.Current());
+                    prop.col = step_utils::ReadColor(aStyle.GetColorSurfRGBA());
+                  }
+              }
+            if(aStyle.IsSetColorCurv())
+              {
+                for(TopExp_Explorer e2(e.Current(), TopAbs_EDGE); e2.More(); e2.Next())
+                  {
+                    auto & prop = OCCGeometry::GetProperties(e2.Current());
+                    prop.col = step_utils::ReadColor(aStyle.GetColorSurfRGBA());
+                  }
+              }
           }
 
         // load names

libsrc/occ/python_occ.cpp

@@ -21,6 +21,9 @@
 #include <XCAFDoc_DocumentTool.hxx>
 #include <XCAFDoc_MaterialTool.hxx>
 #include <XCAFDoc_ShapeTool.hxx>
+#include <TopoDS_Edge.hxx>
+#include <BRepAdaptor_Curve.hxx>
+#include <GCPnts_TangentialDeflection.hxx>
 
 using namespace netgen;
 
@@ -165,10 +168,40 @@ DLL_HEADER void ExportNgOCC(py::module &m)
     {
       ng_geometry = geo;
     })
+    .def_property_readonly("solids", [](shared_ptr<OCCGeometry> geo)
+         {
+           ListOfShapes solids;
+           for (int i = 1; i <= geo->somap.Extent(); i++)
+             solids.push_back(geo->somap(i));
+           return solids;
+         }, "Get solids in order that they will be in the mesh")
+    .def_property_readonly("faces", [](shared_ptr<OCCGeometry> geo)
+         {
+           ListOfShapes faces;
+           for (int i = 1; i <= geo->fmap.Extent(); i++)
+             faces.push_back(geo->fmap(i));
+           return faces;
+         }, "Get faces in order that they will be in the mesh")
+    .def_property_readonly("edges", [](shared_ptr<OCCGeometry> geo)
+         {
+           ListOfShapes edges;
+           for (int i = 1; i <= geo->emap.Extent(); i++)
+             edges.push_back(geo->emap(i));
+           return edges;
+         }, "Get edges in order that they will be in the mesh")
+    .def_property_readonly("vertices", [](shared_ptr<OCCGeometry> geo)
+         {
+           ListOfShapes vertices;
+           for (int i = 1; i <= geo->vmap.Extent(); i++)
+             vertices.push_back(geo->vmap(i));
+           return vertices;
+         }, "Get vertices in order that they will be in the mesh")
     .def("_visualizationData", [] (shared_ptr<OCCGeometry> occ_geo)
          {
            std::vector<float> vertices;
-           std::vector<int> trigs;
+           std::vector<uint32_t> indices;
+           std::vector<float> edges;
+           std::vector<uint32_t> edge_indices;
            std::vector<float> normals;
            std::vector<float> min = {std::numeric_limits<float>::max(),
                                std::numeric_limits<float>::max(),
@@ -176,7 +209,8 @@ DLL_HEADER void ExportNgOCC(py::module &m)
            std::vector<float> max = {std::numeric_limits<float>::lowest(),
                                std::numeric_limits<float>::lowest(),
                                std::numeric_limits<float>::lowest()};
-           std::vector<string> surfnames;
+           std::vector<float> face_colors;
+           std::vector<float> edge_colors;
            auto box = occ_geo->GetBoundingBox();
            for(int i = 0; i < 3; i++)
              {
@@ -188,11 +222,67 @@ DLL_HEADER void ExportNgOCC(py::module &m)
            gp_Pnt pnt;
            gp_Vec n;
            gp_Pnt p[3];
-           int count = 0;
+           for(int edge_index = 1; edge_index <= occ_geo->emap.Extent();
+               edge_index++)
+             {
+               auto edge = TopoDS::Edge(occ_geo->emap(edge_index));
+               if(OCCGeometry::HaveProperties(edge))
+                 {
+                   const auto& props = OCCGeometry::GetProperties(edge);
+                   if(props.col)
+                     edge_colors.insert(edge_colors.end(),
+                                        {float((*props.col)[0]),
+                                         float((*props.col)[1]),
+                                         float((*props.col)[2]),
+                                         float((*props.col)[3])});
+                   else
+                     edge_colors.insert(edge_colors.end(),{0.f,0.f,0.f,1.f});
+                 }
+               else
+                 {
+                   edge_colors.insert(edge_colors.end(),{0.f,0.f,0.f,1.f});
+                 }
+               BRepAdaptor_Curve adapt_crv = BRepAdaptor_Curve(edge);
+               GCPnts_TangentialDeflection discretizer;
+               discretizer.Initialize(adapt_crv, 0.09, 0.01);
+		if (discretizer.NbPoints() > 1)
+		{
+                  for (int j = 1; j <= discretizer.NbPoints()-1; ++j)
+                    {
+                      gp_Pnt p_0 = discretizer.Value(j);
+                      gp_Pnt p_1 = discretizer.Value(j+1);
+                      edges.insert(edges.end(),
+                                   {float(p_0.X()),
+                                    float(p_0.Y()),
+                                    float(p_0.Z()),
+                                    float(p_1.X()),
+                                    float(p_1.Y()),
+                                    float(p_1.Z())});
+                      edge_indices.push_back(uint32_t(edge_index-1));
+                    }
+		}
+             }
            for (int i = 1; i <= occ_geo->fmap.Extent(); i++)
              {
-               surfnames.push_back("occ_surface" + to_string(i));
                auto face = TopoDS::Face(occ_geo->fmap(i));
+               if (OCCGeometry::HaveProperties(face))
+                 {
+                   const auto& props = OCCGeometry::GetProperties(face);
+                   if(props.col)
+                     face_colors.insert(face_colors.end(),
+                                        {float((*props.col)[0]),
+                                         float((*props.col)[1]),
+                                         float((*props.col)[2]),
+                                         float((*props.col)[3])});
+                    else
+                      {
+                        face_colors.insert(face_colors.end(),{0.7,0.7,0.7,1.});
+                      }
+                 }
+               else
+                 {
+                   face_colors.insert(face_colors.end(),{0.7,0.7,0.7,1.});
+                 }
                auto surf = BRep_Tool::Surface(face);
                TopLoc_Location loc;
                BRepAdaptor_Surface sf(face, Standard_False);
@@ -200,7 +290,7 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc);
                if (triangulation.IsNull())
                  cout << "cannot visualize face " << i << endl;
-               trigs.reserve(trigs.size() + triangulation->NbTriangles()*4);
+               indices.reserve(indices.size() + triangulation->NbTriangles());
                vertices.reserve(vertices.size() + triangulation->NbTriangles()*3*3);
                normals.reserve(normals.size() + triangulation->NbTriangles()*3*3);
                for (int j = 1; j < triangulation->NbTriangles()+1; j++)
@@ -208,11 +298,13 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                    auto triangle = triangulation->Triangle(j);
                    for (int k = 1; k < 4; k++)
                      p[k-1] = triangulation->Node(triangle(k)).Transformed(loc);
+                   indices.push_back(uint32_t(i-1));
                    for (int k = 1; k < 4; k++)
                      {
-                       vertices.insert(vertices.end(),{float(p[k-1].X()), float(p[k-1].Y()), float(p[k-1].Z())});
-                       trigs.insert(trigs.end(),{count, count+1, count+2,i});
-                       count += 3;
+                       vertices.insert(vertices.end(),{
+                           float(p[k-1].X()),
+                           float(p[k-1].Y()),
+                           float(p[k-1].Z())});
                        uv = triangulation->UVNode(triangle(k));
                        prop.SetParameters(uv.X(), uv.Y());
                        if (prop.IsNormalDefined())
@@ -231,12 +323,13 @@ DLL_HEADER void ExportNgOCC(py::module &m)
             py::gil_scoped_acquire ac;
             py::dict res;
             py::list snames;
-            for(auto name : surfnames)
-              snames.append(py::cast(name));
             res["vertices"] = MoveToNumpy(vertices);
-            res["triangles"] = MoveToNumpy(trigs);
+            res["edges"] = MoveToNumpy(edges);
+            res["edge_indices"] = MoveToNumpy(edge_indices);
+            res["edge_colors"] = MoveToNumpy(edge_colors);
+            res["indices"] = MoveToNumpy(indices);
             res["normals"] = MoveToNumpy(normals);
-            res["surfnames"] = snames;
+            res["face_colors"] = MoveToNumpy(face_colors);
             res["min"] = MoveToNumpy(min);
             res["max"] = MoveToNumpy(max);
             return res;

libsrc/occ/python_occ_shapes.cpp

@@ -14,6 +14,9 @@
 #pragma clang diagnostic push
 #pragma clang diagnostic ignored "-Wdeprecated-declarations"
 
+#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=6
+#include <BinTools_ShapeWriter.hxx>
+#endif // OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4
 #include <BOPAlgo_Builder.hxx>
 #include <BOPTools_AlgoTools.hxx>
 #include <BRepAlgoAPI_Common.hxx>
@@ -49,6 +52,7 @@
 #include <BRepTools.hxx>
 #include <GCE2d_MakeArcOfCircle.hxx>
 #include <GCE2d_MakeCircle.hxx>
+#include <GCE2d_MakeEllipse.hxx>
 #include <GCE2d_MakeSegment.hxx>
 #include <GC_MakeArcOfCircle.hxx>
 #include <GC_MakeCircle.hxx>
@@ -598,6 +602,19 @@ public:
     return shared_from_this();    
   }
 
+  auto Ellipse(double major, double minor)
+  {
+    Handle(Geom2d_Ellipse) ell_curve = GCE2d_MakeEllipse(localpos, major, minor).Value();
+
+    auto edge = BRepBuilderAPI_MakeEdge(ell_curve, surf).Edge();
+    BRepLib::BuildCurves3d(edge);
+
+    wire_builder.Add(edge);
+    wires.push_back (wire_builder.Wire());
+    wire_builder = BRepBuilderAPI_MakeWire();
+    return shared_from_this();
+  }
+
   auto NameVertex (string name)
   {
     if (!lastvertex.IsNull())
@@ -689,6 +706,32 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
     OCCGeometry::global_shape_properties.clear();
     OCCGeometry::global_shape_property_indices.Clear();
   });
+
+  struct SwigTypeInfo
+  {
+    const char* name;  // SWIG's type name string
+    // Other fields...
+  };
+
+  struct SwigPyObject{
+    PyObject_HEAD
+    void *ptr;
+    SwigTypeInfo* ty; // SWIG type information
+    int own; // ownership flag
+  };
+
+  m.def("From_PyOCC", [](py::object shape)
+  {
+    py::object py_this = shape.attr("this");
+    PyObject* obj = py_this.ptr();
+    SwigPyObject* swig_obj = reinterpret_cast<SwigPyObject*>(obj);
+    if (!swig_obj->ptr || !swig_obj->ty || !swig_obj->ty->name) {
+        throw std::runtime_error("SWIG object does not contain a valid pointer");
+    }
+    if(strcmp(swig_obj->ty->name, "_p_TopoDS_Shape") != 0)
+      throw std::runtime_error("Does not contain TopoDS_Shape from pyocc!");
+    return py::cast(static_cast<TopoDS_Shape*>(swig_obj->ptr));
+  }, py::return_value_policy::reference, py::keep_alive<0,1>());
   
   py::class_<TopoDS_Shape> (m, "TopoDS_Shape")
     .def("__str__", [] (const TopoDS_Shape & shape)
@@ -807,10 +850,33 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
     .def("WriteStep", [](const TopoDS_Shape & shape, string & filename)
             { step_utils::WriteSTEP(shape, filename); }
          , py::arg("filename"), "export shape in STEP - format")
-    .def("WriteBrep", [](const TopoDS_Shape & shape, const string& filename)
+    .def("WriteBrep", [](const TopoDS_Shape & shape, const string& filename,
+                         bool withTriangles, bool withNormals,
+                         optional<int> version, bool binary)
     {
-      BRepTools::Write(shape, filename.c_str());
-    }, py::arg("filename"), "export shape in BREP - format")
+      if(binary)
+        {
+#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=6
+          BinTools_FormatVersion v = version ? BinTools_FormatVersion(*version) : BinTools_FormatVersion_CURRENT;
+          BinTools::Write(shape, filename.c_str(), withTriangles, withNormals, v);
+# else // OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=6
+          throw Exception("Binary BREP export not supported in this version of OpenCascade");
+#endif // OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=6
+        }
+      else
+        {
+#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=6
+          TopTools_FormatVersion v = version ? (TopTools_FormatVersion)(*version) : TopTools_FormatVersion_CURRENT;
+          BRepTools::Write(shape, filename.c_str(), withTriangles, withNormals, v);
+#else // OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=6
+          BRepTools::Write(shape, filename.c_str());
+#endif // OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=6
+        }
+    }, py::arg("filename"), py::arg("withTriangles")=true,
+       py::arg("withNormals")=false,
+         py::arg("version")=py::none(),
+         py::arg("binary")=false,
+       "export shape in BREP - format")
     .def("bc", [](const TopoDS_Shape & shape, const string & name)
          {
            for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next())
@@ -1898,21 +1964,21 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
       })
     .def("Edge", [](Handle(Geom2d_Curve) curve) {
         // static Geom_Plane surf{gp_Ax3()}; // crashes in nbconvert ???
-        static auto surf = new Geom_Plane{gp_Ax3()};          
+        static auto surf = Handle(Geom_Plane)(new Geom_Plane{gp_Ax3()});
         auto edge = BRepBuilderAPI_MakeEdge(curve, surf).Edge();
         BRepLib::BuildCurves3d(edge);
         return edge;
       })
     .def("Wire", [](Handle(Geom2d_Curve) curve) {
         // static Geom_Plane surf{gp_Ax3()}; // crashes in nbconvert ???
-        static auto surf = new Geom_Plane{gp_Ax3()};          
+        static auto surf = Handle(Geom_Plane)(new Geom_Plane{gp_Ax3()});
         auto edge = BRepBuilderAPI_MakeEdge(curve, surf).Edge();
         BRepLib::BuildCurves3d(edge);
         return BRepBuilderAPI_MakeWire(edge).Wire();                
       })
     .def("Face", [](Handle(Geom2d_Curve) curve) {
         // static Geom_Plane surf{gp_Ax3()};  // crashes in nbconvert ???
-        static auto surf = new Geom_Plane{gp_Ax3()};
+        static auto surf = Handle(Geom_Plane)(new Geom_Plane{gp_Ax3()});
         auto edge = BRepBuilderAPI_MakeEdge(curve, surf).Edge();
         BRepLib::BuildCurves3d(edge);        
         auto wire = BRepBuilderAPI_MakeWire(edge).Wire();        
@@ -2046,14 +2112,19 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
           return BRepOffsetAPI_MakePipe (spine, profile).Shape();
         }, py::arg("spine"), py::arg("profile"), py::arg("twist")=nullopt, py::arg("auxspine")=nullopt);
   
-  m.def("PipeShell", [] (const TopoDS_Wire & spine, const TopoDS_Shape & profile, const TopoDS_Wire & auxspine) {
+  m.def("PipeShell", [] (const TopoDS_Wire & spine, variant<TopoDS_Shape, std::vector<TopoDS_Shape>> profile, std::optional<TopoDS_Wire> auxspine) {
       try
         {
           BRepOffsetAPI_MakePipeShell builder(spine);
-          builder.SetMode (auxspine, Standard_True);
-          builder.Add (profile);
-          // builder.Build();
-          // builder.MakeSolid();
+          if(auxspine)
+            builder.SetMode (*auxspine, Standard_True);
+          if(std::holds_alternative<TopoDS_Shape>(profile))
+            builder.Add (std::get<TopoDS_Shape>(profile));
+          else
+            {
+              for(auto s : std::get<std::vector<TopoDS_Shape>>(profile))
+                builder.Add(s);
+            }
           return builder.Shape();
         }
       catch (Standard_Failure & e)
@@ -2062,13 +2133,13 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
           e.Print(errstr);
           throw NgException("cannot create PipeShell: "+errstr.str());
         }
-    }, py::arg("spine"), py::arg("profile"), py::arg("auxspine"));
+    }, py::arg("spine"), py::arg("profile"), py::arg("auxspine")=nullopt);
 
 
   // Handle(Geom2d_Ellipse) anEllipse1 = new Geom2d_Ellipse(anAx2d, aMajor, aMinor);
   m.def("Ellipse", [] (const gp_Ax2d & ax, double major, double minor) -> Handle(Geom2d_Curve)
         {
-          return new Geom2d_Ellipse(ax, major, minor);
+          return Handle(Geom2d_Ellipse) (GCE2d_MakeEllipse(ax, major, minor));
         }, py::arg("axes"), py::arg("major"), py::arg("minor"), "create 2d ellipse curve");
   
   m.def("Segment", [](gp_Pnt2d p1, gp_Pnt2d p2) -> Handle(Geom2d_Curve) {
@@ -2636,6 +2707,8 @@ degen_tol : double
     .def("Circle", [](WorkPlane&wp, double x, double y, double r) {
         return wp.Circle(x,y,r); }, py::arg("h"), py::arg("v"), py::arg("r"), "draw circle with center (h,v) and radius 'r'")
     .def("Circle", [](WorkPlane&wp, double r) { return wp.Circle(r); }, py::arg("r"), "draw circle with center in current position")
+    .def("Ellipse", [](WorkPlane& wp, double major, double minor)
+    { return wp.Ellipse(major, minor); }, py::arg("major"), py::arg("minor"), "draw ellipse with current position as center")
     .def("NameVertex", &WorkPlane::NameVertex, py::arg("name"), "name vertex at current position")
     .def("Offset", &WorkPlane::Offset, py::arg("d"), "replace current wire by offset curve of distance 'd'")
     .def("Reverse", &WorkPlane::Reverse, "revert orientation of current wire")

libsrc/visualization/mvdraw.cpp

@@ -57,6 +57,13 @@ namespace netgen
       return opengl_text_width;
   }
 
+  void MyOpenGLLines(FlatArray<Point<3>> points)
+  {
+    glBegin(GL_LINES);
+    for (auto p : points)
+      glVertex3dv(&p[0]);
+    glEnd();
+  }
 
   // texture for color decoding
   // GLubyte * VisualScene :: colortexture = NULL;

libsrc/visualization/mvdraw.hpp

@@ -90,6 +90,7 @@ namespace netgen
   NGGUI_API extern void Set_OpenGLText_Callback ( void (*fun) (const char * text), int width );
   NGGUI_API extern VisualScene visual_scene_cross;
   NGGUI_API extern VisualScene *visual_scene;
+  NGGUI_API extern void MyOpenGLLines (FlatArray<Point<3>> points);
 
 
 

libsrc/visualization/vssolution.cpp

@@ -841,7 +841,7 @@ namespace netgen
         if (vispar.clipping.enable && clipsolution == 2)      
           {
             mesh->Mutex().unlock();
-            mesh->BuildElementSearchTree();
+            mesh->BuildElementSearchTree(3);
             mesh->Mutex().lock();
           }
 
@@ -2667,6 +2667,8 @@ namespace netgen
                   for (int i=first; i<next; i++)
                     {
                       double val;
+                      if(!VolumeElementActive(sol, *mesh, (*mesh)[ElementIndex(i)]))
+                        continue;
                       bool considerElem = GetValue (sol, i, 0.333, 0.333, 0.333, comp, val);
                       if (considerElem)
                         {
@@ -2698,6 +2700,8 @@ namespace netgen
               // for (int i = 0; i < nse; i++)
               for (SurfaceElementIndex i : mesh->SurfaceElements().Range())
                 {
+                  if(!SurfaceElementActive(sol, *mesh, (*mesh)[i]))
+                    continue;
                   ELEMENT_TYPE type = (*mesh)[i].GetType();
                   double val;
                   bool considerElem = (type == QUAD) 
@@ -4337,7 +4341,7 @@ namespace netgen
             }
 
           double lami[3];
-          int elnr = mesh->GetElementOfPoint (hp, lami,0,cindex,allowindex)-1;
+          int elnr = mesh->GetElementOfPoint (hp, lami,0,cindex,allowindex);
 
           if (elnr != -1)
             {
@@ -4723,7 +4727,7 @@ namespace netgen
 
 
   bool VisualSceneSolution ::
-  SurfaceElementActive(const SolData *data, const Mesh & mesh, const Element2d & el)
+  SurfaceElementActive(const SolData *data, const Mesh & mesh, const Element2d & el) const
   {
     if(data == nullptr) return true;
     bool is_active = true;
@@ -4749,7 +4753,7 @@ namespace netgen
   }
 
   bool VisualSceneSolution ::
-  VolumeElementActive(const SolData *data, const Mesh & mesh, const Element & el)
+  VolumeElementActive(const SolData *data, const Mesh & mesh, const Element & el) const
   {
     bool is_active = true;
     if(data->draw_volumes)
@@ -4858,18 +4862,18 @@ namespace netgen
             if(sol.iscomplex && rcomponent != 0)
             {
               rcomponent = 2 * ((rcomponent-1)/2) + 1;
-              GetValue(&sol, el3d-1,  lami[0], lami[1], lami[2], rcomponent+1,
+              GetValue(&sol, el3d,  lami[0], lami[1], lami[2], rcomponent+1,
                   imag);
               comp = (scalcomp-1)/2 + 1;
             }
-            GetValue(&sol, el3d-1,  lami[0], lami[1], lami[2], rcomponent, val);
+            GetValue(&sol, el3d,  lami[0], lami[1], lami[2], rcomponent, val);
             printScalValue(sol, comp, val, imag, sol.iscomplex && comp > 0);
           }
           if(vecfunction!=-1 && soldata[vecfunction]->draw_volume)
           {
             auto & sol = *soldata[vecfunction];
             ArrayMem<double, 10> values(sol.components);
-            GetValues(&sol, el3d-1,  lami[0], lami[1], lami[2], &values[0]);
+            GetValues(&sol, el3d,  lami[0], lami[1], lami[2], &values[0]);
             printVecValue(sol, values);
           }
           return;
@@ -4880,7 +4884,7 @@ namespace netgen
     double lami[3] = {0.0, 0.0, 0.0};
     // Check if unprojected Point is close to surface element (eps of 1e-3 due to z-Buffer accuracy)
     bool found_2del = false;
-    if(selelement>0 && mesh->PointContainedIn2DElement(p, lami, selelement, false && fabs(lami[2])<1e-3))
+    if(selelement>0 && mesh->PointContainedIn2DElement(p, lami, selelement-1, false && fabs(lami[2])<1e-3))
       {
         // Found it, use coordinates of point projected to surface element
         mesh->GetCurvedElements().CalcSurfaceTransformation({1.0-lami[0]-lami[1], lami[0]}, selelement-1, p);

libsrc/visualization/vssolution.hpp

@@ -253,8 +253,8 @@ private:
   void DrawCone (const Point<3> & p1, const Point<3> & p2, double r);
   void DrawCylinder (const Point<3> & p1, const Point<3> & p2, double r);
 
-  bool SurfaceElementActive(const SolData *data, const Mesh & mesh, const Element2d & sei);
-  bool VolumeElementActive(const SolData *data, const Mesh & mesh, const Element & ei);
+  bool SurfaceElementActive(const SolData *data, const Mesh & mesh, const Element2d & sei) const;
+  bool VolumeElementActive(const SolData *data, const Mesh & mesh, const Element & ei) const;
 
   // Get Function Value, local coordinates lam1, lam2, lam3, 
   bool GetValue (const SolData * data, ElementIndex elnr,