Merge remote-tracking branch 'origin/master' into StefanBruens-optionally_use_system_pybind11

CMakeLists.txt

@@ -250,7 +250,7 @@ target_include_directories(nglib PRIVATE ${ZLIB_INCLUDE_DIRS})
 if(USE_GUI)
     target_include_directories(nggui PRIVATE ${ZLIB_INCLUDE_DIRS})
 endif(USE_GUI)
-target_link_libraries(nglib PUBLIC ${ZLIB_LIBRARIES})
+target_link_libraries(nglib PRIVATE ${ZLIB_LIBRARIES})
 
 #######################################################################
 if(WIN32)
@@ -415,9 +415,9 @@ if (USE_OCC)
       endif()
     endif()
     message(STATUS "OCC DIRS ${OpenCASCADE_INCLUDE_DIR}")
-    if(WIN32)
+    if(WIN32 AND USE_GUI)
       target_link_libraries(nggui PRIVATE ${OCC_LIBRARIES})
-    endif(WIN32)
+    endif(WIN32 AND USE_GUI)
 endif (USE_OCC)
 
 #######################################################################

README.md

@@ -1,3 +1,6 @@
 Netgen mesh generator
 
 NETGEN is an automatic 3d tetrahedral mesh generator. It accepts input from constructive solid geometry (CSG) or boundary representation (BRep) from STL file format. The connection to a geometry kernel allows the handling of IGES and STEP files. NETGEN contains modules for mesh optimization and hierarchical mesh refinement. Netgen 6.x supports scripting via a Python interface. Netgen is open source based on the LGPL license. It is available for Unix/Linux, Windows, and OSX.
+
+Find the Open Source Community on https://ngsolve.org
+Support & Services: https://cerbsim.com

libsrc/core/python_ngcore_export.cpp

@@ -240,7 +240,10 @@ threads : int
   class ParallelContextManager {
       int num_threads;
     public:
-      ParallelContextManager() : num_threads(0) {};
+      ParallelContextManager() : num_threads(0) {
+        TaskManager::SetPajeTrace(0);
+        PajeTrace::SetMaxTracefileSize(0);
+      };
       ParallelContextManager(size_t pajesize) : num_threads(0) {
         TaskManager::SetPajeTrace(pajesize > 0);
         PajeTrace::SetMaxTracefileSize(pajesize);

libsrc/core/simd_avx512.hpp

@@ -143,7 +143,7 @@ namespace ngcore
     }
   };
 
-  NETGEN_INLINE SIMD<double,8> operator- (SIMD<double,8> a) { return -a.Data(); }
+  NETGEN_INLINE SIMD<double,8> operator- (SIMD<double,8> a) { return _mm512_xor_pd(a.Data(), _mm512_set1_pd(-0.0)); } //{ return -a.Data(); }
   NETGEN_INLINE SIMD<double,8> operator+ (SIMD<double,8> a, SIMD<double,8> b) { return _mm512_add_pd(a.Data(),b.Data()); }
   NETGEN_INLINE SIMD<double,8> operator- (SIMD<double,8> a, SIMD<double,8> b) { return _mm512_sub_pd(a.Data(),b.Data()); }
   NETGEN_INLINE SIMD<double,8> operator* (SIMD<double,8> a, SIMD<double,8> b) { return _mm512_mul_pd(a.Data(),b.Data()); }
@@ -154,7 +154,7 @@ namespace ngcore
   NETGEN_INLINE SIMD<double,8> sqrt (SIMD<double,8> a) { return _mm512_sqrt_pd(a.Data()); }
   NETGEN_INLINE SIMD<double,8> floor (SIMD<double,8> a) { return _mm512_floor_pd(a.Data()); }
   NETGEN_INLINE SIMD<double,8> ceil (SIMD<double,8> a) { return _mm512_ceil_pd(a.Data()); }
-  NETGEN_INLINE SIMD<double,8> fabs (SIMD<double,8> a) { return _mm512_max_pd(a.Data(), -a.Data()); }
+  NETGEN_INLINE SIMD<double,8> fabs (SIMD<double,8> a) { return _mm512_max_pd(a.Data(), ( - a).Data()); }
 
   NETGEN_INLINE SIMD<mask64,8> operator<= (SIMD<double,8> a , SIMD<double,8> b)
   { return _mm512_cmp_pd_mask (a.Data(), b.Data(), _CMP_LE_OQ); }
@@ -233,9 +233,9 @@ namespace ngcore
     // sum01  b a b a b a b a
     // sum23  d c d c d c d c
     // __m512 perm = _mm512_permutex2var_pd (sum01.Data(), _mm512_set_epi64(1,2,3,4,5,6,7,8), sum23.Data());
-    __m256d ab =  _mm512_extractf64x4_pd(sum01.Data(),0) + _mm512_extractf64x4_pd(sum01.Data(),1);
+    SIMD<double,4> ab =  _mm512_extractf64x4_pd(sum01.Data(),0) + _mm512_extractf64x4_pd(sum01.Data(),1);
-    __m256d cd =  _mm512_extractf64x4_pd(sum23.Data(),0) + _mm512_extractf64x4_pd(sum23.Data(),1);
+    SIMD<double,4> cd =  _mm512_extractf64x4_pd(sum23.Data(),0) + _mm512_extractf64x4_pd(sum23.Data(),1);
-    return _mm256_add_pd (_mm256_permute2f128_pd (ab, cd, 1+2*16), _mm256_blend_pd (ab, cd, 12));
+    return _mm256_add_pd (_mm256_permute2f128_pd (ab.Data(), cd.Data(), 1 + 2 * 16), _mm256_blend_pd(ab.Data(), cd.Data(), 12));
   }
 
   NETGEN_INLINE SIMD<double,8> FMA (SIMD<double,8> a, SIMD<double,8> b, SIMD<double,8> c)

libsrc/core/table.hpp

@@ -46,8 +46,7 @@ namespace ngcore
     /// Access entry
     NETGEN_INLINE const FlatArray<T> operator[] (IndexType i) const
     {
-      i = i-BASE;
+      return FlatArray<T> (index[i-BASE+1]-index[i-BASE], data+index[i-BASE]);
-      return FlatArray<T> (index[i+1]-index[i], data+index[i]);
     }
 
     NETGEN_INLINE T * Data() const { return data; }

libsrc/csg/csgeom.cpp

@@ -417,7 +417,7 @@ namespace netgen
       & identpoints & boundingbox & isidenticto & ideps
       & filename & spline_surfaces & splinecurves2d & splinecurves3d & surf2prim;
     if(archive.Input())
-      FindIdenticSurfaces(1e-6);
+      FindIdenticSurfaces(1e-8 * MaxSize());
   }
 
   void CSGeometry :: SaveSurfaces (ostream & out) const
@@ -949,6 +949,7 @@ namespace netgen
   {
     int inv;
     int nsurf = GetNSurf();
+    identicsurfaces.DeleteData();
 
 
     isidenticto.SetSize(nsurf);

libsrc/csg/python_csg.cpp

@@ -676,7 +676,7 @@ However, when r = 0, the top part becomes a point(tip) and meshing fails!
     .def("Draw", FunctionPointer
          ([] (shared_ptr<CSGeometry> self)
           {
-             self->FindIdenticSurfaces(1e-6);
+             self->FindIdenticSurfaces(1e-8 * self->MaxSize());
              self->CalcTriangleApproximation(0.01, 20);
              ng_geometry = self;
           })
@@ -706,8 +706,8 @@ However, when r = 0, the top part becomes a point(tip) and meshing fails!
                auto surf = csg_geo->GetSurface(i);
                surfnames.push_back(surf->GetBCName());
              }
-           csg_geo->FindIdenticSurfaces(1e-6);
+           csg_geo->FindIdenticSurfaces(1e-8 * csg_geo->MaxSize());
-           csg_geo->CalcTriangleApproximation(0.01,100);
+           csg_geo->CalcTriangleApproximation(0.01,20);
            auto nto = csg_geo->GetNTopLevelObjects();
            size_t np = 0;
            size_t ntrig = 0;

libsrc/gprim/geomtest3d.hpp

@@ -69,10 +69,10 @@ extern double ComputeCylinderRadius (const Vec3d & n1, const Vec3d & n2,
 				     double h1, double h2);
 
 /// Minimal distance of point p to the line segment [lp1,lp2]
-extern double MinDistLP2 (const Point2d & lp1, const Point2d & lp2, const Point2d & p);
+DLL_HEADER double MinDistLP2 (const Point2d & lp1, const Point2d & lp2, const Point2d & p);
 
 /// Minimal distance of point p to the line segment [lp1,lp2]
-extern double MinDistLP2 (const Point3d & lp1, const Point3d & lp2, const Point3d & p);
+DLL_HEADER double MinDistLP2 (const Point3d & lp1, const Point3d & lp2, const Point3d & p);
 
 /// Minimal distance of point p to the triangle segment [tp1,tp2,pt3]
 DLL_HEADER double MinDistTP2 (const Point3d & tp1, const Point3d & tp2, 

libsrc/gprim/spline.hpp

@@ -188,12 +188,12 @@ namespace netgen
     mutable double proj_latest_t;
   public:
     ///
-    SplineSeg3 (const GeomPoint<D> & ap1, 
+    DLL_HEADER SplineSeg3 (const GeomPoint<D> & ap1,
 		const GeomPoint<D> & ap2, 
 		const GeomPoint<D> & ap3,
                 string bcname="default",
                 double maxh=1e99);
-    SplineSeg3 (const GeomPoint<D> & ap1,
+    DLL_HEADER SplineSeg3 (const GeomPoint<D> & ap1,
 		const GeomPoint<D> & ap2,
 		const GeomPoint<D> & ap3,
                 double aweight,

libsrc/interface/nginterface_v2.cpp

@@ -701,9 +701,9 @@ namespace netgen
       {
 	int hpelnr = -1;
 	if (mesh->GetDimension() == 2)
-	  hpelnr = mesh->SurfaceElement(ei).hp_elnr;
+	  hpelnr = mesh->SurfaceElement(ei).GetHpElnr();
 	else
-	  hpelnr = mesh->VolumeElement(ei).hp_elnr;
+	  hpelnr = mesh->VolumeElement(ei).GetHpElnr();
 
         if (hpelnr < 0)
           throw NgException("Ngx_Mesh::GetHPElementLevel: Wrong hp-element number!");

libsrc/meshing/basegeom.cpp

@@ -785,6 +785,8 @@ namespace netgen
     {
         auto & face = *faces[k];
         FaceDescriptor fd(k+1, face.domin+1, face.domout+1, k+1);
+        if(face.properties.col)
+          fd.SetSurfColour(*face.properties.col);
         mesh.AddFaceDescriptor(fd);
         mesh.SetBCName(k, face.properties.GetName());
         if(face.primary == &face)
@@ -974,6 +976,12 @@ namespace netgen
           }
 
     xbool do_invert = maybe;
+    if(dst.identifications[0].type == Identifications::PERIODIC)
+      {
+        auto other = static_cast<GeometryFace*>(dst.primary);
+        if(dst.domin != other->domout && dst.domout != other->domin)
+          do_invert = true;
+      }
 
     // now insert mapped surface elements
     for(auto sei : mesh.SurfaceElements().Range())
@@ -1050,9 +1058,10 @@ namespace netgen
 
   void NetgenGeometry :: FinalizeMesh(Mesh& mesh) const
   {
-    for (int i = 0; i < mesh.GetNDomains(); i++)
+    if(solids.Size())
-      if (auto name = solids[i]->properties.name)
+      for (int i = 0; i < mesh.GetNDomains(); i++)
-        mesh.SetMaterial (i+1, *name);
+        if (auto name = solids[i]->properties.name)
+          mesh.SetMaterial (i+1, *name);
   }
   
   shared_ptr<NetgenGeometry> GeometryRegisterArray :: LoadFromMeshFile (istream & ist) const

libsrc/meshing/boundarylayer.cpp

@@ -194,7 +194,8 @@ namespace netgen
 
   ArrayMem<Point<3>, 4> BoundaryLayerTool :: GetMappedFace( SurfaceElementIndex sei, int face )
   {
-      if(face == -1) return GetMappedFace(sei);
+      if(face == -1) return GetFace(sei);
+      if(face == -2) return GetMappedFace(sei);
       const auto & sel = mesh[sei];
       auto np = sel.GetNP();
       auto pi0 = sel[face % np];
@@ -210,23 +211,25 @@ namespace netgen
   Vec<3> BoundaryLayerTool :: getEdgeTangent(PointIndex pi, int edgenr)
   {
       Vec<3> tangent = 0.0;
+      ArrayMem<PointIndex,2> pts;
       for(auto segi : topo.GetVertexSegments(pi))
       {
           auto & seg = mesh[segi];
-          if(seg.edgenr != edgenr)
+          if(seg.edgenr != edgenr+1)
               continue;
           PointIndex other = seg[0]+seg[1]-pi;
-          tangent += mesh[other] - mesh[pi];
+          if(!pts.Contains(other))
+            pts.Append(other);
       }
+      if(pts.Size() != 2)
+        throw Exception("Something went wrong in getEdgeTangent!");
+      tangent = mesh[pts[1]] - mesh[pts[0]];
       return tangent.Normalize();
   }
 
   void BoundaryLayerTool :: LimitGrowthVectorLengths()
   {
     static Timer tall("BoundaryLayerTool::LimitGrowthVectorLengths"); RegionTimer rtall(tall);
-    height = 0.0;
-    for (auto h : params.heights)
-       height += h;
 
     limits.SetSize(np);
     limits = 1.0;
@@ -320,7 +323,7 @@ namespace netgen
                 const auto & sel = mesh[sei];
                 if(sel.PNums().Contains(pi))
                     return false;
-                auto face = GetFace(sei);
+                auto face = GetMappedFace(sei, -2);
                 double lam_ = 999;
                 bool is_bl_sel = params.surfid.Contains(sel.GetIndex());
 
@@ -433,7 +436,6 @@ namespace netgen
                   }
               }
 
-              // seg.edgenr = topo.GetEdge(segi)+1;
               segments.Append(seg);
           }
       }
@@ -464,6 +466,8 @@ namespace netgen
     auto np = mesh.GetNP();
     BitArray is_point_on_bl_surface(np+1);
     is_point_on_bl_surface.Clear();
+    BitArray is_point_on_other_surface(np+1);
+    is_point_on_other_surface.Clear();
     Array<Vec<3>, PointIndex> normals(np);
     for(auto pi : Range(growthvectors))
         normals[pi] = growthvectors[pi];
@@ -474,20 +478,33 @@ namespace netgen
           {
             auto facei = mesh[sei].GetIndex();
             if(facei < nfd_old && !params.surfid.Contains(facei))
-                continue;
+              {
-            for(auto pi : mesh[sei].PNums())
+                for(auto pi : mesh[sei].PNums())
-                if(mesh[pi].Type() == SURFACEPOINT)
+                  if(mesh[pi].Type() == SURFACEPOINT)
+                    is_point_on_other_surface.SetBitAtomic(pi);
+              }
+            else
+              {
+                for(auto pi : mesh[sei].PNums())
+                  if(mesh[pi].Type() == SURFACEPOINT)
                     is_point_on_bl_surface.SetBitAtomic(pi);
+              }
           }
       });
 
     Array<PointIndex> points;
     for(PointIndex pi : mesh.Points().Range())
+      {
         if(is_point_on_bl_surface[pi])
         {
             points.Append(pi);
             growthvectors[pi] = 0.0;
         }
+        if(is_point_on_other_surface[pi])
+          {
+            points.Append(pi);
+          }
+      }
 
     // smooth tangential part of growth vectors from edges to surface elements
     RegionTimer rtsmooth(tsmooth);
@@ -511,7 +528,10 @@ namespace netgen
                         auto gw_other = growthvectors[pi1];
                         auto normal_other = getNormal(mesh[sei]);
                         auto tangent_part = gw_other - (gw_other*normal_other)*normal_other;
-                        new_gw += tangent_part;
+                        if(is_point_on_bl_surface[pi])
+                          new_gw += tangent_part;
+                        else
+                          new_gw += gw_other;
                     }
             }
 
@@ -533,6 +553,10 @@ namespace netgen
     //for(auto & seg : mesh.LineSegments())
         //seg.edgenr = seg.epgeominfo[1].edgenr;
 
+    height = 0.0;
+    for (auto h : params.heights)
+      height += h;
+
     max_edge_nr = -1;
     for(const auto& seg : mesh.LineSegments())
       if(seg.edgenr > max_edge_nr)
@@ -777,7 +801,7 @@ namespace netgen
     // interpolate tangential component of growth vector along edge
     for(auto edgenr : Range(max_edge_nr))
       {
-        if(!is_edge_moved[edgenr+1]) continue;
+        // if(!is_edge_moved[edgenr+1]) continue;
 
         // build sorted list of edge
         Array<PointIndex> points;
@@ -785,6 +809,9 @@ namespace netgen
         double edge_len = 0.;
         auto is_end_point = [&] (PointIndex pi)
         {
+          // if(mesh[pi].Type() == FIXEDPOINT)
+          //   return true;
+          // return false;
             auto segs = topo.GetVertexSegments(pi);
             auto first_edgenr = mesh[segs[0]].edgenr;
             for(auto segi : segs)
@@ -792,17 +819,31 @@ namespace netgen
                     return true;
             return false;
         };
+
+        bool any_grows = false;
+        
         for(const auto& seg : segments)
           {
-            if(seg.edgenr-1 == edgenr && is_end_point(seg[0]))
+            if(seg.edgenr-1 == edgenr)
               {
-                points.Append(seg[0]);
+                if(growthvectors[seg[0]].Length2() != 0 ||
-                points.Append(seg[1]);
+                   growthvectors[seg[1]].Length2() != 0)
-                edge_len += (mesh[seg[1]] - mesh[seg[0]]).Length();
+                  any_grows = true;
-                break;
+                if(points.Size() == 0 && is_end_point(seg[0]))
+                  {
+                    points.Append(seg[0]);
+                    points.Append(seg[1]);
+                    edge_len += (mesh[seg[1]] - mesh[seg[0]]).Length();
+                  }
               }
           }
 
+        if(!any_grows)
+          continue;
+
+        if(!points.Size())
+          throw Exception("Could not find startpoint for edge " + ToString(edgenr));
+
         while(true)
           {
             bool point_found = false;
@@ -831,17 +872,28 @@ namespace netgen
               break;
             if(!point_found)
               {
-                cout << "points = " << points << endl;
                 throw Exception(string("Could not find connected list of line segments for edge ") + edgenr);
               }
           }
 
+        if(growthvectors[points[0]].Length2() == 0 &&
+           growthvectors[points.Last()].Length2() == 0)
+          continue;
+
         // tangential part of growth vectors
         auto t1 = (mesh[points[1]]-mesh[points[0]]).Normalize();
         auto gt1 = growthvectors[points[0]] * t1 * t1;
         auto t2 = (mesh[points.Last()]-mesh[points[points.Size()-2]]).Normalize();
         auto gt2 = growthvectors[points.Last()] * t2 * t2;
 
+        if(!is_edge_moved[edgenr+1])
+          {
+            if(growthvectors[points[0]] * (mesh[points[1]] - mesh[points[0]]) < 0)
+              gt1 = 0.;
+            if(growthvectors[points.Last()] * (mesh[points[points.Size()-2]] - mesh[points.Last()]) < 0)
+              gt2 = 0.;
+          }
+
         double len = 0.;
         for(size_t i = 1; i < points.Size()-1; i++)
           {
@@ -1262,7 +1314,9 @@ namespace netgen
 
       auto in_surface_direction = ProjectGrowthVectorsOnSurface();
       InterpolateGrowthVectors();
-      LimitGrowthVectorLengths();
+
+      if(params.limit_growth_vectors)
+        LimitGrowthVectorLengths();
       FixVolumeElements();
       InsertNewElements(segmap, in_surface_direction);
       SetDomInOut();
@@ -1339,20 +1393,24 @@ namespace netgen
      compress = PointIndex::INVALID;
 
      PointIndex cnt = PointIndex::BASE;
-     for(const auto & seg : mesh.LineSegments())
-        if (seg.si == domain)
-           for (auto pi : {seg[0], seg[1]})
-              if (compress[pi]==PointIndex{PointIndex::INVALID})
-              {
-                 meshing.AddPoint(mesh[pi], pi);
-                 compress[pi] = cnt++;
-              }
 
+     auto p2sel = mesh.CreatePoint2SurfaceElementTable();
+     Array<Segment> domain_segments;
      PointGeomInfo gi;
      gi.trignum = domain;
-     for(const auto & seg : mesh.LineSegments())
+     for(auto seg : mesh.LineSegments())
-        if (seg.si == domain)
+     {
-           meshing.AddBoundaryElement (compress[seg[0]], compress[seg[1]], gi, gi);
+         for (auto pi : {seg[0], seg[1]})
+            if (compress[pi]==PointIndex{PointIndex::INVALID})
+            {
+               meshing.AddPoint(mesh[pi], pi);
+               compress[pi] = cnt++;
+            }
+         if(seg.domin == domain)
+             meshing.AddBoundaryElement (compress[seg[0]], compress[seg[1]], gi, gi);
+         if(seg.domout == domain)
+             meshing.AddBoundaryElement (compress[seg[1]], compress[seg[0]], gi, gi);
+     }
 
      auto oldnf = mesh.GetNSE();
      auto res = meshing.GenerateMesh (mesh, mp, mp.maxh, domain);
@@ -1395,17 +1453,20 @@ namespace netgen
      }
 
      mesh.Compress();
+     mesh.CalcSurfacesOfNode();
      mesh.OrderElements();
      mesh.SetNextMajorTimeStamp();
-
   }
 
   int GenerateBoundaryLayer2 (Mesh & mesh, int domain, const Array<double> & thicknesses, bool should_make_new_domain, const Array<int> & boundaries)
   {
+     mesh.GetTopology().SetBuildVertex2Element(true);
+     mesh.UpdateTopology();
+     const auto & line_segments = mesh.LineSegments();
      SegmentIndex first_new_seg = mesh.LineSegments().Range().Next();
 
      int np = mesh.GetNP();
-     int nseg = mesh.GetNSeg();
+     int nseg = line_segments.Size();
      int ne = mesh.GetNSE();
      mesh.UpdateTopology();
 
@@ -1427,23 +1488,10 @@ namespace netgen
 
      auto & meshtopo = mesh.GetTopology();
 
-     Array<SurfaceElementIndex, SegmentIndex> seg2surfel(mesh.GetNSeg());
-     seg2surfel = -1;
-     NgArray<SurfaceElementIndex> temp_els;
-     for(auto si : Range(mesh.LineSegments()))
-     {
-        meshtopo.GetSegmentSurfaceElements ( si+1, temp_els );
-        // NgArray<int> surfeledges;
-        // meshtopo.GetSurfaceElementEdges(si+1, surfeledges);
-        for(auto seli : temp_els)
-           if(mesh[seli].GetIndex() == mesh[si].si)
-              seg2surfel[si] = seli;
-     }
-
      Array<SegmentIndex> segments;
 
     // surface index map
-    Array<int> si_map(mesh.GetNFD()+1);
+    Array<int> si_map(mesh.GetNFD()+2);
     si_map = -1;
 
     int fd_old = mesh.GetNFD();
@@ -1451,12 +1499,14 @@ namespace netgen
     int max_edge_nr = -1;
     int max_domain = -1;
 
-    for(const auto& seg : mesh.LineSegments())
+    for(const auto& seg : line_segments)
     {
       if(seg.epgeominfo[0].edgenr > max_edge_nr)
         max_edge_nr = seg.epgeominfo[0].edgenr;
-      if(seg.si > max_domain)
+      if(seg.domin > max_domain)
-         max_domain = seg.si;
+         max_domain = seg.domin;
+      if(seg.domout > max_domain)
+         max_domain = seg.domout;
     }
 
     int new_domain = max_domain+1;
@@ -1472,95 +1522,43 @@ namespace netgen
        for(auto edgenr : boundaries)
           active_boundaries.SetBit(edgenr);
 
-    for(auto segi : Range(mesh.LineSegments()))
+    for(auto segi : Range(line_segments))
     {
-       const auto seg = mesh[segi];
+       const auto seg = line_segments[segi];
-       if(active_boundaries.Test(seg.epgeominfo[0].edgenr) && seg.si==domain)
+       if(active_boundaries.Test(seg.epgeominfo[0].edgenr) && (seg.domin==domain || seg.domout==domain))
           active_segments.SetBit(segi);
     }
 
-    for(auto segi : Range(mesh.LineSegments()))
     {
-        const auto& seg = mesh[segi];
-        auto si = seg.si;
-
-        if(si_map[si]!=-1)
-           continue;
-
-        if(!active_segments.Test(segi))
-           continue;
-
         FaceDescriptor new_fd(0, 0, 0, -1);
         new_fd.SetBCProperty(new_domain);
         int new_fd_index = mesh.AddFaceDescriptor(new_fd);
-        si_map[si] = new_domain;
         if(should_make_new_domain)
-           mesh.SetBCName(new_domain-1, "mapped_" + mesh.GetBCName(si-1));
+           mesh.SetBCName(new_domain-1, "mapped_" + mesh.GetBCName(domain-1));
     }
 
-    for(auto si : Range(mesh.LineSegments()))
+    for(auto segi : Range(line_segments))
       {
-        if(segs_done[si]) continue;
+        if(segs_done[segi]) continue;
-        segs_done.SetBit(si);
+        segs_done.SetBit(segi);
-        const auto& segi = mesh[si];
+        const auto& seg = line_segments[segi];
-        if(si_map[segi.si] == -1) continue;
+        if(seg.domin != domain && seg.domout != domain) continue;
-        if(!active_boundaries.Test(segi.epgeominfo[0].edgenr))
+        if(!active_boundaries.Test(seg.epgeominfo[0].edgenr))
            continue;
-        moved_segs.Append(si);
+        moved_segs.Append(segi);
       }
 
      // calculate growth vectors (average normal vectors of adjacent segments at each point)
      for (auto si : moved_segs)
      {
-       auto & seg = mesh[si];
+       auto & seg = line_segments[si];
-
-       meshtopo.GetSegmentSurfaceElements ( si+1, temp_els );
-       ArrayMem<int, 10> seg_domains;
-
-       temp_els.SetSize(0);
-       if(seg2surfel[si]!=-1)
-          temp_els.Append(seg2surfel[si]);
-
-       int n_temp_els = temp_els.Size();
-       if(n_temp_els==0)
-          continue;
-
-       int dom0 = mesh[temp_els[0]].GetIndex();
-       int dom1 = n_temp_els==2 ? mesh[temp_els[1]].GetIndex() : 0;
-
-       bool in_dom0 = dom0 == domain;
-       bool in_dom1 = dom1 == domain;
-
-       if(!in_dom0 && !in_dom1)
-          continue;
-
-       int side = in_dom0 ? 0 : 1;
-
-       auto & sel = mesh[ temp_els[side] ];
-
-       int domain = sel.GetIndex();
-       Vec<3> pcenter = 0.0;
-       for(auto i : IntRange(sel.GetNP()))
-       {
-          for(auto d : IntRange(3))
-             pcenter[d] += mesh[sel[i]][d];
-       }
-       pcenter = 1.0/sel.GetNP() * pcenter;
 
        auto n = mesh[seg[1]] - mesh[seg[0]];
        n = {-n[1], n[0], 0};
        n.Normalize();
 
-       Vec<3> p0{mesh[seg[0]]};
+       if(seg.domout == domain)
-       Vec<3> p1{mesh[seg[0]]};
+           n = -n;
-
-
-       auto v = pcenter -0.5*(p0+p1);
-
-       const auto Dot = [](Vec<3> a, Vec<3> b)
-       { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; };
-       if(Dot(n, v)<0)
-          n = -1*n;
 
        AddDirection(growthvectors[seg[0]], n);
        AddDirection(growthvectors[seg[1]], n);
@@ -1573,7 +1571,7 @@ namespace netgen
 
      for(auto si : moved_segs)
      {
-        auto current_seg = mesh[si];
+        auto current_seg = line_segments[si];
         auto current_si = si;
 
         auto first = current_seg[0];
@@ -1711,8 +1709,9 @@ namespace netgen
         const auto & seg0 = mesh[segi0];
         const auto & seg1 = mesh[segi1];
 
-        if(seg0.si != seg1.si)
+        if( (seg0.domin != domain && seg0.domout != domain) ||
-           return;
+            (seg1.domin != domain && seg1.domout != domain) )
+            return;
 
         if(segi0 == segi1)
            return;
@@ -1825,7 +1824,7 @@ namespace netgen
      {
         auto p2 = [](Point<3> p) { return Point<2>{p[0], p[1]}; };
 
-        auto seg = mesh[segi];
+        auto seg = line_segments[segi];
         double alpha,beta;
         intersect( p2(mesh[seg[0]]), p2(mesh[seg[0]]+total_thickness*growthvectors[seg[0]]), p2(mesh[seg[1]]), p2(mesh[seg[1]]+total_thickness*growthvectors[seg[1]]), alpha, beta );
 
@@ -1863,13 +1862,13 @@ namespace netgen
      // insert new elements ( and move old ones )
      for(auto si : moved_segs)
      {
-        auto seg = mesh[si];
+        auto seg = line_segments[si];
 
         bool swap = false;
         auto & pm0 = mapto[seg[0]];
         auto & pm1 = mapto[seg[1]];
 
-        auto newindex = si_map[seg.si];
+        auto newindex = si_map[domain];
 
         Segment s = seg;
         s.geominfo[0] = {};
@@ -1884,10 +1883,6 @@ namespace netgen
         s.si = seg.si;
         mesh.AddSegment(s);
 
-        Swap(s[0], s[1]);
-        s.si =  newindex;
-        mesh.AddSegment(s);
-
         for ( auto i : Range(thicknesses))
         {
            PointIndex pi0, pi1, pi2, pi3;
@@ -1925,14 +1920,14 @@ namespace netgen
            newel[1] = pi1;
            newel[2] = pi2;
            newel[3] = pi3;
-           newel.SetIndex(si_map[seg.si]);
+           newel.SetIndex(new_domain);
            newel.GeomInfo() = PointGeomInfo{};
 
-//            if(swap)
+            if(swap)
-//            {
+            {
-//               Swap(newel[0], newel[1]);
+               Swap(newel[0], newel[1]);
-//               Swap(newel[2], newel[3]);
+               Swap(newel[2], newel[3]);
-//            }
+            }
 
            for(auto i : Range(4))
            {
@@ -1943,7 +1938,10 @@ namespace netgen
 
         }
         // segment now adjacent to new 2d-domain!
-        mesh[si].si = si_map[seg.si];
+        if(line_segments[si].domin == domain)
+            line_segments[si].domin = new_domain;
+        if(line_segments[si].domout == domain)
+            line_segments[si].domout = new_domain;
      }
 
      for(auto pi : Range(mapto))
@@ -1988,7 +1986,12 @@ namespace netgen
         }
 
         for(auto segi : moved_segs)
-           mesh[segi].si = domain;
+        {
+            if(mesh[segi].domin == new_domain)
+                mesh[segi].domin = domain;
+            if(mesh[segi].domout == new_domain)
+                mesh[segi].domout = domain;
+        }
 
         mesh.Compress();
         mesh.CalcSurfacesOfNode();

libsrc/meshing/boundarylayer.hpp

@@ -18,6 +18,7 @@ public:
   BitArray domains;
   bool outside = false; // set the boundary layer on the outside
   bool grow_edges = false;
+  bool limit_growth_vectors = true;
   Array<size_t> project_boundaries;
 };
 

libsrc/meshing/curvedelems.cpp

@@ -1632,7 +1632,7 @@ namespace netgen
     if (mesh.coarsemesh)
       {
 	const HPRefElement & hpref_el =
-	  (*mesh.hpelements) [mesh[elnr].hp_elnr];
+	  (*mesh.hpelements) [mesh[elnr].GetHpElnr()];
 	
 	return mesh.coarsemesh->GetCurvedElements().IsSurfaceElementCurved (hpref_el.coarse_elnr);
       }
@@ -1679,7 +1679,7 @@ namespace netgen
     if (mesh.coarsemesh)
       {
 	const HPRefElement & hpref_el =
-	  (*mesh.hpelements) [mesh[elnr].hp_elnr];
+	  (*mesh.hpelements) [mesh[elnr].GetHpElnr()];
 	
 	// xi umrechnen
 	double lami[4];
@@ -2428,7 +2428,7 @@ namespace netgen
     if (mesh.coarsemesh)
       {
 	const HPRefElement & hpref_el =
-	  (*mesh.hpelements) [mesh[elnr].hp_elnr];
+	  (*mesh.hpelements) [mesh[elnr].GetHpElnr()];
 	
 	return mesh.coarsemesh->GetCurvedElements().IsElementCurved (hpref_el.coarse_elnr);
       }
@@ -2480,7 +2480,7 @@ namespace netgen
     if (mesh.coarsemesh)
       {
 	const HPRefElement & hpref_el =
-	  (*mesh.hpelements) [mesh[elnr].hp_elnr];
+	  (*mesh.hpelements) [mesh[elnr].GetHpElnr()];
 	
 	return mesh.coarsemesh->GetCurvedElements().IsElementHighOrder (hpref_el.coarse_elnr);
       }
@@ -2519,7 +2519,7 @@ namespace netgen
     if (mesh.coarsemesh)
       {
 	const HPRefElement & hpref_el =
-	  (*mesh.hpelements) [mesh[elnr].hp_elnr];
+	  (*mesh.hpelements) [mesh[elnr].GetHpElnr()];
 	  
 	// xi umrechnen
 	double lami[8];
@@ -4065,7 +4065,7 @@ namespace netgen
     if (mesh.coarsemesh)
       {
 	const HPRefElement & hpref_el =
-	  (*mesh.hpelements) [mesh[elnr].hp_elnr];
+	  (*mesh.hpelements) [mesh[elnr].GetHpElnr()];
 	
 	// xi umrechnen
 	T lami[4];
@@ -4529,7 +4529,7 @@ namespace netgen
     if (mesh.coarsemesh)
       {
 	const HPRefElement & hpref_el =
-	  (*mesh.hpelements) [mesh[elnr].hp_elnr];
+	  (*mesh.hpelements) [mesh[elnr].GetHpElnr()];
 	
 	// xi umrechnen
 	T lami[8];

libsrc/meshing/delaunay2d.cpp

@@ -668,6 +668,11 @@ namespace netgen
     for (auto pi : points_range)
       dmesh.AddPoint(pi);
 
+    auto first_new_point = points_range.Next();
+    tempmesh.AddPoint(P3(temp_points[first_new_point]));
+    tempmesh.AddPoint(P3(temp_points[first_new_point+1]));
+    tempmesh.AddPoint(P3(temp_points[first_new_point+2]));
+
     timer_addpoints.Stop();
 
     static Timer taddseg("addseg");

libsrc/meshing/fieldlines.cpp

@@ -77,7 +77,7 @@ namespace netgen
   {
     bool finished = false;
     
-    if(stepcount <= steps)
+    if(stepcount <= steps && stepcount>0)
       {
 	t = startt + c[stepcount-1]*h;
 	val = startval;

libsrc/meshing/hprefinement.cpp

@@ -1403,7 +1403,7 @@ namespace netgen
 		    Element2d el(hpel.np); 
 		    for(int j=0;j<hpel.np;j++) 
 		      el.PNum(j+1) = hpel.pnums[j]; 
-		    el.hp_elnr = i; 
+		    el.SetHpElnr(i);
 		    el.SetIndex(hpel.index);
 		    if(setorders)
 		      el.SetOrder(act_ref+1,act_ref+1,0); 
@@ -1421,7 +1421,7 @@ namespace netgen
 		    for(int j=0;j<hpel.np;j++) 
 		      el.PNum(j+1) = hpel.pnums[j]; 
 		    el.SetIndex(hpel.index); 
-		    el.hp_elnr = i; 
+		    el.SetHpElnr(i);
 		    if(setorders)
 		      el.SetOrder(act_ref+1,act_ref+1,act_ref+1);
                     if((*mesh.coarsemesh)[ElementIndex{hpel.coarse_elnr}].IsCurved())
@@ -1451,7 +1451,7 @@ namespace netgen
 	for(ElementIndex i=0;i<mesh.GetNE(); i++) 
 	  { 
 	    // Element el = mesh[i] ;
-	    HPRefElement & hpel = hpelements[mesh[i].hp_elnr];
+	    HPRefElement & hpel = hpelements[mesh[i].GetHpElnr()];
 	    const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (mesh[i].GetType());
 	    double dist[3] = {0,0,0}; 
 	    int ord_dir[3] = {0,0,0}; 
@@ -1523,7 +1523,7 @@ namespace netgen
 	for(SurfaceElementIndex i=0;i<mesh.GetNSE(); i++) 
 	  { 
 	    // Element2d el = mesh[i] ;
-	    HPRefElement & hpel = hpelements[mesh[i].hp_elnr];
+	    HPRefElement & hpel = hpelements[mesh[i].GetHpElnr()];
 	    const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (mesh[i].GetType());
 	    double dist[3] = {0,0,0}; 
 	    int ord_dir[3] = {0,0,0}; 

libsrc/meshing/improve2.cpp

@@ -339,6 +339,9 @@ namespace netgen
             if (swapped[t1])
               return;
 
+            if(mesh[t1].GetNP() != 3)
+              return;
+
             if (multithread.terminate)
               throw NgException ("Meshing stopped");
 

libsrc/meshing/improve2.hpp

@@ -19,7 +19,7 @@ inline void AppendEdges( const Element & elem, PointIndex pi, Array<std::tuple<P
   { { 0, 1 }, { 0, 2 }, { 0, 3 },
       { 1, 2 }, { 1, 3 }, { 2, 3 } };
 
-  if(elem.flags.fixed)
+  if(elem.Flags().fixed)
       return;
   for (int j = 0; j < 6; j++)
   {

libsrc/meshing/improve3.cpp

@@ -41,7 +41,7 @@ static ArrayMem<Element, 3> SplitElement (Element old, PointIndex pi0, PointInde
   ArrayMem<Element, 3> new_elements;
   // split element by cutting edge pi0,pi1 at pinew
   auto np = old.GetNP();
-  old.flags.illegal_valid = 0;
+  old.Flags().illegal_valid = 0;
   if(np == 4)
   {
     // Split tet into two tets
@@ -232,7 +232,7 @@ double MeshOptimize3d :: CombineImproveEdge (Mesh & mesh,
                   break;
               }
 
-          elem.flags.illegal_valid = 0;
+          elem.Flags().illegal_valid = 0;
           if (!mesh.LegalTet(elem))
               badness_new += 1e4;
       }
@@ -255,7 +255,7 @@ double MeshOptimize3d :: CombineImproveEdge (Mesh & mesh,
               if (elem[l] == pi1)
                   elem[l] = pi0;
 
-          elem.flags.illegal_valid = 0;
+          elem.Flags().illegal_valid = 0;
           if (!mesh.LegalTet (elem))
               (*testout) << "illegal tet " << ei << endl;
       }
@@ -265,7 +265,7 @@ double MeshOptimize3d :: CombineImproveEdge (Mesh & mesh,
 
       for (auto ei : has_both_points)
       {
-          mesh[ei].flags.illegal_valid = 0;
+          mesh[ei].Flags().illegal_valid = 0;
           mesh[ei].Delete();
       }
   }
@@ -505,9 +505,9 @@ double MeshOptimize3d :: SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table
       Element newel1 = oldel;
       Element newel2 = oldel;
 
-      oldel.flags.illegal_valid = 0;
+      oldel.Flags().illegal_valid = 0;
-      newel1.flags.illegal_valid = 0;
+      newel1.Flags().illegal_valid = 0;
-      newel2.flags.illegal_valid = 0;
+      newel2.Flags().illegal_valid = 0;
 
       for (int l = 0; l < 4; l++)
         {
@@ -536,11 +536,11 @@ double MeshOptimize3d :: SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table
           Element newel1 = oldel;
           Element newel2 = oldel;
 
-          oldel.flags.illegal_valid = 0;
+          oldel.Flags().illegal_valid = 0;
           oldel.Delete();
 
-          newel1.flags.illegal_valid = 0;
+          newel1.Flags().illegal_valid = 0;
-          newel2.flags.illegal_valid = 0;
+          newel2.Flags().illegal_valid = 0;
 
           for (int l = 0; l < 4; l++)
             {
@@ -799,9 +799,9 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
           CalcBad (mesh.Points(), el32, 0) +
           CalcBad (mesh.Points(), el33, 0);
 
-      el31.flags.illegal_valid = 0;
+      el31.Flags().illegal_valid = 0;
-      el32.flags.illegal_valid = 0;
+      el32.Flags().illegal_valid = 0;
-      el33.flags.illegal_valid = 0;
+      el33.Flags().illegal_valid = 0;
 
       if (!mesh.LegalTet(el31) ||
               !mesh.LegalTet(el32) ||
@@ -823,8 +823,8 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
       bad2 = CalcBad (mesh.Points(), el21, 0) +
           CalcBad (mesh.Points(), el22, 0);
 
-      el21.flags.illegal_valid = 0;
+      el21.Flags().illegal_valid = 0;
-      el22.flags.illegal_valid = 0;
+      el22.Flags().illegal_valid = 0;
 
       if (!mesh.LegalTet(el21) ||
               !mesh.LegalTet(el22))
@@ -866,8 +866,8 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
           mesh[hasbothpoints[1]].Delete();
           mesh[hasbothpoints[2]].Delete();
 
-          el21.flags.illegal_valid = 0;
+          el21.Flags().illegal_valid = 0;
-          el22.flags.illegal_valid = 0;
+          el22.Flags().illegal_valid = 0;
           mesh.AddVolumeElement(el21);
           mesh.AddVolumeElement(el22);
         }
@@ -942,10 +942,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
           CalcBad (mesh.Points(), el4, 0);
 
 
-      el1.flags.illegal_valid = 0;
+      el1.Flags().illegal_valid = 0;
-      el2.flags.illegal_valid = 0;
+      el2.Flags().illegal_valid = 0;
-      el3.flags.illegal_valid = 0;
+      el3.Flags().illegal_valid = 0;
-      el4.flags.illegal_valid = 0;
+      el4.Flags().illegal_valid = 0;
 
 
       if (goal != OPT_CONFORM)
@@ -978,10 +978,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
           CalcBad (mesh.Points(), el3, 0) +
           CalcBad (mesh.Points(), el4, 0);
 
-      el1.flags.illegal_valid = 0;
+      el1.Flags().illegal_valid = 0;
-      el2.flags.illegal_valid = 0;
+      el2.Flags().illegal_valid = 0;
-      el3.flags.illegal_valid = 0;
+      el3.Flags().illegal_valid = 0;
-      el4.flags.illegal_valid = 0;
+      el4.Flags().illegal_valid = 0;
 
       if (goal != OPT_CONFORM)
         {
@@ -1014,10 +1014,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
           CalcBad (mesh.Points(), el3b, 0) +
           CalcBad (mesh.Points(), el4b, 0);
 
-      el1b.flags.illegal_valid = 0;
+      el1b.Flags().illegal_valid = 0;
-      el2b.flags.illegal_valid = 0;
+      el2b.Flags().illegal_valid = 0;
-      el3b.flags.illegal_valid = 0;
+      el3b.Flags().illegal_valid = 0;
-      el4b.flags.illegal_valid = 0;
+      el4b.Flags().illegal_valid = 0;
 
       if (goal != OPT_CONFORM)
         {
@@ -1054,10 +1054,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
           for (auto i : IntRange(4))
               mesh[hasbothpoints[i]].Delete();
 
-          el1.flags.illegal_valid = 0;
+          el1.Flags().illegal_valid = 0;
-          el2.flags.illegal_valid = 0;
+          el2.Flags().illegal_valid = 0;
-          el3.flags.illegal_valid = 0;
+          el3.Flags().illegal_valid = 0;
-          el4.flags.illegal_valid = 0;
+          el4.Flags().illegal_valid = 0;
           mesh.AddVolumeElement (el1);
           mesh.AddVolumeElement (el2);
           mesh.AddVolumeElement (el3);
@@ -1068,10 +1068,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
           for (auto i : IntRange(4))
               mesh[hasbothpoints[i]].Delete();
 
-          el1b.flags.illegal_valid = 0;
+          el1b.Flags().illegal_valid = 0;
-          el2b.flags.illegal_valid = 0;
+          el2b.Flags().illegal_valid = 0;
-          el3b.flags.illegal_valid = 0;
+          el3b.Flags().illegal_valid = 0;
-          el4b.flags.illegal_valid = 0;
+          el4b.Flags().illegal_valid = 0;
           mesh.AddVolumeElement (el1b);
           mesh.AddVolumeElement (el2b);
           mesh.AddVolumeElement (el3b);
@@ -1174,7 +1174,7 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
               hel[3] = pi2;
 
               bad2 += CalcBad (mesh.Points(), hel, 0);
-              hel.flags.illegal_valid = 0;
+              hel.Flags().illegal_valid = 0;
               if (!mesh.LegalTet(hel)) bad2 += 1e4;
 
               hel[2] = suroundpts[k % nsuround];
@@ -1183,7 +1183,7 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
 
               bad2 += CalcBad (mesh.Points(), hel, 0);
 
-              hel.flags.illegal_valid = 0;
+              hel.Flags().illegal_valid = 0;
               if (!mesh.LegalTet(hel)) bad2 += 1e4;
             }
           // (*testout) << "bad2," << l << " = " << bad2 << endl;
@@ -1253,7 +1253,7 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
               hel[1] = suroundpts[k % nsuround];
               hel[2] = suroundpts[(k+1) % nsuround];
               hel[3] = pi2;
-              hel.flags.illegal_valid = 0;
+              hel.Flags().illegal_valid = 0;
 
               /*
                  (*testout) << nsuround << "-swap, new el,top = "
@@ -1309,7 +1309,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
         for (ElementIndex eli : myrange)
           {
             const auto & el = mesh[eli];
-            if(el.flags.fixed)
+            if(el.Flags().fixed)
               continue;
 
             for (auto pi : el.PNums())
@@ -2412,9 +2412,9 @@ double MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementI
                 CalcBad (mesh.Points(), el33, 0);
 
 
-              el31.flags.illegal_valid = 0;
+              el31.Flags().illegal_valid = 0;
-              el32.flags.illegal_valid = 0;
+              el32.Flags().illegal_valid = 0;
-              el33.flags.illegal_valid = 0;
+              el33.Flags().illegal_valid = 0;
 
               if (!mesh.LegalTet(el31) ||
                   !mesh.LegalTet(el32) ||
@@ -2433,9 +2433,9 @@ double MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementI
 
               if (d_badness<0.0)
               {
-                  el31.flags.illegal_valid = 0;
+                  el31.Flags().illegal_valid = 0;
-                  el32.flags.illegal_valid = 0;
+                  el32.Flags().illegal_valid = 0;
-                  el33.flags.illegal_valid = 0;
+                  el33.Flags().illegal_valid = 0;
 
                   mesh[eli1].Delete();
                   mesh[eli2].Delete();
@@ -2655,7 +2655,7 @@ double MeshOptimize3d :: SplitImprove2Element (Mesh & mesh,
   if(badness_after<badness_before)
   {
     PointIndex pinew = mesh.AddPoint (center);
-    el.flags.illegal_valid = 0;
+    el.Flags().illegal_valid = 0;
     el.Delete();
 
     for (auto ei1 : has_both_points0)

libsrc/meshing/meshclass.cpp

@@ -599,9 +599,9 @@ namespace netgen
       {
         volelements.Append (el);
       }
-    volelements.Last().flags.illegal_valid = 0;
+    volelements.Last().Flags().illegal_valid = 0;
-    volelements.Last().flags.fixed = 0;
+    volelements.Last().Flags().fixed = 0;
-    volelements.Last().flags.deleted = 0;
+    volelements.Last().Flags().deleted = 0;
 
     // while (volelements.Size() > eltyps.Size())
     // eltyps.Append (FREEELEMENT);
@@ -622,9 +622,9 @@ namespace netgen
     */
 
     volelements[ei]  = el;
-    volelements[ei].flags.illegal_valid = 0;
+    volelements[ei].Flags().illegal_valid = 0;
-    volelements[ei].flags.fixed = 0;
+    volelements[ei].Flags().fixed = 0;
-    volelements[ei].flags.deleted = 0;
+    volelements[ei].Flags().deleted = 0;
   }
 
 
@@ -1660,10 +1660,6 @@ namespace netgen
 	clusters -> Update();
       }
 
-    auto geo = geometryregister.LoadFromMeshFile (infile);
-    if(geo)
-       geometry = geo;
-
     SetNextMajorTimeStamp();
     //  PrintMemInfo (cout);
   }
@@ -3245,7 +3241,7 @@ namespace netgen
           if (dist[el[j]] < elmin)
             elmin = dist[el[j]];
 
-        el.flags.fixed = elmin > layers;
+        el.Flags().fixed = elmin > layers;
         // eltyps.Elem(i) = (elmin <= layers) ? 
         // FREEELEMENT : FIXEDELEMENT;
         if (elmin <= layers)
@@ -4388,7 +4384,7 @@ namespace netgen
     for (i = 1; i <= ne; i++)
       {
         Element & el = (Element&) VolumeElement(i);
-        el.flags.badel = 0;
+        el.Flags().badel = 0;
         int nip = el.GetNIP();
         for (j = 1; j <= nip; j++)
           {
@@ -4397,7 +4393,7 @@ namespace netgen
             if (det > 0)
               {
                 PrintError ("Element ", i , " has wrong orientation");
-                el.flags.badel = 1;
+                el.Flags().badel = 1;
               }
           }
       }
@@ -6476,7 +6472,7 @@ namespace netgen
 
         if (el.GetType() != TET)
           {
-            VolumeElement(i).flags.badel = 0;
+            VolumeElement(i).Flags().badel = 0;
             continue;
           }
 
@@ -6558,7 +6554,7 @@ namespace netgen
           }
 
 
-        VolumeElement(i).flags.badel = badel;
+        VolumeElement(i).Flags().badel = badel;
         if (badel) badtets++;
       }
 

libsrc/meshing/meshclass.hpp

@@ -365,7 +365,7 @@ namespace netgen
     Element & operator[] (ElementIndex ei) { return volelements[ei]; }
 
     ELEMENTTYPE ElementType (ElementIndex i) const 
-    { return (volelements[i].flags.fixed) ? FIXEDELEMENT : FREEELEMENT; }
+    { return (volelements[i].Flags().fixed) ? FIXEDELEMENT : FREEELEMENT; }
 
     const auto & VolumeElements() const { return volelements; }
     auto & VolumeElements() { return volelements; }
@@ -384,7 +384,7 @@ namespace netgen
     DLL_HEADER int GetNDomains() const;
     ///
     int GetDimension() const { return dimension; }
-    void SetDimension (int dim); //  { dimension = dim; }
+    DLL_HEADER void SetDimension (int dim); //  { dimension = dim; }
 
     /// sets internal tables
     DLL_HEADER void CalcSurfacesOfNode ();
@@ -478,7 +478,7 @@ namespace netgen
         return lochfunc[0];
       return lochfunc[layer-1];
     }
-    void SetLocalH(shared_ptr<LocalH> loch, int layer=1);
+    DLL_HEADER void SetLocalH(shared_ptr<LocalH> loch, int layer=1);
 
     ///
     bool LocalHFunctionGenerated(int layer=1) const { return (lochfunc[layer-1] != NULL); }

libsrc/meshing/meshfunc.cpp

@@ -201,6 +201,7 @@ namespace netgen
           return;
 
       NgArray<int, PointIndex::BASE> map;
+      std::set<std::tuple<int,int,int>> hex_faces;
       for(auto identnr : Range(1,nmax+1))
       {
           if(identifications.GetType(identnr) != Identifications::CLOSESURFACES)
@@ -211,6 +212,15 @@ namespace netgen
 
           for(auto & sel : mesh.OpenElements())
           {
+              // For quads: check if this open element is already closed by a hex
+              // this happends when we have identifications in two directions
+              if(sel.GetNP() == 4)
+              {
+                  Element2d face = sel;
+                  face.NormalizeNumbering();
+                  if(hex_faces.count({face[0], face[1], face[2]}))
+                      continue;
+              }
               bool is_mapped = true;
               for(auto pi : sel.PNums())
                   if(!PointIndex(map[pi]).IsValid())
@@ -235,21 +245,26 @@ namespace netgen
               if(pis.size() < 2*np)
                   continue;
 
-              bool is_domout = md.domain == mesh.GetFaceDescriptor(sel.GetIndex()).DomainOut();
-
               // check if new element is inside current domain
               auto p0 = mesh[sel[0]];
-              Vec<3> n = Cross(mesh[sel[1]] - p0, mesh[sel[2]] - p0 );
+              Vec<3> n = -Cross(mesh[sel[1]] - p0, mesh[sel[2]] - p0 );
-              n = is_domout ? n : -n;
 
               if(n*(mesh[el[np]]-p0) < 0.0)
                   continue;
 
-              if(is_domout)
-                  el.Invert();
-
               el.SetIndex(md.domain);
               mesh.AddVolumeElement(el);
+              if(el.NP()==8)
+              {
+                  // remember all adjacent faces of the new hex (to skip corresponding openelements accordingly)
+                  for(auto facei : Range(1,7))
+                  {
+                      Element2d face;
+                      el.GetFace(facei, face);
+                      face.NormalizeNumbering();
+                      hex_faces.insert({face[0], face[1], face[2]});
+                  }
+              }
           }
       }
   }
@@ -571,6 +586,8 @@ namespace netgen
 
      auto md = DivideMesh(mesh3d, mp);
 
+     try
+       {
      ParallelFor( md.Range(), [&](int i)
        {
          if (mp.checkoverlappingboundary)
@@ -578,10 +595,16 @@ namespace netgen
              throw NgException ("Stop meshing since boundary mesh is overlapping");
 
          if(md[i].mesh->GetGeometry()->GetGeomType() == Mesh::GEOM_OCC)
-           FillCloseSurface( md[i] );
+            FillCloseSurface( md[i] );
          CloseOpenQuads( md[i] );
          MeshDomain(md[i]);
        }, md.Size());
+       }
+     catch(...)
+       {
+         MergeMeshes(mesh3d, md);
+         return MESHING3_GIVEUP;
+       }
 
      MergeMeshes(mesh3d, md);
 

libsrc/meshing/meshtype.cpp

@@ -178,62 +178,6 @@ namespace netgen
     */
   }    
 
-  Segment::Segment (const Segment & other)
-    : 
-    edgenr(other.edgenr),
-    singedge_left(other.singedge_left),
-    singedge_right(other.singedge_right),
-    seginfo(other.seginfo),
-    si(other.si),
-    domin(other.domin),
-    domout(other.domout),
-    tlosurf(other.tlosurf),
-    geominfo(),
-    surfnr1(other.surfnr1),
-    surfnr2(other.surfnr2),
-    epgeominfo(),
-    meshdocval(other.meshdocval),
-    is_curved(other.is_curved),
-    hp_elnr(other.hp_elnr)
-  {
-    for (int j = 0; j < 3; j++)
-      pnums[j] = other.pnums[j];
-
-    geominfo[0] = other.geominfo[0];
-    geominfo[1] = other.geominfo[1];
-    epgeominfo[0] = other.epgeominfo[0];
-    epgeominfo[1] = other.epgeominfo[1];
-  }
-
-  Segment& Segment::operator=(const Segment & other)
-  {
-    if (&other != this)
-      {
-	pnums[0] = other[0];
-	pnums[1] = other[1];
-	edgenr = other.edgenr;
-	singedge_left = other.singedge_left;
-	singedge_right = other.singedge_right;
-	seginfo = other.seginfo;
-	si = other.si;
-	domin = other.domin;
-	domout = other.domout;
-	tlosurf = other.tlosurf;
-	geominfo[0] = other.geominfo[0];
-	geominfo[1] = other.geominfo[1];
-	surfnr1 = other.surfnr1;
-	surfnr2 = other.surfnr2;
-	epgeominfo[0] = other.epgeominfo[0];
-	epgeominfo[1] = other.epgeominfo[1];
-	pnums[2] = other.pnums[2];
-	meshdocval = other.meshdocval;
-	hp_elnr = other.hp_elnr;
-        is_curved = other.is_curved;
-      }
-    
-    return *this;
-  }
-  
   void Segment :: DoArchive (Archive & ar)
   {
     string * bcname_dummy = nullptr;

libsrc/meshing/meshtype.hpp

@@ -128,15 +128,7 @@ namespace netgen
     EdgePointGeomInfo ()
       : edgenr(-1), body(0), dist(0.0), u(0.0), v(0.0) { ; }
 
-
+    EdgePointGeomInfo & operator= (const EdgePointGeomInfo & gi2) = default;
-    EdgePointGeomInfo & operator= (const EdgePointGeomInfo & gi2)
-    {
-      edgenr = gi2.edgenr;  
-      body = gi2.body;
-      dist = gi2.dist;
-      u = gi2.u; v = gi2.v;
-      return *this;
-    }
   };
 
   inline ostream & operator<< (ostream & ost, const EdgePointGeomInfo & gi)
@@ -430,8 +422,19 @@ namespace netgen
 
     /// a linked list for all segments in the same face
     SurfaceElementIndex next;
+    ///
+    int hp_elnr;
 
   public:
+    static auto GetDataLayout()
+    {
+      return std::map<string, int>({
+          { "pnum", offsetof(Element2d, pnum)},
+          { "index", offsetof(Element2d, index) },
+          { "np", offsetof(Element2d, np) }
+        });
+    }
+
     ///
     DLL_HEADER Element2d ();
     Element2d (const Element2d &) = default;
@@ -587,6 +590,8 @@ namespace netgen
     void SetOrder (int ox, int oy, int  /* oz */) { orderx = ox; ordery = oy;}
     void SetOrder (int ox, int oy) { orderx = ox; ordery = oy;}
 
+    int GetHpElnr() const { return hp_elnr; }
+    void SetHpElnr(int _hp_elnr) { hp_elnr = _hp_elnr; }
 
     ///
     void GetBox (const T_POINTS & points, Box3d & box) const;
@@ -679,10 +684,6 @@ namespace netgen
     bool operator==(const Element2d & el2) const;
 
     int HasFace(const Element2d& el) const;
-    ///
-    int meshdocval;
-    ///
-    int hp_elnr;
   };
 
   ostream & operator<<(ostream  & s, const Element2d & el);
@@ -719,20 +720,6 @@ namespace netgen
     ELEMENT_TYPE typ;
     /// number of points (4..tet, 5..pyramid, 6..prism, 8..hex, 10..quad tet, 12..quad prism)
     int8_t np;
-    ///
-    class flagstruct { 
-    public:
-      bool marked:1;  // marked for refinement
-      bool badel:1;   // angles worse then limit
-      bool reverse:1; // for refinement a la Bey
-      bool illegal:1; // illegal, will be split or swapped 
-      bool illegal_valid:1; // is illegal-flag valid ?
-      bool badness_valid:1; // is badness valid ?
-      bool refflag:1;     // mark element for refinement
-      bool strongrefflag:1;
-      bool deleted:1;   // element is deleted, will be removed from array
-      bool fixed:1;     // don't change element in optimization
-    };
 
     /// sub-domain index
     int index;
@@ -747,8 +734,32 @@ namespace netgen
     float badness;
     bool is_curved:1;   // element is (high order) curved
 
-  public:
+    class flagstruct {
+    public:
+      bool marked:1;  // marked for refinement
+      bool badel:1;   // angles worse then limit
+      bool reverse:1; // for refinement a la Bey
+      bool illegal:1; // illegal, will be split or swapped
+      bool illegal_valid:1; // is illegal-flag valid ?
+      bool badness_valid:1; // is badness valid ?
+      bool refflag:1;     // mark element for refinement
+      bool strongrefflag:1;
+      bool deleted:1;   // element is deleted, will be removed from array
+      bool fixed:1;     // don't change element in optimization
+    };
+
     flagstruct flags;
+    int hp_elnr;
+  public:
+
+    static auto GetDataLayout()
+    {
+      return std::map<string, int>({
+          { "pnum", offsetof(Element, pnum)},
+          { "index", offsetof(Element, index) },
+          { "np", offsetof(Element, np) }
+        });
+    }
 
     ///
     DLL_HEADER Element () = default;
@@ -764,6 +775,9 @@ namespace netgen
     ///
     // Element & operator= (const Element & el2);
 
+    const flagstruct& Flags() const { return flags; }
+    flagstruct& Flags() { return flags; }
+  
     ///
     DLL_HEADER void SetNP (int anp);
     ///
@@ -909,6 +923,9 @@ namespace netgen
     ///
     DLL_HEADER void Invert ();
 
+    int GetHpElnr() const { return hp_elnr; }
+    void SetHpElnr(int _hp_elnr) { hp_elnr = _hp_elnr; }
+
     /// split into 4 node tets
     void GetTets (NgArray<Element> & locels) const;
     /// split into 4 node tets, local point nrs
@@ -993,7 +1010,6 @@ namespace netgen
     bool IsCurved () const { return is_curved; }
     void SetCurved (bool acurved) { is_curved = acurved; }
 
-    int hp_elnr;
   };
 
   ostream & operator<<(ostream  & s, const Element & el);
@@ -1011,10 +1027,7 @@ namespace netgen
   public:
     ///
     DLL_HEADER Segment();
-    DLL_HEADER Segment (const Segment& other);
+    Segment (const Segment& other) = default;
-
-    ~Segment()
-    { ; }
 
     // friend ostream & operator<<(ostream  & s, const Segment & seg);
 
@@ -1050,10 +1063,8 @@ namespace netgen
     ///
     int meshdocval;
 
-  private:
     bool is_curved;
-
+    int hp_elnr;
-  public:
     /*
       PointIndex operator[] (int i) const
       { return (i == 0) ? p1 : p2; }
@@ -1062,10 +1073,9 @@ namespace netgen
       { return (i == 0) ? p1 : p2; }
     */
 
-    Segment& operator=(const Segment & other);
+    Segment& operator=(const Segment & other) = default;
 
   
-    int hp_elnr;
 
     int GetNP() const
     {
@@ -1172,7 +1182,7 @@ namespace netgen
     void SetDomainIn (int di) { domin = di; }
     void SetDomainOut (int dom) { domout = dom; }
     void SetBCProperty (int bc) { bcprop = bc; }
-    void SetBCName (string * bcn); //  { bcname = bcn; }
+    DLL_HEADER void SetBCName (string * bcn); //  { bcname = bcn; }
     void SetBCName (const string & bcn) { bcname = bcn; }    
     // Philippose - 06/07/2009
     // Set the surface colour

libsrc/meshing/netrule3.cpp

@@ -77,6 +77,7 @@ void vnetrule :: SetFreeZoneTransformation (const Vector & allp, int tolclass)
   fzbox.SetPoint (transfreezone.Elem(1));
   for (i = 2; i <= freezone.Size(); i++)
     fzbox.AddPoint (transfreezone.Elem(i));
+  fzbox.IncreaseRel(1e-8);
   
   
   // MARK(setfz3);

libsrc/meshing/python_mesh.cpp

@@ -434,6 +434,27 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
                                    }))
     ;
 
+  if(ngcore_have_numpy)
+  {
+    auto data_layout = Element::GetDataLayout();
+
+    py::detail::npy_format_descriptor<Element>::register_dtype({
+        py::detail::field_descriptor {
+          "nodes", data_layout["pnum"],
+          ELEMENT_MAXPOINTS * sizeof(PointIndex),
+          py::format_descriptor<int[ELEMENT_MAXPOINTS]>::format(),
+          py::detail::npy_format_descriptor<int[ELEMENT_MAXPOINTS]>::dtype() },
+        py::detail::field_descriptor {
+          "index", data_layout["index"], sizeof(int),
+          py::format_descriptor<int>::format(),
+          py::detail::npy_format_descriptor<int>::dtype() },
+        py::detail::field_descriptor {
+          "np", data_layout["np"], sizeof(int8_t),
+          py::format_descriptor<signed char>::format(),
+          pybind11::dtype("int8") }
+      });
+  }
+
   py::class_<Element2d>(m, "Element2D")
     .def(py::init ([](int index, std::vector<PointIndex> vertices)
                    {
@@ -493,8 +514,29 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
                                    }))
     ;
 
+  if(ngcore_have_numpy)
+  {
+    auto data_layout = Element2d::GetDataLayout();
+    py::detail::npy_format_descriptor<Element2d>::register_dtype({
+        py::detail::field_descriptor {
+          "nodes", data_layout["pnum"],
+          ELEMENT2D_MAXPOINTS * sizeof(PointIndex),
+          py::format_descriptor<int[ELEMENT2D_MAXPOINTS]>::format(),
+          py::detail::npy_format_descriptor<int[ELEMENT2D_MAXPOINTS]>::dtype() },
+        py::detail::field_descriptor {
+          "index", data_layout["index"], sizeof(int),
+          py::format_descriptor<int>::format(),
+          py::detail::npy_format_descriptor<int>::dtype() },
+        py::detail::field_descriptor {
+          "np", data_layout["np"], sizeof(int8_t),
+          py::format_descriptor<signed char>::format(),
+          pybind11::dtype("int8") }
+      });
+  }
+
   py::class_<Segment>(m, "Element1D")
-    .def(py::init([](py::list vertices, py::list surfaces, int index, int edgenr)
+    .def(py::init([](py::list vertices, py::list surfaces, int index, int edgenr,
+                     py::list trignums)
                   {
                     Segment * newel = new Segment();
                     for (int i = 0; i < 2; i++)
@@ -505,6 +547,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
                     newel -> epgeominfo[1].edgenr = edgenr;
                     // needed for codim2 in 3d
                     newel -> edgenr = index;
+                    for(auto i : Range(len(trignums)))
+                      newel->geominfo[i].trignum = py::cast<int>(trignums[i]);
                     if (len(surfaces))
                       {
                         newel->surfnr1 = py::extract<int>(surfaces[0])();
@@ -516,6 +560,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
            py::arg("surfaces")=py::list(),
            py::arg("index")=1,
            py::arg("edgenr")=1,
+           py::arg("trignums")=py::list(), // for stl
          "create segment element"
          )
     .def("__repr__", &ToString<Segment>)
@@ -553,6 +598,20 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
 						     }))
     ;
 
+  if(ngcore_have_numpy)
+  {
+    py::detail::npy_format_descriptor<Segment>::register_dtype({
+        py::detail::field_descriptor {
+          "nodes", offsetof(Segment, pnums),
+          3 * sizeof(PointIndex),
+          py::format_descriptor<int[3]>::format(),
+          py::detail::npy_format_descriptor<int[3]>::dtype() },
+        py::detail::field_descriptor {
+          "index", offsetof(Segment, edgenr), sizeof(int),
+          py::format_descriptor<int>::format(),
+          py::detail::npy_format_descriptor<int>::dtype() },
+      });
+  }
 
   py::class_<Element0d>(m, "Element0D")
     .def(py::init([](PointIndex vertex, int index)
@@ -1191,19 +1250,22 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     
     .def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard<py::gil_scoped_release>())
 
+    .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> boundary,
                               variant<double, py::list> thickness,
                               string material,
                               variant<string, int> domain, bool outside,
                               optional<string> project_boundaries,
-                              bool grow_edges)
+                              bool grow_edges, bool limit_growth_vectors)
            {
              BoundaryLayerParameters blp;
+             BitArray boundaries(self.GetNFD()+1);
+             boundaries.Clear();
              if(int* bc = get_if<int>(&boundary); bc)
                {
                  for (int i = 1; i <= self.GetNFD(); i++)
                    if(self.GetFaceDescriptor(i).BCProperty() == *bc)
-                     blp.surfid.Append (i);
+                     boundaries.SetBit(i);
                }
              else
                {
@@ -1213,19 +1275,29 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
                      auto& fd = self.GetFaceDescriptor(i);
                      if(regex_match(fd.GetBCName(), pattern))
                        {
+                         boundaries.SetBit(i);
                          auto dom_pattern = get_if<string>(&domain);
                          // only add if adjacent to domain
                          if(dom_pattern)
                            {
                              regex pattern(*dom_pattern);
-                             if((fd.DomainIn() > 0 && regex_match(self.GetMaterial(fd.DomainIn()), pattern)) || (fd.DomainOut() > 0 && regex_match(self.GetMaterial(fd.DomainOut()), pattern)))
+                             bool mat1_match = fd.DomainIn() > 0 && regex_match(self.GetMaterial(fd.DomainIn()), pattern);
-                               blp.surfid.Append(i);
+                             bool mat2_match = fd.DomainOut() > 0 && regex_match(self.GetMaterial(fd.DomainOut()), pattern);
+                             // if boundary is inner or outer remove from list
+                             if(mat1_match == mat2_match)
+                               boundaries.Clear(i);
+                             // if((fd.DomainIn() > 0 && regex_match(self.GetMaterial(fd.DomainIn()), pattern)) || (fd.DomainOut() > 0 && regex_match(self.GetMaterial(fd.DomainOut()), pattern)))
+                             // boundaries.Clear(i);
+                             // blp.surfid.Append(i);
                            }
-                         else
+                         // else
-                           blp.surfid.Append(i);
+                         //   blp.surfid.Append(i);
                        }
                      }
                }
+             for(int i = 1; i<=self.GetNFD(); i++)
+               if(boundaries.Test(i))
+                 blp.surfid.Append(i);
              blp.new_mat = material;
 
              if(project_boundaries.has_value())
@@ -1265,12 +1337,13 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
 
              blp.outside = outside;
              blp.grow_edges = grow_edges;
+             blp.limit_growth_vectors = limit_growth_vectors;
 
              GenerateBoundaryLayer (self, blp);
              self.UpdateTopology();
            }, py::arg("boundary"), py::arg("thickness"), py::arg("material"),
           py::arg("domains") = ".*", py::arg("outside") = false,
-          py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true,
+          py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true, py::arg("limit_growth_vectors") = true,
           R"delimiter(
 Add boundary layer to mesh.
 

libsrc/meshing/refine.cpp

@@ -414,7 +414,7 @@ namespace netgen
 	       { 2, 4, 9 },
 	       { 3, 4, 10 } };
 
-	     int elrev = el.flags.reverse;
+	     int elrev = el.Flags().reverse;
 
 	     for (int j = 1; j <= 4; j++)
                pnums.Elem(j) = el.PNum(j);
@@ -480,10 +480,10 @@ namespace netgen
 	      for (int k = 1; k <= 4; k++)
 	        nel.PNum(k) = pnums.Get(reftab[j][k-1]);
 	      nel.SetIndex(ind);
-	      nel.flags.reverse = reverse[j];
+	      nel.Flags().reverse = reverse[j];
 	      if (elrev)
 	      {
-		nel.flags.reverse = !nel.flags.reverse;
+		nel.Flags().reverse = !nel.Flags().reverse;
 		swap (nel.PNum(3), nel.PNum(4));
 	      }
 
@@ -794,10 +794,10 @@ namespace netgen
       if (mesh.VolumeElement(i).Volume(mesh.Points()) < 0)
 	{
 	  wrongels++;
-	  mesh.VolumeElement(i).flags.badel = 1;
+	  mesh.VolumeElement(i).Flags().badel = 1;
 	}
       else
-	mesh.VolumeElement(i).flags.badel = 0;
+	mesh.VolumeElement(i).Flags().badel = 0;
 
     if (wrongels)
       {
@@ -900,14 +900,14 @@ namespace netgen
 		    if (mesh.VolumeElement(i).Volume(mesh.Points()) < 0)
 		      {
 			wrongels++;
-			mesh.VolumeElement(i).flags.badel = 1;
+			mesh.VolumeElement(i).Flags().badel = 1;
 			(*testout) << "wrong el: ";
 			for (int j = 1; j <= 4; j++)
 			  (*testout) << mesh.VolumeElement(i).PNum(j) << " ";
 			(*testout) << endl;
 		      }
 		    else
-		      mesh.VolumeElement(i).flags.badel = 0;
+		      mesh.VolumeElement(i).Flags().badel = 0;
 		  }
 		cout << "wrongels = " << wrongels << endl;
 	      }

libsrc/meshing/secondorder.cpp

@@ -473,10 +473,10 @@ namespace netgen
       if (mesh.VolumeElement(i).CalcJacobianBadness (mesh.Points()) > 1e10)
 	{
 	  wrongels++;
-	  mesh.VolumeElement(i).flags.badel = 1;
+	  mesh.VolumeElement(i).Flags().badel = 1;
 	}
       else
-	mesh.VolumeElement(i).flags.badel = 0;
+	mesh.VolumeElement(i).Flags().badel = 0;
 
     double facok = 0;
     double factry;
@@ -559,7 +559,7 @@ namespace netgen
 		      {
 			wrongels++;
 			Element & el = mesh.VolumeElement(i);
-			el.flags.badel = 1;
+			el.Flags().badel = 1;
 		     
 		      
 			if (lam < 1e-4)
@@ -574,7 +574,7 @@ namespace netgen
 			*/
 		      }
 		    else
-		      mesh.VolumeElement(i).flags.badel = 0;
+		      mesh.VolumeElement(i).Flags().badel = 0;
 		  }
 		cout << "wrongels = " << wrongels << endl;
 	      }
@@ -597,10 +597,10 @@ namespace netgen
 	if (illegalels.Test(i))
 	  {
 	    cout << "illegal element: " << i << endl;
-	    mesh.VolumeElement(i).flags.badel = 1;	
+	    mesh.VolumeElement(i).Flags().badel = 1;
 	  }
 	else
-	  mesh.VolumeElement(i).flags.badel = 0;	
+	  mesh.VolumeElement(i).Flags().badel = 0;
       }
   
     /*

libsrc/meshing/smoothing2.cpp

@@ -756,10 +756,17 @@ namespace netgen
       {
         for (auto & se : mesh.SurfaceElements())
           if (se.GetNP() != 3)
-            {
+          {
-              mixed = true;
+              for(auto pi : se.PNums())
-              break;
+                  if(mesh[pi].Type() == SURFACEPOINT)
-            }
+                  {
+                      mixed = true;
+                      break;
+                  }
+              if(mixed)
+                  break;
+          }
+
         const auto & getDofs = [&] (int i)
           {
             return elementsonpoint[i+PointIndex::BASE];

libsrc/meshing/topology.cpp

@@ -1320,7 +1320,7 @@ namespace netgen
 				  if (mesh->coarsemesh && mesh->hpelements->Size() == mesh->GetNE() )
 				    {
 				      const HPRefElement & hpref_el =
-					(*mesh->hpelements) [ (*mesh)[vertels[k]].hp_elnr];
+					(*mesh->hpelements) [ (*mesh)[vertels[k]].GetHpElnr()];
 				      (*testout) << "coarse eleme = " << hpref_el.coarse_elnr << endl;
 				    }
 

libsrc/occ/occ_edge.cpp

@@ -9,7 +9,6 @@ namespace netgen
 {
     OCCEdge::OCCEdge(TopoDS_Shape edge_, GeometryVertex & start_, GeometryVertex & end_)
         : GeometryEdge(start_, end_),
-          tedge(edge_.TShape()),
           edge(TopoDS::Edge(edge_))
     {
         curve = BRep_Tool::Curve(edge, s0, s1);
@@ -51,7 +50,7 @@ namespace netgen
 
     size_t OCCEdge::GetHash() const
     {
-        return reinterpret_cast<size_t>(tedge.get());
+      return edge.HashCode(std::numeric_limits<Standard_Integer>::max());
     }
 
     void OCCEdge::ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const

libsrc/occ/occ_edge.hpp

@@ -15,7 +15,6 @@ namespace netgen
     class OCCEdge : public GeometryEdge 
     {
         public:
-        T_Shape tedge;
         TopoDS_Edge edge;
         Handle(Geom_Curve) curve;
         double s0, s1;
@@ -25,7 +24,6 @@ namespace netgen
         OCCEdge(TopoDS_Shape edge_, GeometryVertex & start_, GeometryVertex & end_);
 
         auto Shape() const { return edge; }
-        T_Shape TShape() const { return tedge; }
 
         double GetLength() const override;
         Point<3> GetCenter() const override;

libsrc/occ/occ_face.cpp

@@ -9,8 +9,7 @@
 namespace netgen
 {
     OCCFace::OCCFace(TopoDS_Shape dshape)
-        : tface(dshape.TShape()),
+        : face(TopoDS::Face(dshape))
-          face(TopoDS::Face(dshape))
     {
         BRepGProp::SurfaceProperties (dshape, props);
         bbox = ::netgen::GetBoundingBox(face);
@@ -27,7 +26,7 @@ namespace netgen
 
     size_t OCCFace::GetHash() const
     {
-        return reinterpret_cast<size_t>(tface.get());
+      return face.HashCode(std::numeric_limits<Standard_Integer>::max());
     }
 
     Point<3> OCCFace::GetCenter() const
@@ -59,9 +58,9 @@ namespace netgen
         for(auto edge_ : GetEdges(face))
         {
             auto edge = TopoDS::Edge(edge_);
-            if(geom.edge_map.count(edge.TShape())==0)
+            if(geom.edge_map.count(edge)==0)
                 continue;
-            auto edgenr = geom.edge_map[edge.TShape()];
+            auto edgenr = geom.edge_map[edge];
             auto & orientation = edge_orientation[edgenr];
             double s0, s1;
             auto cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1);

libsrc/occ/occ_face.hpp

@@ -13,7 +13,6 @@ namespace netgen
 {
     class OCCFace : public GeometryFace 
     {
-        T_Shape tface;
         TopoDS_Face face;
         GProp_GProps props;
         Box<3> bbox;
@@ -26,7 +25,6 @@ namespace netgen
         OCCFace(TopoDS_Shape dshape);
 
         const TopoDS_Face Shape() const { return face; }
-        T_Shape TShape() { return tface; }
 
         size_t GetHash() const override;
         Point<3> GetCenter() const override;

libsrc/occ/occ_solid.hpp

@@ -10,17 +10,14 @@ namespace netgen
 {
     class OCCSolid : public GeometrySolid 
     {
-        T_Shape tsolid;
         TopoDS_Solid solid;
 
         public:
         OCCSolid(TopoDS_Shape dshape)
-            : tsolid(dshape.TShape()),
+            : solid(TopoDS::Solid(dshape))
-              solid(TopoDS::Solid(dshape))
         { }
 
-        T_Shape TShape() { return tsolid; }
+        size_t GetHash() const override { return solid.HashCode(std::numeric_limits<Standard_Integer>::max()); }
-        size_t GetHash() const override { return reinterpret_cast<size_t>(tsolid.get()); }
     };
 }
 

libsrc/occ/occ_utils.cpp

@@ -6,9 +6,11 @@
 
 namespace netgen
 {
-    Point<3> occ2ng (Handle(TopoDS_TShape) shape)
+    Point<3> occ2ng (const TopoDS_Shape& shape)
     {
-        return occ2ng( Handle(BRep_TVertex)::DownCast(shape)->Pnt() );
+      if(shape.ShapeType() != TopAbs_VERTEX)
+        throw Exception("Try to convert non vertex to point!");
+      return occ2ng( BRep_Tool::Pnt(TopoDS::Vertex(shape)) );
     }
 
     Transformation<3> occ2ng (const gp_Trsf & occ_trafo)

libsrc/occ/occ_utils.hpp

@@ -22,8 +22,6 @@
 
 namespace netgen
 {
-    typedef Handle(TopoDS_TShape) T_Shape;
-
     inline Point<3> occ2ng (const gp_Pnt & p)
     {
         return Point<3> (p.X(), p.Y(), p.Z());
@@ -39,12 +37,7 @@ namespace netgen
         return Vec<3> (v.X(), v.Y(), v.Z());
     }
 
-    DLL_HEADER Point<3> occ2ng (T_Shape shape);
+    DLL_HEADER Point<3> occ2ng (const TopoDS_Shape & s);
-
-    inline Point<3> occ2ng (const TopoDS_Shape & s)
-    {
-        return occ2ng(s.TShape());
-    }
 
     inline Point<3> occ2ng (const TopoDS_Vertex & v)
     {
@@ -70,8 +63,8 @@ namespace netgen
     class OCCIdentification
     {
     public:
-      T_Shape from;
+      TopoDS_Shape from;
-      T_Shape to;
+      TopoDS_Shape to;
       Transformation<3> trafo;
       string name;
       Identifications::ID_TYPE type;
@@ -134,14 +127,6 @@ namespace netgen
       return IndexMapIterator(indmap);
     }
 
-    struct ShapeLess
-    {
-      bool operator() (const TopoDS_Shape& s1, const TopoDS_Shape& s2) const
-      {
-        return s1.TShape() < s2.TShape();
-      }
-    };
-
     class ListOfShapes : public std::vector<TopoDS_Shape>
     {
     public:
@@ -297,7 +282,12 @@ namespace netgen
     GProp_GProps props;
     switch (shape.ShapeType())
       {
+      case TopAbs_SOLID:
+      case TopAbs_COMPOUND:
+      case TopAbs_COMPSOLID:
+        BRepGProp::VolumeProperties (shape, props); break;
       case TopAbs_FACE:
+      case TopAbs_SHELL:
         BRepGProp::SurfaceProperties (shape, props); break;
       default:
         BRepGProp::LinearProperties(shape, props);

libsrc/occ/occ_vertex.cpp

@@ -7,8 +7,7 @@ namespace netgen
 {
 
     OCCVertex::OCCVertex( TopoDS_Shape s )
-        : vertex(TopoDS::Vertex(s)),
+        : vertex(TopoDS::Vertex(s))
-          tvertex(s.TShape())
     {
         p = occ2ng(vertex);
     }
@@ -20,6 +19,6 @@ namespace netgen
 
     size_t OCCVertex::GetHash() const
     {
-        return reinterpret_cast<size_t>(tvertex.get());
+      return vertex.HashCode(std::numeric_limits<Standard_Integer>::max());
     }
 }

libsrc/occ/occ_vertex.hpp

@@ -12,7 +12,6 @@ namespace netgen
     class OCCVertex : public GeometryVertex
     {
         TopoDS_Vertex vertex;
-        T_Shape tvertex;
         Point<3> p;
 
         public:
@@ -21,7 +20,6 @@ namespace netgen
         ~OCCVertex() {}
         Point<3> GetPoint() const override;
         size_t GetHash() const override;
-        T_Shape TShape() { return tvertex; }
     };
 }
 

libsrc/occ/occgenmesh.cpp

@@ -252,7 +252,6 @@ namespace netgen
     FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
     auto face = TopoDS::Face(geom.fmap(k));
     const auto& occface = dynamic_cast<const OCCFace&>(geom.GetFace(k-1));
-    auto fshape = face.TShape();
 
     int oldnf = mesh.GetNSE();
 
@@ -403,11 +402,11 @@ namespace netgen
 
 
     // Philippose - 15/01/2009
-    double maxh = min2(geom.face_maxh[k-1], OCCGeometry::global_shape_properties[TopoDS::Face(geom.fmap(k)).TShape()].maxh);
+    double maxh = min2(geom.face_maxh[k-1], OCCGeometry::global_shape_properties[geom.fmap(k)].maxh);
     //double maxh = mparam.maxh;
     //      int noldpoints = mesh->GetNP();
     int noldsurfel = mesh.GetNSE();
-    int layer = OCCGeometry::global_shape_properties[TopoDS::Face(geom.fmap(k)).TShape()].layer;
+    int layer = OCCGeometry::global_shape_properties[geom.fmap(k)].layer;
 
     static Timer tsurfprop("surfprop");
     tsurfprop.Start();
@@ -475,8 +474,8 @@ namespace netgen
     int dom = 0;
     for (TopExp_Explorer e(geom.GetShape(), TopAbs_SOLID); e.More(); e.Next(), dom++)
     {
-      maxhdom[dom] = min2(maxhdom[dom], OCCGeometry::global_shape_properties[e.Current().TShape()].maxh);
+      maxhdom[dom] = min2(maxhdom[dom], OCCGeometry::global_shape_properties[e.Current()].maxh);
-      maxlayer = max2(maxlayer, OCCGeometry::global_shape_properties[e.Current().TShape()].layer);
+      maxlayer = max2(maxlayer, OCCGeometry::global_shape_properties[e.Current()].layer);
     }
 
 
@@ -519,7 +518,7 @@ namespace netgen
         for (int i = 1; i <= nedges && !multithread.terminate; i++)
           {
             TopoDS_Edge e = TopoDS::Edge (geom.emap(i));
-            int layer = OCCGeometry::global_shape_properties[e.TShape()].layer;
+            int layer = OCCGeometry::global_shape_properties[e].layer;
             multithread.percent = 100 * (i-1)/double(nedges);
             if (BRep_Tool::Degenerated(e)) continue;
 
@@ -535,7 +534,7 @@ namespace netgen
 
             bool is_identified_edge = false;
             // TODO: change to use hash value
-            const auto& gedge = geom.GetEdge(geom.edge_map.at(e.TShape()));
+            const auto& gedge = geom.GetEdge(geom.edge_map.at(e));
             auto& v0 = gedge.GetStartVertex();
             auto& v1 = gedge.GetEndVertex();
             for(auto & ident : v0.identifications)
@@ -565,12 +564,12 @@ namespace netgen
                 int face_index = geom.fmap.FindIndex(parent_face);
 
                 if(face_index >= 1) localh = min(localh,geom.face_maxh[face_index - 1]);
-                localh = min2(localh, OCCGeometry::global_shape_properties[parent_face.TShape()].maxh);
+                localh = min2(localh, OCCGeometry::global_shape_properties[parent_face].maxh);
               }
 
             Handle(Geom_Curve) c = BRep_Tool::Curve(e, s0, s1);
 
-            localh = min2(localh, OCCGeometry::global_shape_properties[e.TShape()].maxh);
+            localh = min2(localh, OCCGeometry::global_shape_properties[e].maxh);
             maxedgelen = max (maxedgelen, len);
             minedgelen = min (minedgelen, len);
             int maxj = max((int) ceil(len/localh), 2);
@@ -593,7 +592,7 @@ namespace netgen
             double maxcur = 0;
             multithread.percent = 100 * (i-1)/double(nedges);
             TopoDS_Edge edge = TopoDS::Edge (geom.emap(i));
-            int layer = OCCGeometry::global_shape_properties[edge.TShape()].layer;
+            int layer = OCCGeometry::global_shape_properties[edge].layer;
             if (BRep_Tool::Degenerated(edge)) continue;
             double s0, s1;
             Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1);
@@ -628,7 +627,7 @@ namespace netgen
           {
             multithread.percent = 100 * (i-1)/double(nfaces);
             TopoDS_Face face = TopoDS::Face(geom.fmap(i));
-            int layer = OCCGeometry::global_shape_properties[face.TShape()].layer;
+            int layer = OCCGeometry::global_shape_properties[face].layer;
             TopLoc_Location loc;
             Handle(Geom_Surface) surf = BRep_Tool::Surface (face);
             Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc);
@@ -694,7 +693,7 @@ namespace netgen
             for (int i = 1; i <= nedges && !multithread.terminate; i++)
               {
                 TopoDS_Edge edge = TopoDS::Edge (geom.emap(i));
-                int layer = OCCGeometry::global_shape_properties[edge.TShape()].layer;
+                int layer = OCCGeometry::global_shape_properties[edge].layer;
                 if (BRep_Tool::Degenerated(edge)) continue;
 
                 double s0, s1;

libsrc/occ/occgeom.cpp

@@ -71,8 +71,8 @@ namespace netgen
   void LoadOCCInto(OCCGeometry* occgeo, const filesystem::path & filename);
   void PrintContents (OCCGeometry * geom);
 
-  std::map<Handle(TopoDS_TShape), ShapeProperties> OCCGeometry::global_shape_properties;
+  std::map<TopoDS_Shape, ShapeProperties> OCCGeometry::global_shape_properties;
-  std::map<Handle(TopoDS_TShape), std::vector<OCCIdentification>> OCCGeometry::identifications;
+  std::map<TopoDS_Shape, std::vector<OCCIdentification>> OCCGeometry::identifications;
   
   TopoDS_Shape ListOfShapes::Max(gp_Vec dir)
   {
@@ -125,7 +125,7 @@ namespace netgen
 
   ListOfShapes ListOfShapes::SubShapes(TopAbs_ShapeEnum type) const
   {
-    std::set<TopoDS_Shape, ShapeLess> unique_shapes;
+    std::set<TopoDS_Shape> unique_shapes;
     for(const auto& shape : *this)
       for(TopExp_Explorer e(shape, type); e.More(); e.Next())
         unique_shapes.insert(e.Current());
@@ -200,7 +200,7 @@ namespace netgen
   {
     MeshingParameters local_mp = mparam;
     auto face = TopoDS::Face(fmap(nr+1));
-    if(auto quad_dominated = OCCGeometry::global_shape_properties[face.TShape()].quad_dominated; quad_dominated.has_value())
+    if(auto quad_dominated = OCCGeometry::global_shape_properties[face].quad_dominated; quad_dominated.has_value())
       local_mp.quad = *quad_dominated;
 
     bool failed = OCCMeshFace(*this, mesh, glob2loc, local_mp, nr, PARAMETERSPACE, true);
@@ -376,9 +376,9 @@ namespace netgen
     
     for (TopExp_Explorer e(shape, TopAbs_SOLID); e.More(); e.Next())
       {
-        if (auto name = OCCGeometry::global_shape_properties[e.Current().TShape()].name)
+        if (auto name = OCCGeometry::global_shape_properties[e.Current()].name)
           for (auto mods : history->Modified(e.Current()))
-            OCCGeometry::global_shape_properties[mods.TShape()].name = *name;
+            OCCGeometry::global_shape_properties[mods].name = *name;
       }
 #endif // OCC_HAVE_HISTORY    
     
@@ -444,7 +444,7 @@ namespace netgen
          for (exp0.Init (shape, TopAbs_FACE); exp0.More(); exp0.Next())
          {
             TopoDS_Face face = TopoDS::Face (exp0.Current());
-            auto props = global_shape_properties[face.TShape()];
+            auto props = global_shape_properties[face];
 
             sff = new ShapeFix_Face (face);
             sff->FixAddNaturalBoundMode() = Standard_True;
@@ -475,7 +475,7 @@ namespace netgen
 
             // Set the original properties of the face to the newly created 
             // face (after the healing process)
-            global_shape_properties[face.TShape()];
+            global_shape_properties[face];
          }
          shape = rebuild->Apply(shape);
       }
@@ -1122,54 +1122,51 @@ namespace netgen
       for(auto i1 : Range(1, vmap.Extent()+1))
       {
           auto v = vmap(i1);
-          auto tshape = v.TShape();
+          if(vertex_map.count(v)!=0)
-          if(vertex_map.count(tshape)!=0)
               continue;
           auto occ_vertex = make_unique<OCCVertex>(TopoDS::Vertex(v));
           occ_vertex->nr = vertices.Size();
-          vertex_map[tshape] = occ_vertex->nr;
+          vertex_map[v] = occ_vertex->nr;
 
-          if(global_shape_properties.count(tshape)>0)
+          if(global_shape_properties.count(v)>0)
-              occ_vertex->properties = global_shape_properties[tshape];
+              occ_vertex->properties = global_shape_properties[v];
           vertices.Append(std::move(occ_vertex));
       }
 
       for(auto i1 : Range(1, emap.Extent()+1))
       {
           auto e = emap(i1);
-          auto tshape = e.TShape();
           auto edge = TopoDS::Edge(e);
-          if(edge_map.count(tshape)!=0)
+          if(edge_map.count(e)!=0)
               continue;
-          edge_map[tshape] = edges.Size();
+          edge_map[e] = edges.Size();
           auto verts = GetVertices(e);
-          auto occ_edge = make_unique<OCCEdge>(edge, *vertices[vertex_map[verts[0].TShape()]], *vertices[vertex_map[verts[1].TShape()]] );
+          auto occ_edge = make_unique<OCCEdge>(edge, *vertices[vertex_map[verts[0]]], *vertices[vertex_map[verts[1]]] );
-          occ_edge->properties = global_shape_properties[tshape];
+          occ_edge->properties = global_shape_properties[e];
           edges.Append(std::move(occ_edge));
       }
 
       for(auto i1 : Range(1, fmap.Extent()+1))
       {
           auto f = fmap(i1);
-          auto tshape = f.TShape();
+          if(face_map.count(f)==0)
-          if(face_map.count(tshape)==0)
           {
 
               auto k = faces.Size();
-              face_map[tshape] = k;
+              face_map[f] = k;
               auto occ_face = make_unique<OCCFace>(f);
 
               for(auto e : GetEdges(f))
-                  occ_face->edges.Append( edges[edge_map[e.TShape()]].get() );
+                  occ_face->edges.Append( edges[edge_map[e]].get() );
 
-              if(global_shape_properties.count(tshape)>0)
+              if(global_shape_properties.count(f)>0)
-                  occ_face->properties = global_shape_properties[tshape];
+                  occ_face->properties = global_shape_properties[f];
               faces.Append(std::move(occ_face));
 
               if(dimension==2)
                   for(auto e : GetEdges(f))
                   {
-                      auto & edge = *edges[edge_map[e.TShape()]];
+                      auto & edge = *edges[edge_map[e]];
                       if(e.Orientation() == TopAbs_REVERSED)
                           edge.domout = k;
                       else
@@ -1182,21 +1179,20 @@ namespace netgen
       for(auto i1 : Range(1, somap.Extent()+1))
       {
           auto s = somap(i1);
-          auto tshape = s.TShape();
           int k;
-          if(solid_map.count(tshape)==0)
+          if(solid_map.count(s)==0)
           {
               k = solids.Size();
-              solid_map[tshape] = k;
+              solid_map[s] = k;
               auto occ_solid = make_unique<OCCSolid>(s);
-              if(global_shape_properties.count(tshape)>0)
+              if(global_shape_properties.count(s)>0)
-                  occ_solid->properties = global_shape_properties[tshape];
+                  occ_solid->properties = global_shape_properties[s];
               solids.Append(std::move(occ_solid));
           }
 
           for(auto f : GetFaces(s))
           {
-              auto face_nr = face_map[f.TShape()];
+              auto face_nr = face_map[f];
               auto & face = faces[face_nr];
               if(face->domin==-1)
                   face->domin = k;
@@ -1208,9 +1204,9 @@ namespace netgen
       // Add identifications
       auto add_identifications = [&](auto & shapes, auto & shape_map)
       {
-          for(auto &[tshape, nr] : shape_map)
+          for(auto &[shape, nr] : shape_map)
-            if(identifications.count(tshape))
+            if(identifications.count(shape))
-                for(auto & ident : identifications[tshape])
+                for(auto & ident : identifications[shape])
                 {
                     if(shape_map.count(ident.from)==0 || shape_map.count(ident.to)==0)
                         continue;
@@ -1321,7 +1317,7 @@ namespace netgen
     Array<GeometryVertex*> verts;
     const auto& occface = dynamic_cast<const OCCFace&>(face);
     for(auto& vert : GetVertices(occface.Shape()))
-      verts.Append(vertices[vertex_map.at(vert.TShape())].get());
+      verts.Append(vertices[vertex_map.at(vert)].get());
     return move(verts);
   }
 
@@ -1589,7 +1585,11 @@ namespace netgen
     if(ar.Output())
       {
         std::stringstream ss;
+#if OCC_VERSION_HEX < 0x070600
         BRepTools::Write(shape, ss);
+#else
+        BRepTools::Write(shape, ss, false, false, TopTools_FormatVersion_VERSION_1);
+#endif
         ar << ss.str();
       }
     else
@@ -1611,34 +1611,33 @@ namespace netgen
         auto occ_hash = key.HashCode(1<<31UL);
         return std::hash<decltype(occ_hash)>()(occ_hash);
     };
-    std::map<Handle(TopoDS_TShape), int> tshape_map;
+    std::map<TopoDS_Shape, int> shape_map;
-    Array<Handle(TopoDS_TShape)> tshape_list;
+    Array<TopoDS_Shape> shape_list;
 
     ar & dimension;
     for (auto typ : { TopAbs_SOLID, TopAbs_FACE,  TopAbs_EDGE })
       for (TopExp_Explorer e(shape, typ); e.More(); e.Next())
         {
           auto ds = e.Current(); 
-          auto ts = ds.TShape();
+          if(shape_map.count(ds)==0)
-          if(tshape_map.count(ts)==0)
             {
-              tshape_map[ts] = tshape_list.Size();
+              shape_map[ds] = shape_list.Size();
-              tshape_list.Append(ts);
+              shape_list.Append(ds);
             }
         }
 
-    for (auto ts : tshape_list)
+    for (auto s : shape_list)
       {
-        bool has_properties = global_shape_properties.count(ts);
+        bool has_properties = global_shape_properties.count(s);
         ar & has_properties;
         if(has_properties)
-            ar & global_shape_properties[ts];
+            ar & global_shape_properties[s];
 
-        bool has_identifications = identifications.count(ts);
+        bool has_identifications = identifications.count(s);
         ar & has_identifications;
         if(has_identifications)
           {
-            auto & idents = identifications[ts];
+            auto & idents = identifications[s];
             auto n_idents = idents.size();
             ar & n_idents;
             idents.resize(n_idents);
@@ -1648,14 +1647,14 @@ namespace netgen
                 int id_from, id_to;
                 if(ar.Output())
                 {
-                    id_from = tshape_map[id.from];
+                    id_from = shape_map[id.from];
-                    id_to = tshape_map[id.to];
+                    id_to = shape_map[id.to];
                 }
                 ar & id_from & id_to & id.trafo & id.name;
                 if(ar.Input())
                 {
-                    id.from = tshape_list[id_from];
+                    id.from = shape_list[id_from];
-                    id.to = tshape_list[id_to];
+                    id.to = shape_list[id_to];
                 }
               }
           }
@@ -1977,7 +1976,7 @@ namespace netgen
       if(tree.GetTolerance() < Dist(trafo(c_me), c_you))
           return false;
 
-      std::map<T_Shape, T_Shape> vmap;
+      std::map<TopoDS_Shape, optional<TopoDS_Shape>> vmap;
 
       auto verts_me = GetVertices(me);
       auto verts_you = GetVertices(you);
@@ -1987,21 +1986,21 @@ namespace netgen
 
       for (auto i : Range(verts_me.size()))
       {
-          auto s = verts_me[i].TShape();
+          auto s = verts_me[i];
           if(vmap.count(s)>0)
               continue;
           auto p = trafo(occ2ng(s));
           tree.Insert( p, i );
-          vmap[s] = nullptr;
+          vmap[s] = nullopt;
       }
           
       for (auto vert : verts_you)
       {
-          auto s = vert.TShape();
+          auto s = vert;
           auto p = occ2ng(s);
           bool vert_mapped = false;
           tree.GetFirstIntersecting( p, p, [&](size_t i ) {
-                  vmap[verts_me[i].TShape()] = s;
+                  vmap[verts_me[i]] = s;
                   vert_mapped = true;
                   return true;
           });
@@ -2057,8 +2056,8 @@ namespace netgen
             if(!IsMappedShape(trafo, shape_me, shape_you))
                 continue;
 
-            OCCGeometry::identifications[shape_me.TShape()].push_back
+            OCCGeometry::identifications[shape_me].push_back
-                (OCCIdentification { shape_me.TShape(), shape_you.TShape(), trafo, name, type });
+                (OCCIdentification { shape_me, shape_you, trafo, name, type });
         }
   }
 
@@ -2120,7 +2119,7 @@ namespace netgen
             XCAFPrs_Style aStyle;
             set.FindFromKey(e.Current(), aStyle);
 
-            auto & prop = OCCGeometry::global_shape_properties[e.Current().TShape()];
+            auto & prop = OCCGeometry::global_shape_properties[e.Current()];
             if(aStyle.IsSetColorSurf())
                 prop.col = step_utils::ReadColor(aStyle.GetColorSurfRGBA());
           }
@@ -2140,7 +2139,7 @@ namespace netgen
             if (!transProc->IsBound(item))
               continue;
 
-            OCCGeometry::global_shape_properties[shape.TShape()].name = name;
+            OCCGeometry::global_shape_properties[shape].name = name;
           }
 
 
@@ -2164,7 +2163,7 @@ namespace netgen
             if(name != "netgen_geometry_properties")
                 continue;
 
-            auto & prop = OCCGeometry::global_shape_properties[shape.TShape()];
+            auto & prop = OCCGeometry::global_shape_properties[shape];
 
             auto nprops = item->NbItemElement();
 
@@ -2190,7 +2189,7 @@ namespace netgen
           Handle(StepRepr_RepresentationItem) item = STEPConstruct::FindEntity(finder, shape);
           if(!item)
               return;
-          auto prop = OCCGeometry::global_shape_properties[shape.TShape()];
+          auto prop = OCCGeometry::global_shape_properties[shape];
 
           if(auto n = prop.name)
               item->SetName(MakeName(*n));
@@ -2219,7 +2218,7 @@ namespace netgen
       void WriteIdentifications(const Handle(Interface_InterfaceModel) model, const TopoDS_Shape & shape, const Handle(Transfer_FinderProcess) finder)
       {
           Handle(StepRepr_RepresentationItem) item = STEPConstruct::FindEntity(finder, shape);
-          auto & identifications = OCCGeometry::identifications[shape.TShape()];
+          auto & identifications = OCCGeometry::identifications[shape];
           if(identifications.size()==0)
               return;
           auto n = identifications.size();
@@ -2270,7 +2269,7 @@ namespace netgen
 
               result.push_back(ident);
           }
-          OCCGeometry::identifications[shape_origin.TShape()] = result;
+          OCCGeometry::identifications[shape_origin] = result;
       }
 
       void WriteSTEP(const TopoDS_Shape & shape, const filesystem::path & filename)
@@ -2294,7 +2293,7 @@ namespace netgen
           for (auto typ : { TopAbs_SOLID, TopAbs_FACE,  TopAbs_EDGE })
             for (TopExp_Explorer e(shape, typ); e.More(); e.Next())
               {
-                auto prop = OCCGeometry::global_shape_properties[e.Current().TShape()];
+                auto prop = OCCGeometry::global_shape_properties[e.Current()];
                 if(auto col = prop.col)
                     colortool->SetColor(e.Current(), step_utils::MakeColor(*col), XCAFDoc_ColorGen);
               }

libsrc/occ/occgeom.hpp

@@ -31,6 +31,20 @@
 #define OCC_HAVE_HISTORY
 #endif
 
+namespace std
+{
+  template<>
+  struct less<TopoDS_Shape>
+  {
+    bool operator() (const TopoDS_Shape& s1, const TopoDS_Shape& s2) const
+    {
+      return s1.HashCode(std::numeric_limits<Standard_Integer>::max()) <
+        s2.HashCode(std::numeric_limits<Standard_Integer>::max());
+    }
+  };
+}
+
+
 namespace netgen
 {
 
@@ -135,15 +149,15 @@ namespace netgen
     Point<3> center;
     OCCParameters occparam;
   public:
-    static std::map<T_Shape, ShapeProperties> global_shape_properties;
+    static std::map<TopoDS_Shape, ShapeProperties> global_shape_properties;
-    static std::map<T_Shape, std::vector<OCCIdentification>> identifications;
+    static std::map<TopoDS_Shape, std::vector<OCCIdentification>> identifications;
 
     TopoDS_Shape shape;
     TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap; // legacy maps
     NgArray<bool> fsingular, esingular, vsingular;
     Box<3> boundingbox;
 
-    std::map<T_Shape, int> edge_map, vertex_map, face_map, solid_map;
+    std::map<TopoDS_Shape, int> solid_map, face_map, edge_map, vertex_map;
 
     mutable int changed;
     mutable NgArray<int> facemeshstatus;
@@ -188,6 +202,12 @@ namespace netgen
     Mesh::GEOM_TYPE GetGeomType() const override
     { return Mesh::GEOM_OCC; }
 
+    void SetDimension(int dim)
+    {
+        dimension = dim;
+        BuildFMap();
+    }
+
     void SetOCCParameters(const OCCParameters& par)
     { occparam = par; }
 
@@ -376,8 +396,8 @@ namespace netgen
   template <class TBuilder>
   void PropagateIdentifications (TBuilder & builder, TopoDS_Shape shape, std::optional<Transformation<3>> trafo = nullopt)
   {
-    std::map<T_Shape, std::set<T_Shape>> mod_map;
+    std::map<TopoDS_Shape, std::set<TopoDS_Shape>> mod_map;
-    std::map<T_Shape, bool> tshape_handled;
+    std::map<TopoDS_Shape, bool> shape_handled;
     Transformation<3> trafo_inv;
     if(trafo)
         trafo_inv = trafo->CalcInverse();
@@ -385,35 +405,34 @@ namespace netgen
     for (auto typ : { TopAbs_SOLID, TopAbs_FACE,  TopAbs_EDGE, TopAbs_VERTEX })
       for (TopExp_Explorer e(shape, typ); e.More(); e.Next())
       {
-          auto tshape = e.Current().TShape();
+          auto s = e.Current();
-          mod_map[tshape].insert(tshape);
+          mod_map[s].insert(s);
-          tshape_handled[tshape] = false;
+          shape_handled[s] = false;
       }
   
     for (auto typ : { TopAbs_SOLID, TopAbs_FACE,  TopAbs_EDGE, TopAbs_VERTEX })
       for (TopExp_Explorer e(shape, typ); e.More(); e.Next())
         {
-          auto tshape = e.Current().TShape();
+          auto s = e.Current();
-  
+            for (auto mods : builder.Modified(s))
-          for (auto mods : builder.Modified(e.Current()))
+              mod_map[s].insert(mods);
-              mod_map[tshape].insert(mods.TShape());
         }
   
     for (auto typ : { TopAbs_SOLID, TopAbs_FACE,  TopAbs_EDGE, TopAbs_VERTEX })
       for (TopExp_Explorer e(shape, typ); e.More(); e.Next())
       {
-          auto tshape = e.Current().TShape();
+          auto s = e.Current();
   
-          if(tshape_handled[tshape])
+          if(shape_handled[s])
               continue;
-          tshape_handled[tshape] = true;
+          shape_handled[s] = true;
   
-          if(OCCGeometry::identifications.count(tshape)==0)
+          if(OCCGeometry::identifications.count(s)==0)
               continue;
   
-          auto tshape_mapped = mod_map[tshape];
+          auto shape_mapped = mod_map[s];
   
-          for(auto ident : OCCGeometry::identifications[tshape])
+          for(auto ident : OCCGeometry::identifications[s])
           {
               // nothing happened
               if(mod_map[ident.to].size()==1 && mod_map[ident.from].size() ==1)
@@ -425,11 +444,8 @@ namespace netgen
               for(auto from_mapped : mod_map[from])
                   for(auto to_mapped : mod_map[to])
                   {
-                      if(from==from_mapped && to==to_mapped)
+                      if(from.IsSame(from_mapped) && to.IsSame(to_mapped))
-                          continue;
+                        continue;
-  
-                      TopoDS_Shape s_from; s_from.TShape(from_mapped);
-                      TopoDS_Shape s_to; s_to.TShape(to_mapped);
   
                       Transformation<3> trafo_mapped = ident.trafo;
                       if(trafo)
@@ -439,14 +455,14 @@ namespace netgen
                           trafo_mapped.Combine(*trafo, trafo_temp);
                       }
   
-                      if(!IsMappedShape(trafo_mapped, s_from, s_to))
+                      if(!IsMappedShape(trafo_mapped, from_mapped, to_mapped))
                           continue;
   
                       OCCIdentification id_new = ident;
                       id_new.to = to_mapped;
                       id_new.from = from_mapped;
                       id_new.trafo = trafo_mapped;
-                      auto id_owner = from == tshape ? from_mapped : to_mapped;
+                      auto id_owner = from.IsSame(s) ? from_mapped : to_mapped;
                       OCCGeometry::identifications[id_owner].push_back(id_new);
                   }
           }
@@ -461,11 +477,11 @@ namespace netgen
     for (auto typ : { TopAbs_SOLID, TopAbs_FACE,  TopAbs_EDGE })
       for (TopExp_Explorer e(shape, typ); e.More(); e.Next())
         {
-          auto tshape = e.Current().TShape();
+          auto s = e.Current();
-          auto & prop = OCCGeometry::global_shape_properties[tshape];
+          auto & prop = OCCGeometry::global_shape_properties[s];
-          for (auto mods : builder.Modified(e.Current()))
+          for (auto mods : builder.Modified(s))
-            OCCGeometry::global_shape_properties[mods.TShape()].Merge(prop);
+            OCCGeometry::global_shape_properties[mods].Merge(prop);
-          have_identifications |= OCCGeometry::identifications.count(tshape) > 0;
+          have_identifications |= OCCGeometry::identifications.count(s) > 0;
         }
     if(have_identifications)
         PropagateIdentifications(builder, shape, trafo);

libsrc/occ/python_occ.cpp

@@ -26,6 +26,7 @@ using namespace netgen;
 namespace netgen
 {
   extern std::shared_ptr<NetgenGeometry> ng_geometry;
+  extern std::shared_ptr<Mesh> mesh;
 }
 
 static string occparameter_description = R"delimiter(
@@ -117,11 +118,11 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                     
                     for (auto & s : shapes)
                       for (TopExp_Explorer e(s, TopAbs_SOLID); e.More(); e.Next())
-                        if (auto name = OCCGeometry::global_shape_properties[e.Current().TShape()].name)
+                        if (auto name = OCCGeometry::global_shape_properties[e.Current()].name)
                           {
                             TopTools_ListOfShape modlist = history->Modified(e.Current());
                             for (auto mods : modlist)
-                              OCCGeometry::global_shape_properties[mods.TShape()].name = *name;
+                              OCCGeometry::global_shape_properties[mods].name = *name;
                           }
 #endif // OCC_HAVE_HISTORY
 
@@ -133,7 +134,7 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                   }), py::arg("shape"),
          "Create Netgen OCCGeometry from existing TopoDS_Shape")
     
-    .def(py::init([] (const string& filename)
+    .def(py::init([] (const string& filename, int dim)
                   {
                     shared_ptr<OCCGeometry> geo;
                     if(EndsWith(filename, ".step") || EndsWith(filename, ".stp"))
@@ -144,9 +145,11 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                       geo.reset(LoadOCC_IGES(filename));
                     else
                       throw Exception("Cannot load file " + filename + "\nValid formats are: step, stp, brep, iges");
+                    if(dim<3)
+                      geo->SetDimension(dim);
                     ng_geometry = geo;
                     return geo;
-                  }), py::arg("filename"),
+                  }), py::arg("filename"), py::arg("dim")=3,
         "Load OCC geometry from step, brep or iges file")
     .def(NGSPickle<OCCGeometry>())
     .def("Glue", &OCCGeometry::GlueGeometry)
@@ -247,7 +250,8 @@ DLL_HEADER void ExportNgOCC(py::module &m)
             return res;
          }, py::call_guard<py::gil_scoped_release>())
     .def("GenerateMesh", [](shared_ptr<OCCGeometry> geo,
-                            MeshingParameters* pars, NgMPI_Comm comm, py::kwargs kwargs)
+                            MeshingParameters* pars, NgMPI_Comm comm,
+                            shared_ptr<Mesh> mesh, py::kwargs kwargs)
                          {
                            MeshingParameters mp;
                            OCCParameters occparam;
@@ -263,7 +267,8 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                              CreateMPfromKwargs(mp, kwargs);
                            }
                            geo->SetOCCParameters(occparam);
-                           auto mesh = make_shared<Mesh>();
+                           if(!mesh)
+                             mesh = make_shared<Mesh>();
                            mesh->SetCommunicator(comm);
                            mesh->SetGeometry(geo);
 
@@ -272,7 +277,10 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                                SetGlobalMesh(mesh);
                                auto result = geo->GenerateMesh(mesh, mp);
                                if(result != 0)
-                                 throw Exception("Meshing failed!");
+                                 {
+                                   netgen::mesh = mesh;   // keep mesh for debugging
+                                   throw Exception("Meshing failed!");
+                                 }
                                ng_geometry = geo;
                                if (comm.Size() > 1)
                                  mesh->Distribute();
@@ -283,7 +291,7 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                              }
                            return mesh;
                          }, py::arg("mp") = nullptr, py::arg("comm")=NgMPI_Comm{},
-      py::call_guard<py::gil_scoped_release>(),
+         py::arg("mesh")=nullptr, py::call_guard<py::gil_scoped_release>(),
          (meshingparameter_description + occparameter_description).c_str())
     .def_property_readonly("shape", [](const OCCGeometry & self) { return self.GetShape(); })
     ;
@@ -319,8 +327,6 @@ DLL_HEADER void ExportNgOCC(py::module &m)
       Handle(XCAFDoc_MaterialTool) material_tool = XCAFDoc_DocumentTool::MaterialTool(doc->Main());
       // Handle(XCAFDoc_VisMaterialTool) vismaterial_tool = XCAFDoc_DocumentTool::VisMaterialTool(doc->Main());
 
-      cout << "handle(shape) = " << *(void**)(void*)(&(shape.TShape())) << endl;
-      
       // TDF_LabelSequence doc_shapes;
       // shape_tool->GetShapes(doc_shapes);
       // cout << "shape tool nbentities: " << doc_shapes.Size() << endl;

libsrc/occ/python_occ_basic.cpp

@@ -76,6 +76,8 @@ DLL_HEADER void ExportNgOCCBasic(py::module &m)
     .def_property("x", [](gp_Vec&p) { return p.X(); }, [](gp_Vec&p,double x) { p.SetX(x); })
     .def_property("y", [](gp_Vec&p) { return p.Y(); }, [](gp_Vec&p,double y) { p.SetY(y); })
     .def_property("z", [](gp_Vec&p) { return p.Z(); }, [](gp_Vec&p,double z) { p.SetZ(z); })
+    .def("Norm", [](const gp_Vec& v)
+    { return v.Magnitude(); })
     .def("__str__", [] (const gp_Vec & p) {
         stringstream str;
         str << "(" << p.X() << ", " << p.Y() << ", " << p.Z() << ")";

libsrc/occ/python_occ_shapes.cpp

@@ -304,7 +304,7 @@ public:
 
     // auto edge = BRepBuilderAPI_MakeEdge(curve).Edge();
     if (name)
-      OCCGeometry::global_shape_properties[edge.TShape()].name = name;
+      OCCGeometry::global_shape_properties[edge].name = name;
     wire_builder.Add(edge);
 
     if (closing) Finish();
@@ -591,7 +591,7 @@ public:
   auto NameVertex (string name)
   {
     if (!lastvertex.IsNull())
-      OCCGeometry::global_shape_properties[lastvertex.TShape()].name = name;
+      OCCGeometry::global_shape_properties[lastvertex].name = name;
     return shared_from_this();
   }
 
@@ -769,7 +769,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
            trafo.SetTranslation(v);
            BRepBuilderAPI_Transform builder(shape, trafo, true);
            PropagateProperties(builder, shape, occ2ng(trafo));
-           return builder.Shape();
+           return CastShape(builder.Shape());
            // version 2: change location
            // ...
          }, py::arg("v"), "copy shape, and translate copy by vector 'v'")
@@ -822,37 +822,37 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
     .def("bc", [](const TopoDS_Shape & shape, const string & name)
          {
            for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next())
-             OCCGeometry::global_shape_properties[e.Current().TShape()].name = name;
+             OCCGeometry::global_shape_properties[e.Current()].name = name;
            return shape;
          }, py::arg("name"), "sets 'name' property for all faces of shape")
 
     .def("mat", [](const TopoDS_Shape & shape, const string & name)
          {
            for (TopExp_Explorer e(shape, TopAbs_SOLID); e.More(); e.Next())
-             OCCGeometry::global_shape_properties[e.Current().TShape()].name = name;
+             OCCGeometry::global_shape_properties[e.Current()].name = name;
            return shape;
          }, py::arg("name"), "sets 'name' property to all solids of shape")
     
     .def_property("name", [](const TopoDS_Shape & self) -> optional<string> {
-        if (auto name = OCCGeometry::global_shape_properties[self.TShape()].name)
+        if (auto name = OCCGeometry::global_shape_properties[self].name)
           return *name;
         else
           return nullopt;
       }, [](const TopoDS_Shape & self, optional<string> name) {
-        OCCGeometry::global_shape_properties[self.TShape()].name = name;            
+        OCCGeometry::global_shape_properties[self].name = name;
       }, "'name' of shape")
     
     .def_property("maxh",
                   [](const TopoDS_Shape& self)
                   {
-                    return OCCGeometry::global_shape_properties[self.TShape()].maxh;
+                    return OCCGeometry::global_shape_properties[self].maxh;
                   },
                   [](TopoDS_Shape& self, double val)
                   {
                     for (auto typ : { TopAbs_SOLID, TopAbs_FACE,  TopAbs_EDGE, TopAbs_VERTEX })
                       for (TopExp_Explorer e(self, typ); e.More(); e.Next())
                       {
-                        auto & maxh = OCCGeometry::global_shape_properties[e.Current().TShape()].maxh;
+                        auto & maxh = OCCGeometry::global_shape_properties[e.Current()].maxh;
                         maxh = min2(val, maxh);
                       }
                   }, "maximal mesh-size for shape")
@@ -860,16 +860,16 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
     .def_property("hpref",
                   [](const TopoDS_Shape& self)
                   {
-                    return OCCGeometry::global_shape_properties[self.TShape()].hpref;
+                    return OCCGeometry::global_shape_properties[self].hpref;
                   },
                   [](TopoDS_Shape& self, double val)
                   {
-                    auto & hpref = OCCGeometry::global_shape_properties[self.TShape()].hpref;
+                    auto & hpref = OCCGeometry::global_shape_properties[self].hpref;
                     hpref = max2(val, hpref);
                   }, "number of refinement levels for geometric refinement")
     
     .def_property("col", [](const TopoDS_Shape & self) {
-        auto it = OCCGeometry::global_shape_properties.find(self.TShape());
+        auto it = OCCGeometry::global_shape_properties.find(self);
         Vec<4> col(0.2, 0.2, 0.2);
         if (it != OCCGeometry::global_shape_properties.end() && it->second.col)
           col = *it->second.col;
@@ -878,7 +878,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
         Vec<4> col(c[0], c[1], c[2], 1.0);
         if(c.size() == 4)
           col[3] = c[3];
-        OCCGeometry::global_shape_properties[self.TShape()].col = col;    
+        OCCGeometry::global_shape_properties[self].col = col;
       }, "color of shape as RGB - tuple")
     .def("UnifySameDomain", [](const TopoDS_Shape& shape,
                                bool edges, bool faces,
@@ -891,9 +891,9 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
       for (auto typ : { TopAbs_SOLID, TopAbs_FACE,  TopAbs_EDGE })
         for (TopExp_Explorer e(shape, typ); e.More(); e.Next())
           {
-            auto prop = OCCGeometry::global_shape_properties[e.Current().TShape()];
+            auto prop = OCCGeometry::global_shape_properties[e.Current()];
             for (auto mods : history->Modified(e.Current()))
-              OCCGeometry::global_shape_properties[mods.TShape()].Merge(prop);
+              OCCGeometry::global_shape_properties[mods].Merge(prop);
           }
       return unify.Shape();
     }, py::arg("unifyEdges")=true, py::arg("unifyFaces")=true,
@@ -919,9 +919,9 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
           for (auto & s : { shape1, shape2 })
             for (TopExp_Explorer e(s, typ); e.More(); e.Next())
               {
-                auto prop = OCCGeometry::global_shape_properties[e.Current().TShape()];
+                auto prop = OCCGeometry::global_shape_properties[e.Current()];
                 for (auto mods : history->Modified(e.Current()))
-                  OCCGeometry::global_shape_properties[mods.TShape()].Merge(prop);
+                  OCCGeometry::global_shape_properties[mods].Merge(prop);
               }
 #endif        
         */
@@ -945,9 +945,9 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
             for (auto typ : { TopAbs_SOLID, TopAbs_FACE,  TopAbs_EDGE })
               for (TopExp_Explorer e(fused, typ); e.More(); e.Next())
                 {
-                  auto prop = OCCGeometry::global_shape_properties[e.Current().TShape()];
+                  auto prop = OCCGeometry::global_shape_properties[e.Current()];
                   for (auto mods : history->Modified(e.Current()))
-                    OCCGeometry::global_shape_properties[mods.TShape()].Merge(prop);
+                    OCCGeometry::global_shape_properties[mods].Merge(prop);
                 }
             // #endif        
             // PropagateProperties (unify, fused);
@@ -971,9 +971,9 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
           for (auto & s : { shape1, shape2 })
             for (TopExp_Explorer e(s, typ); e.More(); e.Next())
               {
-                auto prop = OCCGeometry::global_shape_properties[e.Current().TShape()];
+                auto prop = OCCGeometry::global_shape_properties[e.Current()];
                 for (auto mods : history->Modified(e.Current()))
-                  OCCGeometry::global_shape_properties[mods.TShape()].Merge(prop);
+                  OCCGeometry::global_shape_properties[mods].Merge(prop);
               }
 #endif // OCC_HAVE_HISTORY
         */
@@ -994,9 +994,9 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
           for (auto & s : { shape1, shape2 })
             for (TopExp_Explorer e(s, typ); e.More(); e.Next())
               {
-                auto prop = OCCGeometry::global_shape_properties[e.Current().TShape()];
+                auto prop = OCCGeometry::global_shape_properties[e.Current()];
                 for (auto mods : history->Modified(e.Current()))
-                  OCCGeometry::global_shape_properties[mods.TShape()].Merge(prop);
+                  OCCGeometry::global_shape_properties[mods].Merge(prop);
               }
 #endif // OCC_HAVE_HISTORY
         */
@@ -1005,6 +1005,12 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
         
         return builder.Shape();        
       }, "cut of shapes")
+    .def("__eq__", [] (const TopoDS_Shape& shape1, const TopoDS_Shape& shape2) {
+      return shape1.IsSame(shape2);
+    })
+    .def("__hash__", [] (const TopoDS_Shape& shape) {
+      return shape.HashCode(std::numeric_limits<Standard_Integer>::max());
+    })
 
     .def("Reversed", [](const TopoDS_Shape & shape) {
         return CastShape(shape.Reversed()); })
@@ -1029,9 +1035,9 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
             for (auto typ : { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX })
               for (TopExp_Explorer e(shape, typ); e.More(); e.Next())
                 {
-                  auto prop = OCCGeometry::global_shape_properties[e.Current().TShape()];
+                  auto prop = OCCGeometry::global_shape_properties[e.Current()];
                   for (auto mods : builder.Generated(e.Current()))
-                    OCCGeometry::global_shape_properties[mods.TShape()].Merge(prop);
+                    OCCGeometry::global_shape_properties[mods].Merge(prop);
                 }
             
             return builder.Shape();
@@ -1052,9 +1058,9 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
           for (auto typ : { TopAbs_EDGE, TopAbs_VERTEX })
             for (TopExp_Explorer e(shape, typ); e.More(); e.Next())
               {
-                auto prop = OCCGeometry::global_shape_properties[e.Current().TShape()];
+                auto prop = OCCGeometry::global_shape_properties[e.Current()];
                 for (auto mods : builder.Generated(e.Current()))
-                  OCCGeometry::global_shape_properties[mods.TShape()].Merge(prop);
+                  OCCGeometry::global_shape_properties[mods].Merge(prop);
               }
 
           return builder.Shape();          
@@ -1070,7 +1076,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
         PropagateProperties (mkFillet, shape);
         for (auto e : edges)
           for (auto gen : mkFillet.Generated(e))
-            OCCGeometry::global_shape_properties[gen.TShape()].name = "fillet";
+            OCCGeometry::global_shape_properties[gen].name = "fillet";
         return mkFillet.Shape();
       }, py::arg("edges"), py::arg("r"), "make fillets for edges 'edges' of radius 'r'")
   
@@ -1083,7 +1089,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
         PropagateProperties (mkChamfer, shape);
         for (auto e : edges)
           for (auto gen : mkChamfer.Generated(e))
-            OCCGeometry::global_shape_properties[gen.TShape()].name = "chamfer";
+            OCCGeometry::global_shape_properties[gen].name = "chamfer";
         return mkChamfer.Shape();
 #else
         throw Exception("MakeChamfer not available for occ-version < 7.4");
@@ -1177,17 +1183,21 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
 
            std::vector<double> p[3];
            std::vector<double> n[3];
-           py::list names, colors;
+           py::list names, colors, solid_names;
+           std::vector<std::vector<int>> solid_face_map;
 
            int index = 0;
 
            Box<3> box(Box<3>::EMPTY_BOX);
+           TopTools_IndexedMapOfShape fmap;
            for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next())
            {
                TopoDS_Face face = TopoDS::Face(e.Current());
+               if(fmap.Contains(face)) continue;
                // Handle(TopoDS_Face) face = e.Current();
+               fmap.Add(face);
                ExtractFaceData(face, index, p, n, box);
-               auto & props = OCCGeometry::global_shape_properties[face.TShape()];
+               auto & props = OCCGeometry::global_shape_properties[face];
                if(props.col)
                {
                  auto & c = *props.col;
@@ -1204,6 +1214,19 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
                index++;
            }
 
+           for(auto& solid : GetSolids(shape))
+             {
+               std::vector<int> faces;
+               for(auto& face : GetFaces(solid))
+                 faces.push_back(fmap.FindIndex(face)-1);
+               solid_face_map.push_back(move(faces));
+               auto& props = OCCGeometry::global_shape_properties[solid];
+               if(props.name)
+                 solid_names.append(*props.name);
+               else
+                 solid_names.append("");
+             }
+
            std::vector<double> edge_p[2];
            py::list edge_names, edge_colors;
            index = 0;
@@ -1211,7 +1234,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
            {
                TopoDS_Edge edge = TopoDS::Edge(e.Current());
                ExtractEdgeData(edge, index, edge_p, box);
-               auto & props = OCCGeometry::global_shape_properties[edge.TShape()];
+               auto & props = OCCGeometry::global_shape_properties[edge];
                if(props.col)
                {
                  auto & c = *props.col;
@@ -1263,6 +1286,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
            data["autoscale"] = false;
            data["colors"] = colors;
            data["names"] = names;
+           data["solid_names"] = solid_names;
 
            py::list edges;
            edges.append(edge_p[0]);
@@ -1270,6 +1294,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
            data["edges"] = edges;
            data["edge_names"] = edge_names;
            data["edge_colors"] = edge_colors;
+           data["solid_face_map"] = solid_face_map;
            return data;
          })
     ;
@@ -1437,11 +1462,11 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
         }))
     .def_property("quad_dominated", [](const TopoDS_Face& self) -> optional<bool>
                   {
-                    return OCCGeometry::global_shape_properties[self.TShape()].quad_dominated;
+                    return OCCGeometry::global_shape_properties[self].quad_dominated;
                   },
                   [](TopoDS_Face& self, optional<bool> quad_dominated)
                   {
-                    OCCGeometry::global_shape_properties[self.TShape()].quad_dominated = quad_dominated;
+                    OCCGeometry::global_shape_properties[self].quad_dominated = quad_dominated;
                   })
     .def_property_readonly("surf", [] (TopoDS_Face face) -> Handle(Geom_Surface)
          {
@@ -1492,13 +1517,13 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
             {
                 auto & props = OCCGeometry::global_shape_properties;
                 for(auto & s : GetSolids(shapes[i]))
-                    props[s.TShape()].layer = i+1;
+                    props[s].layer = i+1;
                 for(auto & s : GetFaces(shapes[i]))
-                    props[s.TShape()].layer = i+1;
+                    props[s].layer = i+1;
                 for(auto & s : GetEdges(shapes[i]))
-                    props[s.TShape()].layer = i+1;
+                    props[s].layer = i+1;
                 for(auto & s : GetVertices(shapes[i]))
-                    props[s.TShape()].layer = i+1;
+                    props[s].layer = i+1;
             }
           }
 
@@ -1544,6 +1569,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
   };
   
   py::class_<ListOfShapes> (m, "ListOfShapes")
+    .def(py::init<vector<TopoDS_Shape>>())
     .def("__iter__", [](ListOfShapes &s) {
         return py::make_iterator(ListOfShapesIterator(&*s.begin()),
                                  ListOfShapesIterator(&*s.end()));
@@ -1579,7 +1605,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
            ListOfShapes selected;
            std::regex pattern(name);
            for (auto s : self)
-             if (auto sname = OCCGeometry::global_shape_properties[s.TShape()].name)
+             if (auto sname = OCCGeometry::global_shape_properties[s].name)
                if (std::regex_match(*sname, pattern))
                  selected.push_back(s);
            return selected;
@@ -1602,7 +1628,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
 
     .def("Sorted",[](ListOfShapes self, gp_Vec dir)
          {
-           std::map<Handle(TopoDS_TShape), double> sortval;
+           std::map<TopoDS_Shape, double> sortval;
            for (auto shape : self)
              {
                GProp_GProps props;
@@ -1622,12 +1648,12 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
                  }
                
                double val = center.X()*dir.X() + center.Y()*dir.Y() + center.Z() * dir.Z();
-               sortval[shape.TShape()] = val;
+               sortval[shape] = val;
              }
 
            std::sort (std::begin(self), std::end(self),
-                      [&](TopoDS_Shape a, TopoDS_Shape b)
+                      [&](const TopoDS_Shape& a, const TopoDS_Shape& b)
-                      { return sortval[a.TShape()] < sortval[b.TShape()]; });
+                      { return sortval[a] < sortval[b]; });
            return self;
          }, py::arg("dir"), "returns list of shapes, where center of gravity is sorted in direction of 'dir'")
     
@@ -1654,7 +1680,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
       {
         for(auto& shape : shapes)
           {
-            OCCGeometry::global_shape_properties[shape.TShape()].name = name;
+            OCCGeometry::global_shape_properties[shape].name = name;
           }
       }, "set name for all elements of list")
     .def_property("col", [](ListOfShapes& shapes) {
@@ -1664,7 +1690,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
         if(c.size() == 4)
           col[3] = c[3];
         for(auto& shape : shapes)
-          OCCGeometry::global_shape_properties[shape.TShape()].col = col;
+          OCCGeometry::global_shape_properties[shape].col = col;
       }, "set col for all elements of list")
     
     .def_property("maxh", [](ListOfShapes& shapes)
@@ -1676,13 +1702,13 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
         for(auto& shape : shapes)
           {
             for(auto& s : GetSolids(shape))
-              OCCGeometry::global_shape_properties[s.TShape()].maxh = maxh;
+              OCCGeometry::global_shape_properties[s].maxh = maxh;
             for(auto& s : GetFaces(shape))
-              OCCGeometry::global_shape_properties[s.TShape()].maxh = maxh;
+              OCCGeometry::global_shape_properties[s].maxh = maxh;
             for(auto& s : GetEdges(shape))
-              OCCGeometry::global_shape_properties[s.TShape()].maxh = maxh;
+              OCCGeometry::global_shape_properties[s].maxh = maxh;
             for(auto& s : GetVertices(shape))
-              OCCGeometry::global_shape_properties[s.TShape()].maxh = maxh;
+              OCCGeometry::global_shape_properties[s].maxh = maxh;
           }
       }, "set maxh for all elements of list")
     .def_property("hpref", [](ListOfShapes& shapes)
@@ -1693,7 +1719,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
       {
         for(auto& shape : shapes)
           {
-            auto& val = OCCGeometry::global_shape_properties[shape.TShape()].hpref;
+            auto& val = OCCGeometry::global_shape_properties[shape].hpref;
             val = max2(hpref, val);
           }
       }, "set hpref for all elements of list")
@@ -1704,7 +1730,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
                   [](ListOfShapes& shapes, optional<bool> quad_dominated)
                   {
                     for(auto& shape : shapes)
-                      OCCGeometry::global_shape_properties[shape.TShape()].quad_dominated = quad_dominated;
+                      OCCGeometry::global_shape_properties[shape].quad_dominated = quad_dominated;
                   })
     
     .def("Identify", [](const ListOfShapes& me,
@@ -1802,7 +1828,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
                         optional<string> bot, optional<string> top, optional<string> mantle) {
     auto builder = BRepPrimAPI_MakeCylinder (gp_Ax2(cpnt, cdir), r, h);
     if(mantle)
-      OCCGeometry::global_shape_properties[builder.Face().TShape()].name = *mantle;
+      OCCGeometry::global_shape_properties[builder.Face()].name = *mantle;
     auto pyshape = py::cast(builder.Solid());
     gp_Vec v = cdir;
     if(bot)
@@ -2009,6 +2035,8 @@ tangents : Dict[int, gp_Vec2d]
   
   m.def("Glue", [] (const std::vector<TopoDS_Shape> shapes) -> TopoDS_Shape
         {
+          if(shapes.size() == 1)
+            return shapes[0];
           BOPAlgo_Builder builder;
           for (auto & s : shapes)
             {
@@ -2053,9 +2081,9 @@ tangents : Dict[int, gp_Vec2d]
             for (auto & s : shapes)
               for (TopExp_Explorer e(s, typ); e.More(); e.Next())
                 {
-                  auto prop = OCCGeometry::global_shape_properties[e.Current().TShape()];
+                  auto prop = OCCGeometry::global_shape_properties[e.Current()];
                   for (auto mods : history->Modified(e.Current()))
-                    OCCGeometry::global_shape_properties[mods.TShape()].Merge(prop);
+                    OCCGeometry::global_shape_properties[mods].Merge(prop);
                 }
 #endif // OCC_HAVE_HISTORY
           */
@@ -2084,9 +2112,9 @@ tangents : Dict[int, gp_Vec2d]
 
           for (TopExp_Explorer e(shape, TopAbs_SOLID); e.More(); e.Next())
             {
-              auto prop = OCCGeometry::global_shape_properties[e.Current().TShape()];
+              auto prop = OCCGeometry::global_shape_properties[e.Current()];
               for (auto mods : history->Modified(e.Current()))
-                OCCGeometry::global_shape_properties[mods.TShape()].Merge(prop);
+                OCCGeometry::global_shape_properties[mods].Merge(prop);
             }
 #endif // OCC_HAVE_HISTORY
           */

libsrc/occ/vsocc.cpp

@@ -548,7 +548,7 @@ namespace netgen
 
          if (!occgeometry->fvispar[i-1].IsHighlighted())
          {
-            auto c = OCCGeometry::global_shape_properties[face.TShape()].col.value_or(Vec<4>(0,1,0,1) );
+            auto c = OCCGeometry::global_shape_properties[face].col.value_or(Vec<4>(0,1,0,1) );
             for(auto j : Range(4))
                mat_col[j] = c[j];
          }

libsrc/stlgeom/meshstlsurface.cpp

@@ -106,7 +106,7 @@ static void STLFindEdges (STLGeometry & geom, Mesh & mesh,
 		   << ", trig2 = " << trig2
 		   << ", trig2b = " << trig2b << endl;
 
-	  if (trig1 <= 0 || trig2 <= 0 || trig1b <= 0 || trig2b <= 0)
+	  if (trig1 <= 0 || trig2 < 0 || trig1b <= 0 || trig2b < 0)
 	    {
 	      cout << "negative trigs, "
 		   << ", trig1 = " << trig1
@@ -177,10 +177,13 @@ static void STLFindEdges (STLGeometry & geom, Mesh & mesh,
 	  mesh.AddSegment (seg);
 
 
+          if(trig2 != 0)
+            {
 	  Segment seg2;
 	  seg2[0] = p2 + PointIndex::BASE-1;;
 	  seg2[1] = p1 + PointIndex::BASE-1;;
 	  seg2.si = geom.GetTriangle(trig2).GetFaceNum();
+
 	  seg2.edgenr = i;
 
 	  seg2.epgeominfo[0].edgenr = i;
@@ -219,8 +222,8 @@ static void STLFindEdges (STLGeometry & geom, Mesh & mesh,
 	      (*testout) << "Get GeomInfo PROBLEM" << endl;
 	    }
 	  */	  
-
 	  mesh.AddSegment (seg2);
+            }
 	}
     }
 

libsrc/stlgeom/python_stl.cpp

@@ -9,7 +9,7 @@
 using namespace netgen;
 namespace netgen
 {
-  //extern shared_ptr<Mesh> mesh;
+  extern shared_ptr<Mesh> mesh;
   extern shared_ptr<NetgenGeometry> ng_geometry;
 }
 
@@ -125,11 +125,12 @@ NGCORE_API_EXPORT void ExportSTL(py::module & m)
 {
   py::class_<STLGeometry,shared_ptr<STLGeometry>, NetgenGeometry> (m,"STLGeometry")
     .def(py::init<>())
-    .def(py::init<>([](const string& filename)
+    .def(py::init<>([](const string& filename, bool surface)
                     {
                       ifstream ist(filename);
-                      return shared_ptr<STLGeometry>(STLGeometry::Load(ist));
+                      return shared_ptr<STLGeometry>(STLGeometry::Load(ist,
-                    }), py::arg("filename"),
+                                                                       surface));
+                    }), py::arg("filename"), py::arg("surface")=false,
       py::call_guard<py::gil_scoped_release>())
     .def(NGSPickle<STLGeometry>())
     .def("_visualizationData", [](shared_ptr<STLGeometry> stl_geo)
@@ -182,7 +183,8 @@ NGCORE_API_EXPORT void ExportSTL(py::module & m)
             return res;
          }, py::call_guard<py::gil_scoped_release>())
     .def("GenerateMesh", [] (shared_ptr<STLGeometry> geo,
-                             MeshingParameters* pars, py::kwargs kwargs)
+                             MeshingParameters* pars,
+                             shared_ptr<Mesh> mesh, py::kwargs kwargs)
                          {
                            MeshingParameters mp;
                            STLParameters stlparam;
@@ -197,16 +199,22 @@ NGCORE_API_EXPORT void ExportSTL(py::module & m)
                              CreateSTLParametersFromKwargs(stlparam, kwargs);
                              CreateMPfromKwargs(mp, kwargs); // this will throw if any kwargs are not passed
                            }
-                           auto mesh = make_shared<Mesh>();
+                           if(!mesh)
+                             {
+                               mesh = make_shared<Mesh>();
+                             }
                            mesh->SetGeometry(geo);
                            ng_geometry = geo;
                            SetGlobalMesh(mesh);
                            auto result = STLMeshingDummy(geo.get(), mesh, mp, stlparam);
                            if(result != 0)
-                             throw Exception("Meshing failed!");
+                             {
+                               netgen::mesh = mesh;
+                               throw Exception("Meshing failed!");
+                             }
 
                            return mesh;
-                         }, py::arg("mp") = nullptr,
+                         }, py::arg("mp") = nullptr, py::arg("mesh") = nullptr,
       py::call_guard<py::gil_scoped_release>(),
          (meshingparameter_description + stlparameter_description).c_str())
     .def("Draw", FunctionPointer

libsrc/stlgeom/stlgeom.cpp

@@ -2578,7 +2578,7 @@ void STLGeometry :: CalcEdgeDataAngles()
   for (int i = 1; i <= GetNTE(); i++)
     {
       STLTopEdge & edge = GetTopEdge (i);
-      double cosang = 
+      double cosang = edge.TrigNum(2) == 0 ? 1. :
 	GetTriangle(edge.TrigNum(1)).Normal() *
 	GetTriangle(edge.TrigNum(2)).Normal();
       edge.SetCosAngle (cosang);
@@ -2611,6 +2611,8 @@ void STLGeometry :: FindEdgesFromAngles(const STLParameters& stlparam)
   for (int i = 1; i <= edgedata->Size(); i++)
     {
       STLTopEdge & sed = edgedata->Elem(i);
+      if(sed.TrigNum(2) == 0)
+        sed.SetStatus(ED_CONFIRMED);
       if (sed.GetStatus() == ED_CANDIDATE || 
 	  sed.GetStatus() == ED_UNDEFINED)
 	{
@@ -3187,7 +3189,7 @@ void STLGeometry :: BuildSmoothEdges ()
       ng1 = trig.GeomNormal(points);
       ng1 /= (ng1.Length() + 1e-24);
 
-      for (int j = 1; j <= 3; j++)
+      for (int j = 1; j <= NONeighbourTrigs(i); j++)
 	{ 
 	  int nbt = NeighbourTrig (i, j);
 
@@ -3261,7 +3263,7 @@ void STLGeometry :: AddConeAndSpiralEdges(const STLParameters& stlparam)
 	  STLTrigId t = chart.GetChartTrig1(j); 
 	  const STLTriangle& tt = GetTriangle(t);
 
-	  for (int k = 1; k <= 3; k++)
+	  for (int k = 1; k <= NONeighbourTrigs(t); k++)
 	    {
 	      STLTrigId nt = NeighbourTrig(t,k); 
 	      if (GetChartNr(nt) != i && !TrigIsInOC(nt,i))

libsrc/stlgeom/stlgeomchart.cpp

@@ -238,7 +238,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh, const MeshingParameters& mparam, cons
 			      */
 			      
 			      //find overlapping charts exacter (fast, too): 
-			      for (int k = 1; k <= 3; k++) 
+			      for (int k = 1; k <= NONeighbourTrigs(nt); k++)
 				{ 
 				  int nnt = NeighbourTrig(nt,k);
 				  if (GetMarker(nnt) != chartnum)
@@ -387,7 +387,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh, const MeshingParameters& mparam, cons
                       // NgProfiler::StartTimer (timer4a);
                       
 		      if (spiralcheckon && !isdirtytrig)
-			for (int k = 1; k <= 3; k++) 
+			for (int k = 1; k <= NONeighbourTrigs(nt); k++)
 			  {
                             // NgProfiler::StartTimer (timer4b);                            
 			    STLTrigId nnt = NeighbourTrig(nt,k);
@@ -695,7 +695,7 @@ void STLGeometry :: GetInnerChartLimes(NgArray<twoint>& limes, ChartId chartnum)
     {
       STLTrigId t = chart.GetChartTrig1(j); 
       const STLTriangle& tt = GetTriangle(t);
-      for (int k = 1; k <= 3; k++)
+      for (int k = 1; k <= NONeighbourTrigs(t); k++)
 	{
 	  STLTrigId nt = NeighbourTrig(t,k); 
 	  if (GetChartNr(nt) != chartnum)
@@ -756,7 +756,7 @@ void STLGeometry :: GetDirtyChartTrigs(int chartnum, STLChart& chart,
       STLTrigId t = chart.GetChartTrig1(j); 
       const STLTriangle& tt = GetTriangle(t);
       
-      for (int k = 1; k <= 3; k++)
+      for (int k = 1; k <= NONeighbourTrigs(t); k++)
 	{
 	  STLTrigId nt = NeighbourTrig(t,k); 
 	  if (GetChartNr(nt) != chartnum && outercharttrigs[nt] != chartnum)

libsrc/stlgeom/stlgeommesh.cpp

@@ -1169,7 +1169,7 @@ void STLGeometry :: RestrictHChartDistOneChart(ChartId chartnum, NgArray<int>& a
     {
       int t = chart.GetChartTrig1(j); 
       tt = GetTriangle(t);
-      for (int k = 1; k <= 3; k++)
+      for (int k = 1; k <= NONeighbourTrigs(t); k++)
 	{
 	  int nt = NeighbourTrig(t,k); 
 	  if (GetChartNr(nt) != chartnum)
@@ -1495,6 +1495,9 @@ int STLMeshingDummy (STLGeometry* stlgeometry, shared_ptr<Mesh> & mesh, const Me
       if (multithread.terminate)
 	return 0;
 
+      if(stlgeometry->IsSurfaceSTL())
+        return 0;
+
       if (mparam.perfstepsstart <= MESHCONST_MESHVOLUME && 
 	  mparam.perfstepsend >= MESHCONST_MESHVOLUME)
 	{

libsrc/stlgeom/stltopology.cpp

@@ -338,7 +338,7 @@ void STLTopology :: Save (const filesystem::path & filename) const
 }
 
 
-STLGeometry *  STLTopology ::Load (istream & ist)
+STLGeometry *  STLTopology ::Load (istream & ist, bool surface)
 {
   // Check if the file starts with "solid". If not, the file is binary
   {
@@ -457,6 +457,7 @@ STLGeometry *  STLTopology ::Load (istream & ist)
       PrintWarning("File has normal vectors which differ extremely from geometry->correct with stldoctor!!!");
     }
 
+  geom->surface = surface;
   geom->InitSTLGeometry(readtrigs);
   return geom;
 }
@@ -650,6 +651,7 @@ void STLTopology :: FindNeighbourTrigs()
 	  }
       }
 
+  if(!surface)
   for (int i = 1; i <= ne; i++)
     {
       const STLTopEdge & edge = GetTopEdge (i);
@@ -668,9 +670,12 @@ void STLTopology :: FindNeighbourTrigs()
 	  const STLTriangle & t = GetTriangle (i);
 	  for (int j = 1; j <= 3; j++)
 	    {
-	      const STLTriangle & nbt = GetTriangle (t.NBTrigNum(j));
+              if(t.NBTrigNum(j) != 0)
-	      if (!t.IsNeighbourFrom (nbt))
+                {
-		orientation_ok = 0;
+                  const STLTriangle & nbt = GetTriangle (t.NBTrigNum(j));
+                  if (!t.IsNeighbourFrom (nbt))
+                    orientation_ok = 0;
+                }
 	    }
 	}
     }
@@ -801,7 +806,8 @@ void STLTopology :: FindNeighbourTrigs()
       neighbourtrigs.SetSize(GetNT());
       for (int i = 1; i <= GetNT(); i++)
 	for (int k = 1; k <= 3; k++)
-	  AddNeighbourTrig (i, GetTriangle(i).NBTrigNum(k));
+          if(GetTriangle(i).NBTrigNum(k) != 0)
+            AddNeighbourTrig (i, GetTriangle(i).NBTrigNum(k));
     }
   else
     {

libsrc/stlgeom/stltopology.hpp

@@ -120,7 +120,12 @@ public:
 
 
   STLTriangle (const STLPointId * apts);
-  STLTriangle () {pts[0]=0;pts[1]=0;pts[2]=0;}
+  STLTriangle ()
+  {
+    pts[0]=0;pts[1]=0;pts[2]=0;
+    nbtrigs[0][0] = nbtrigs[0][1] = nbtrigs[0][2] = 0.;
+    nbtrigs[1][0] = nbtrigs[1][1] = nbtrigs[1][2] = 0.;
+  }
 
   void DoArchive(Archive& ar)
   {
@@ -282,6 +287,7 @@ protected:
   Array<STLTriangle, STLTrigId> trias;
   NgArray<STLTopEdge> topedges;
   Array<Point<3>, STLPointId> points;
+  bool surface = false;
 
   // mapping of sorted pair of points to topedge
   INDEX_2_HASHTABLE<int> * ht_topedges;
@@ -313,13 +319,15 @@ public:
   virtual ~STLTopology();
 
   static STLGeometry * LoadNaomi (istream & ist);
-  static STLGeometry * Load (istream & ist);
+  DLL_HEADER static STLGeometry * Load (istream & ist, bool surface=false);
   static STLGeometry * LoadBinary (istream & ist);
 
   void Save (const filesystem::path & filename) const;
   void SaveBinary (const filesystem::path & filename, const char* aname) const;
   void SaveSTLE (const filesystem::path & filename) const; // stores trigs and edges
 
+  bool IsSurfaceSTL() const { return surface; }
+
   virtual void DoArchive(Archive& ar)
   {
     ar & trias & points & boundingbox & pointtol;

libsrc/visualization/vsfieldlines.cpp

@@ -260,7 +260,7 @@ namespace netgen
    
     double phaser=1.0;
     double phasei=0.0;
-    std::function eval_func = [&](int elnr, const double * lami, Vec<3> & vec)
+    std::function<bool(int, const double *, Vec<3> &)> eval_func = [&](int elnr, const double * lami, Vec<3> & vec)
     {
         double values[6] = {0., 0., 0., 0., 0., 0.};
         bool drawelem;

libsrc/visualization/vsmesh.cpp

@@ -597,8 +597,8 @@ namespace netgen
 
 	for (int i = 1; i <= mesh->GetNE(); i++)
 	  {
-            if (mesh->VolumeElement(i).flags.badel ||
+            if (mesh->VolumeElement(i).Flags().badel ||
-		mesh->VolumeElement(i).flags.illegal ||
+		mesh->VolumeElement(i).Flags().illegal ||
 		(i == vispar.drawelement))
 	      {
 		// copy to be thread-safe
@@ -636,7 +636,7 @@ namespace netgen
 
 	for (ElementIndex ei : mesh->VolumeElements().Range())
 	  {
-            if ((*mesh)[ei].flags.badel)
+            if ((*mesh)[ei].Flags().badel)
 	      {
 		// copy to be thread-safe
 		Element el = (*mesh)[ei];

libsrc/visualization/vssolution.cpp

@@ -795,12 +795,24 @@ namespace netgen
             static double oldrad = 0;
           
             mesh->GetBox (pmin, pmax, -1);
-            center = Center (pmin, pmax);
+            if(vispar.use_center_coords && zoomall == 2)
+              {
+                center.X() = vispar.centerx;
+                center.Y() = vispar.centery;
+                center.Z() = vispar.centerz;
+              }
+            else if(selpoint >= 1 && zoomall == 2)
+              center = mesh->Point(selpoint);
+            else if(vispar.centerpoint >= 1 && zoomall == 2)
+              center = mesh->Point(vispar.centerpoint);
+            else
+              center = Center (pmin, pmax);
             rad = 0.5 * Dist (pmin, pmax);
+            if(rad == 0) rad = 1e-6;
           
             glEnable (GL_NORMALIZE);
           
-            if (rad > 1.5 * oldrad ||
+            if (rad > 1.2 * oldrad ||
                 mesh->GetMajorTimeStamp() > surfeltimestamp ||
                 zoomall)
               {

ng/occgeom.tcl

@@ -101,7 +101,7 @@ proc occdialogbuildtree {} {
 	set nrfaces [expr [llength $faces]]
 	if {$nrfaces >= 2} {
 	    #$hlist add ErrorFaces -itemtype text -text "Faces with surface meshing error"
-        $w.tree insert {} -id ErrorFaces -text "Faces with surface meshing error"
+        $w.tree insert {} end -id "ErrorFaces" -text "Faces with surface meshing error"
 	    #$w.mtre open ErrorFaces
         $w.tree item ErrorFaces -open true
 	    set i [expr 0]
@@ -109,12 +109,12 @@ proc occdialogbuildtree {} {
 		set entity [lindex $faces [expr $i]]
         set myroot [string range $entity 0 [string last / $entity]-1]
         if { [string length $myroot] == 0 } {
-            set myroot ErrorFaces
+            set myroot "ErrorFaces"
-        }
+                }
 		incr i 1
 		set entityname [lindex $faces [expr $i]]
 		#$hlist add ErrorFaces/$entity -text $entityname -data $entityname
-        $w.tree insert {myroot} end -id $entity -text $entityname -value 0
+        $w.tree insert $myroot end -id $entity -text $entityname -value 0
 		incr i 1
 	    }
 	}

ng/onetcl.cpp

@@ -3992,18 +3992,18 @@ DLL_HEADER const char * ngscript[] = {""
 ,"set faces [Ng_OCCCommand getunmeshedfaceinfo]\n"
 ,"set nrfaces [expr [llength $faces]]\n"
 ,"if {$nrfaces >= 2} {\n"
-,"$w.tree insert {} -id ErrorFaces -text \"Faces with surface meshing error\"\n"
+,"$w.tree insert {} end -id \"ErrorFaces\" -text \"Faces with surface meshing error\"\n"
 ,"$w.tree item ErrorFaces -open true\n"
 ,"set i [expr 0]\n"
 ,"while {$i < $nrfaces} {\n"
 ,"set entity [lindex $faces [expr $i]]\n"
 ,"set myroot [string range $entity 0 [string last / $entity]-1]\n"
 ,"if { [string length $myroot] == 0 } {\n"
-,"set myroot ErrorFaces\n"
+,"set myroot \"ErrorFaces\"\n"
 ,"}\n"
 ,"incr i 1\n"
 ,"set entityname [lindex $faces [expr $i]]\n"
-,"$w.tree insert {myroot} end -id $entity -text $entityname -value 0\n"
+,"$w.tree insert $myroot end -id $entity -text $entityname -value 0\n"
 ,"incr i 1\n"
 ,"}\n"
 ,"}\n"

python/CMakeLists.txt

@@ -32,7 +32,7 @@ if(pybind11_stubgen AND NOT ${pybind11_stubgen} EQUAL 0)
 for better autocompletion support install it with pip.")
 else()
   message("-- Found pybind11-stubgen: ${stubgen_path}")
-  install(CODE "execute_process(COMMAND ${PYTHON_EXECUTABLE} -m pybind11_stubgen --no-setup-py netgen)")
+  install(CODE "execute_process(COMMAND ${PYTHON_EXECUTABLE} -m pybind11_stubgen --no-setup-py --ignore-invalid=all netgen)")
   install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/../stubs/netgen-stubs/ DESTINATION ${NG_INSTALL_DIR_PYTHON}/netgen/ COMPONENT netgen)
 endif()
 endif(BUILD_STUB_FILES)

python/__init__.py

@@ -6,7 +6,8 @@ _netgen_bin_dir=os.path.realpath(os.path.join(os.path.dirname(__file__),'..',con
 _netgen_lib_dir=os.path.realpath(os.path.join(os.path.dirname(__file__),'..',config.NETGEN_PYTHON_RPATH))
 
 if sys.platform.startswith('win'):
-    if sys.version >= '3.8':
+    v = sys.version_info
+    if v.major == 3 and v.minor >= 8:
         os.add_dll_directory(_netgen_bin_dir)
     else:
         os.environ['PATH'] += ';'+_netgen_bin_dir

python/gui.py

@@ -30,6 +30,13 @@ def StartGUI():
 
     win.tk.eval( netgen.libngpy._meshing._ngscript)
 
+    try:
+        from IPython import get_ipython
+        ipython = get_ipython()
+        ipython.magic('gui tk')
+    except:
+        pass
+
     def _Redraw(*args, **kwargs):
         if libngpy._meshing._Redraw(*args, **kwargs):
             import netgen

rules/pyramidrules2.rls

@@ -163,6 +163,69 @@ freeset
 endrule
 
 
+rule "pyramid with one trig, left"
+
+quality 20
+
+mappoints
+(0, 0, 0);
+(1, 0, 0);
+(1, 1, 0);
+(0, 1, 0);
+(0.5, 0.5, -0.5);
+
+mapfaces
+(1, 2, 3, 4) del;
+(3, 2, 5) del;
+
+newpoints
+
+newfaces
+(1, 2, 5);
+(3, 4, 5);
+(4, 1, 5);
+
+elements
+(1, 2, 3, 4, 5);
+
+freezone2
+{ 1 P1 };
+{ 1 P2 };
+{ 1 P3 };
+{ 1 P4 };
+{ 1 P5 };
+{ 0.34 P1, 0.34 P2, 0.34 P5, -0.01 P3 };
+{ 0.34 P3, 0.34 P4, 0.34 P5, -0.02 P1 };
+{ 0.34 P1, 0.34 P4, 0.34 P5, -0.02 P2 };
+
+freezonelimit
+{ 1 P1 };
+{ 1 P2 };
+{ 1 P3 };
+{ 1 P4 };
+{ 1 P5 };
+{ 0.333 P1, 0.333 P2, 0.334 P5, 0 P3 };
+{ 0.333 P3, 0.333 P4, 0.334 P5, 0 P1 };
+{ 0.333 P1, 0.333 P4, 0.334 P5, 0 P2 };
+
+orientations
+(1, 2, 3, 5);
+(1, 3, 4, 5);
+
+
+freeset
+1 2 3 5;
+freeset
+1 3 4 5;
+freeset
+1 2 5 6;
+freeset
+3 4 5 7;
+freeset 
+1 4 5 8;
+endrule
+
+
 
 
 

tests/build_pip.ps1

@@ -9,6 +9,6 @@ $env:NETGEN_CCACHE = 1
 
 $pydir=$args[0]
 & $pydir\python.exe --version
-& $pydir\python.exe -m pip install scikit-build wheel numpy twine
+& $pydir\python.exe -m pip install scikit-build wheel numpy twine pybind11-stubgen
 & $pydir\python setup.py bdist_wheel -G"Visual Studio 16 2019"
 & $pydir\python -m twine upload dist\*.whl

tests/build_pip.sh

@@ -11,15 +11,15 @@ for pyversion in 38 39 310
 do
     export PYDIR="/opt/python/cp${pyversion}-cp${pyversion}/bin"
     echo $PYDIR
-    $PYDIR/pip install -U pytest-check numpy wheel scikit-build
+    $PYDIR/pip install -U pytest-check numpy wheel scikit-build pybind11-stubgen
 
     rm -rf _skbuild
-    $PYDIR/pip wheel --use-feature=in-tree-build .
+    $PYDIR/pip wheel .
     auditwheel repair netgen_mesher*-cp${pyversion}-*.whl
     rm netgen_mesher-*.whl
 
     rm -rf _skbuild
-    NETGEN_ARCH=avx2 $PYDIR/pip wheel --use-feature=in-tree-build .
+    NETGEN_ARCH=avx2 $PYDIR/pip wheel .
     auditwheel repair netgen_mesher_avx2*-cp${pyversion}-*.whl
     rm netgen_mesher_avx2-*.whl
 

tests/build_pip_mac.sh

@@ -6,7 +6,7 @@ export NETGEN_CCACHE=1
 
 export PYDIR=/Library/Frameworks/Python.framework/Versions/$1/bin
 $PYDIR/python3 --version
-$PYDIR/pip3 install --user numpy twine scikit-build wheel
+$PYDIR/pip3 install --user numpy twine scikit-build wheel pybind11-stubgen
 
 export CMAKE_OSX_ARCHITECTURES='arm64;x86_64'
 $PYDIR/python3 setup.py bdist_wheel --plat-name macosx-10.15-universal2 -j10

tests/dockerfile

@@ -1,4 +1,4 @@
-FROM ubuntu:21.10
+FROM ubuntu:22.04
 ENV DEBIAN_FRONTEND=noninteractive
 MAINTAINER Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>
 RUN apt-get update && apt-get -y install \
@@ -29,5 +29,5 @@ RUN apt-get update && apt-get -y install \
         tcl-dev \
         tk-dev
 
-RUN python3 -m pip install pytest-check
+RUN python3 -m pip install pytest-check pybind11-stubgen
 ADD . /root/src/netgen

tests/dockerfile_mpi

@@ -23,5 +23,5 @@ RUN apt-get update && apt-get -y install \
         tcl-dev \
         tk-dev
 
-RUN python3 -m pip install pytest-mpi pytest-check pytest
+RUN python3 -m pip install pytest-mpi pytest-check pytest pybind11-stubgen
 ADD . /root/src/netgen

tests/pytest/test_array.py

@@ -13,5 +13,21 @@ def test_array_numpy():
     a.NumPy().sort()
     assert(all(a == array([0,0,2,2,5])))
 
+def test_mesh_elements_numpy_array_access():
+    from netgen.csg import unit_cube
+    mesh = unit_cube.GenerateMesh()
+    np_els = mesh.Elements3D().NumPy()
+    vol_nodes = np_els["nodes"]
+    indices = np_els["index"]
+    nps = np_els["np"]
+    for nodes, el, index, np in zip(vol_nodes, mesh.Elements3D(), indices, nps):
+        for n1, n2 in zip(nodes, el.vertices):
+            assert n1 == n2
+        for n in nodes[len(el.vertices):]:
+            assert n == 0
+        assert el.index == index
+        assert len(el.vertices) == np
+
 if __name__ == "__main__":
     test_array_numpy()
+    test_mesh_elements_numpy_array_access()

tests/pytest/test_boundarylayer.py

@@ -124,8 +124,9 @@ def test_pyramids(outside):
     assert ngs.Integrate(1, mesh.Materials("layer")) == pytest.approx(0.0016)
     assert ngs.Integrate(1, mesh.Materials("air")) == pytest.approx(0.9664 if outside else 0.968)
 
+# not working yet
 @pytest.mark.parametrize("outside", [True, False])
-def test_with_inner_corner(outside, capfd):
+def _test_with_inner_corner(outside, capfd):
     geo = CSGeometry()
 
     core_thickness = 0.1