Merge branch 'master' into occ_spline_tools

# Conflicts:
#	libsrc/occ/python_occ_shapes.cpp

CMakeLists.txt

@@ -254,6 +254,13 @@ if (USE_GUI)
     add_definitions(-DTCL -DOPENGL -DUSE_TOGL_2 -DUSE_TCL_STUBS -DUSE_TK_STUBS)
     include_directories(${TCL_INCLUDE_PATH})
     include_directories(${TK_INCLUDE_PATH})
+    if(NOT EXISTS ${TK_INCLUDE_PATH}/tkWin.h AND EXISTS ${TK_INCLUDE_PATH}/../win/tkWin.h)
+        include_directories(${TK_INCLUDE_PATH}/../win)
+    endif()
+    if(NOT EXISTS ${TK_INCLUDE_PATH}/x11/Xlib.h AND EXISTS ${TK_INCLUDE_PATH}/../xlib/X11/Xlib.h)
+        include_directories(${TK_INCLUDE_PATH}/../xlib)
+    endif()
+
     set(LIBTOGL togl)
 
     if(WIN32)

cmake/SuperBuild.cmake

@@ -50,7 +50,6 @@ set_vars(SUBPROJECT_CMAKE_ARGS CMAKE_C_COMPILER)
 set_vars(SUBPROJECT_CMAKE_ARGS CMAKE_CXX_COMPILER)
 set_vars(SUBPROJECT_CMAKE_ARGS CMAKE_BUILD_TYPE)
 
-set(SUBPROJECT_CMAKE_ARGS "${SUBPROJECT_CMAKE_ARGS};-DCMAKE_OSX_ARCHITECTURES=${CMAKE_OSX_ARCHITECTURES}" CACHE INTERNAL "")
 set(SUBPROJECT_CMAKE_ARGS "${SUBPROJECT_CMAKE_ARGS};-DCMAKE_POSITION_INDEPENDENT_CODE=ON" CACHE INTERNAL "")
 
 if(USE_CCACHE)
@@ -101,8 +100,8 @@ if(BUILD_OCC)
 
   ExternalProject_Add(project_occ
     DEPENDS project_freetype
-    URL https://github.com/Open-Cascade-SAS/OCCT/archive/refs/tags/V7_5_0.zip
+    URL https://github.com/Open-Cascade-SAS/OCCT/archive/refs/tags/V7_6_0.zip
-    URL_MD5 a24e6d3cf2d24bf9347d2d4aee9dd80a
+    URL_MD5 37519251c99cb3469ccfa82a9241d528
     DOWNLOAD_DIR ${CMAKE_CURRENT_SOURCE_DIR}/external_dependencies
     ${SUBPROJECT_ARGS}
     CMAKE_ARGS
@@ -116,7 +115,7 @@ if(BUILD_OCC)
          -DBUILD_MODULE_DataExchange:BOOL=ON
          -DBUILD_MODULE_ApplicationFramework:BOOL=OFF
          -DBUILD_MODULE_Draw:BOOL=OFF
-         -DUSE_FREETYPE=OFF
+         -DUSE_FREETYPE=ON
          ${SUBPROJECT_CMAKE_ARGS}
     UPDATE_COMMAND ""
     )

cmake/external_projects/tcltk.cmake

@@ -35,7 +35,7 @@ ExternalProject_Add(project_tk
 )
 
 set(TCL_INCLUDE_PATH ${TCL_DIR}/generic)
-set(TK_INCLUDE_PATH ${TK_DIR}/generic;${TK_DIR}/xlib;${TK_DIR}/win)
+set(TK_INCLUDE_PATH ${TK_DIR}/generic)
 list(APPEND NETGEN_DEPENDENCIES project_tcl project_tk)
 
 if(APPLE OR WIN32)

external_dependencies/pybind11

@@ -1 +1 @@
-Subproject commit c16da993094988141101ac4d96a9b2c92f9ac714
+Subproject commit e7e2c79f3f520f78ffc39fcb34f7919003102733

libsrc/core/python_ngcore_export.cpp

@@ -37,12 +37,14 @@ PYBIND11_MODULE(pyngcore, m) // NOLINT
     .def("__len__", &BitArray::Size)
     .def("__getitem__", [] (BitArray & self, int i)
                                          {
+                                           if (i < 0) i+=self.Size();
                                            if (i < 0 || i >= self.Size())
                                              throw py::index_error();
                                            return self.Test(i);
                                          }, py::arg("pos"), "Returns bit from given position")
     .def("__setitem__", [] (BitArray & self, int i, bool b)
                                          {
+                                           if (i < 0) i+=self.Size();
                                            if (i < 0 || i >= self.Size())
                                              throw py::index_error();
                                            if (b) self.SetBit(i); else self.Clear(i);

libsrc/geom2d/csg2d.cpp

@@ -2185,6 +2185,7 @@ shared_ptr<netgen::SplineGeometry2d> CSG2d :: GenerateSplineGeometry()
     if(!is_solid_degenerated)
     {
       geo->SetMaterial(dom, s.name);
+      geo->SetDomainMaxh(dom, s.maxh);
       if(s.layer != 1)
           geo->SetDomainLayer(dom, s.layer);
     }

libsrc/geom2d/csg2d.hpp

@@ -632,6 +632,7 @@ struct Solid2d
 
   int layer = 1;
   string name = MAT_DEFAULT;
+  double maxh = MAXH_DEFAULT;
 
   Solid2d() = default;
   Solid2d(string name_) : name(name_) {}
@@ -697,6 +698,7 @@ struct Solid2d
 
   Solid2d & Maxh(double maxh)
   {
+    this->maxh = maxh;
     for(auto & p : polys)
       for(auto v : p.Vertices(ALL))
         v->info.maxh = maxh;

libsrc/geom2d/genmesh2d.cpp

@@ -429,7 +429,7 @@ namespace netgen
     PrintMessage (1, "Generate Mesh from spline geometry");
 
     Box<2> bbox = geometry.GetBoundingBox ();
-
+    bbox.Increase (1e-2*bbox.Diam());
     t_h.Start();
     if (bbox.Diam() < mp.maxh) 
       mp.maxh = bbox.Diam();

libsrc/gprim/adtree.hpp

@@ -817,6 +817,8 @@ public:
       : BoxTree(box.PMin(), box.PMax())
     { }
 
+  double GetTolerance() { return tol; }
+
   size_t GetNLeaves()
     {
       return n_leaves;

libsrc/meshing/basegeom.cpp

@@ -1,9 +1,33 @@
+#include <set>
+
 #include <mystdlib.h>
 #include "meshing.hpp"
 #include <core/register_archive.hpp>
 
 namespace netgen
 {
+  struct PointTree
+  {
+      BoxTree<3> tree;
+
+      PointTree( Box<3> bb ) : tree(bb) {}
+
+      void Insert(Point<3> p, PointIndex n)
+      {
+          tree.Insert(p, p, n);
+      }
+
+      PointIndex Find(Point<3> p) const
+      {
+          ArrayMem<int, 1> points;
+          tree.GetIntersecting(p, p, points);
+          if(points.Size()==0)
+              throw Exception("cannot find mapped point");
+          return points[0];
+      }
+
+      double GetTolerance() { return tree.GetTolerance(); }
+  };
 
   DLL_HEADER GeometryRegisterArray geometryregister;
   //DLL_HEADER NgArray<GeometryRegister*> geometryregister;
@@ -11,6 +35,78 @@ namespace netgen
   GeometryRegister :: ~GeometryRegister()
   { ; }
 
+  bool GeometryShape :: IsMappedShape( const GeometryShape & other_, const Transformation<3> & trafo, double tol ) const
+  {
+      throw Exception("GeometryShape::IsMappedShape not implemented for class " + Demangle(typeid(this).name()));
+  }
+
+  bool GeometryVertex :: IsMappedShape( const GeometryShape & other_, const Transformation<3> & trafo, double tol ) const
+  {
+      const auto other_ptr = dynamic_cast<const GeometryVertex*>(&other_);
+      if(!other_ptr)
+          return false;
+
+      return Dist(trafo(GetPoint()), other_ptr->GetPoint()) < tol;
+  }
+
+  bool GeometryEdge :: IsMappedShape( const GeometryShape & other_, const Transformation<3> & trafo, double tol ) const
+  {
+      const auto other_ptr = dynamic_cast<const GeometryEdge*>(&other_);
+      if(!other_ptr)
+          return false;
+      auto & e = *other_ptr;
+
+      if(tol < Dist(trafo(GetCenter()), e.GetCenter()))
+          return false;
+
+      auto v0 = trafo(GetStartVertex().GetPoint());
+      auto v1 = trafo(GetEndVertex().GetPoint());
+      auto w0 = e.GetStartVertex().GetPoint();
+      auto w1 = e.GetEndVertex().GetPoint();
+
+      // have two closed edges, use midpoints to compare
+      if(Dist(v0,v1) < tol && Dist(w0,w1) < tol)
+      {
+          v1 = trafo(GetPoint(0.5));
+          w1 = other_ptr->GetPoint(0.5);
+      }
+
+      return( (Dist(v0, w0) < tol && Dist(v1, w1) < tol) ||
+              (Dist(v0, w1) < tol && Dist(v1, w0) < tol) );
+  }
+
+  bool GeometryFace :: IsMappedShape( const GeometryShape & other_, const Transformation<3> & trafo, double tol ) const
+  {
+      const auto other_ptr = dynamic_cast<const GeometryFace*>(&other_);
+      if(!other_ptr)
+          return false;
+      auto & f = *other_ptr;
+
+      if(tol < Dist(GetCenter(), f.GetCenter()))
+          return false;
+
+      // simple check: check if there is a bijective mapping of mapped edges
+      auto & other_edges = f.edges;
+      if(edges.Size() != other_edges.Size())
+          return false;
+
+      auto nedges = edges.Size();
+      Array<bool> is_mapped(nedges);
+      is_mapped = false;
+
+      for(auto e : edges)
+      {
+          int found_mapping = 0;
+          for(auto other_e : other_edges)
+              if(e->IsMappedShape(*other_e, trafo, tol))
+                  found_mapping++;
+          if(found_mapping != 1)
+              return false;
+      }
+
+      return true;
+  }
+
   void GeometryFace :: RestrictHTrig(Mesh& mesh,
                                      const PointGeomInfo& gi0,
                                      const PointGeomInfo& gi1,
@@ -103,6 +199,123 @@ namespace netgen
     }
   };
 
+  void NetgenGeometry :: Clear()
+  {
+      vertex_map.clear();
+      edge_map.clear();
+      face_map.clear();
+      solid_map.clear();
+
+      vertices.SetSize0();
+      edges.SetSize0();
+      faces.SetSize0();
+      solids.SetSize0();
+  }
+
+  void NetgenGeometry :: ProcessIdentifications()
+  {
+      for(auto i : Range(vertices))
+          vertices[i]->nr = i;
+      for(auto i : Range(edges))
+          edges[i]->nr = i;
+      for(auto i : Range(faces))
+          faces[i]->nr = i;
+      for(auto i : Range(solids))
+          solids[i]->nr = i;
+
+      auto mirror_identifications = [&] ( auto & shapes )
+      {
+          for(auto i : Range(shapes))
+          {
+              auto &s = shapes[i];
+              s->nr = i;
+              for(auto & ident : s->identifications)
+                  if(s.get() == ident.from)
+                      ident.to->identifications.Append(ident);
+          }
+      };
+
+      auto tol = 1e-8 * bounding_box.Diam();
+      for(auto & f : faces)
+        for(auto & ident: f->identifications)
+          for(auto e : static_cast<GeometryFace*>(ident.from)->edges)
+            for(auto e_other : static_cast<GeometryFace*>(ident.to)->edges)
+              if(e->IsMappedShape(*e_other, ident.trafo, tol))
+                e->identifications.Append( {e, e_other, ident.trafo, ident.type, ident.name} );
+
+      for(auto & e : edges)
+        for(auto & ident: e->identifications)
+          {
+              auto & from = static_cast<GeometryEdge&>(*ident.from);
+              auto & to = static_cast<GeometryEdge&>(*ident.to);
+
+              GeometryVertex * pfrom[] = { &from.GetStartVertex(), &from.GetEndVertex() };
+              GeometryVertex * pto[] = { &to.GetStartVertex(), &to.GetEndVertex() };
+
+              // swap points of other edge if necessary
+              Point<3> p_from0 = ident.trafo(from.GetStartVertex().GetPoint());
+              Point<3> p_from1 = ident.trafo(from.GetEndVertex().GetPoint());
+              Point<3> p_to0 = to.GetStartVertex().GetPoint();
+
+              if(Dist(p_from1, p_to0) < Dist(p_from0, p_to0))
+                  swap(pto[0], pto[1]);
+
+              for(auto i : Range(2))
+                  pfrom[i]->identifications.Append( {pfrom[i], pto[i], ident.trafo, ident.type, ident.name} );
+          }
+
+      mirror_identifications(vertices);
+      mirror_identifications(edges);
+      mirror_identifications(faces);
+
+
+      auto find_primary = [&] (auto & shapes)
+      {
+          for(auto &s : shapes)
+          {
+              s->primary = s.get();
+              s->primary_to_me = Transformation<3>{ Vec<3> {0,0,0} }; // init with identity
+          }
+
+          bool changed = true;
+
+          while(changed) {
+            changed = false;
+            for(auto &s : shapes)
+            {
+              auto current = s->primary;
+              for(auto & ident : current->identifications)
+              {
+                  bool need_inverse = ident.from == s.get();
+                  auto other = need_inverse ? ident.to : ident.from;
+                  if(other->nr < s->primary->nr)
+                  {
+                      auto trafo = ident.trafo;
+                      if(need_inverse)
+                          trafo = trafo.CalcInverse();
+                      s->primary = other;
+                      s->primary_to_me.Combine(trafo, s->primary_to_me);
+                      changed = true;
+                  }
+                  if(other->primary->nr < s->primary->nr)
+                  {
+                      auto trafo = ident.trafo;
+                      if(need_inverse)
+                          trafo = trafo.CalcInverse();
+                      s->primary = other->primary;
+                      s->primary_to_me.Combine(trafo, other->primary_to_me);
+                      changed = true;
+                  }
+              }
+            }
+          }
+      };
+
+      find_primary(vertices);
+      find_primary(edges);
+      find_primary(faces);
+  }
+
   void NetgenGeometry :: Analyse(Mesh& mesh,
                                  const MeshingParameters& mparam) const
   {
@@ -241,60 +454,13 @@ namespace netgen
     mesh.LoadLocalMeshSize(mparam.meshsizefilename);
   }
 
-  void NetgenGeometry :: FindEdges(Mesh& mesh,
+  void DivideEdge(GeometryEdge * edge, const MeshingParameters & mparam, const Mesh & mesh, Array<Point<3>> & points, Array<double> & params)
-                                   const MeshingParameters& mparam) const
   {
-    static Timer t1("MeshEdges"); RegionTimer regt(t1);
-    static Timer tdivide("Divide Edges");
       static Timer tdivedgesections("Divide edge sections");
-    const char* savetask = multithread.task;
+      static Timer tdivide("Divide Edges");
-    multithread.task = "Mesh Edges";
+      RegionTimer rt(tdivide);
-
-    // create face descriptors and set bc names
-    mesh.SetNBCNames(faces.Size());
-    for(auto i : Range(faces.Size()))
-      {
-        mesh.SetBCName(i, faces[i]->GetName());
-        // todo find attached solids
-        FaceDescriptor fd(i+1, 1, 0, i+1);
-        fd.SetBCName(mesh.GetBCNamePtr(i));
-        mesh.AddFaceDescriptor(fd);
-      }
-
-    std::map<size_t, PointIndex> vert2meshpt;
-    for(auto i : Range(vertices))
-      {
-        const auto& vert = *vertices[i];
-        MeshPoint mp(vert.GetPoint());
-        vert2meshpt[vert.GetHash()] = mesh.AddPoint(mp);
-      }
-
-    size_t segnr = 0;
-    for(auto facenr : Range(faces.Size()))
-      {
-        const auto& face = *faces[facenr];
-        for(auto facebndnr : Range(face.GetNBoundaries()))
-          {
-            auto boundary = face.GetBoundary(facebndnr);
-            for(auto enr : Range(boundary))
-              {
-                multithread.percent = 100. * ((double(enr)/boundary.Size() + facebndnr)/face.GetNBoundaries() + facenr)/faces.Size();
-                const auto& oriented_edge = *boundary[enr];
-                auto edgenr = GetEdgeIndex(oriented_edge);
-                const auto& edge = edges[edgenr];
-                PointIndex startp, endp;
-                // throws if points are not found
-                startp = vert2meshpt.at(edge->GetStartVertex().GetHash());
-                endp = vert2meshpt.at(edge->GetEndVertex().GetHash());
-
-                // ignore collapsed edges
-                if(startp == endp && edge->GetLength() < 1e-10 * bounding_box.Diam())
-                  continue;
-                Array<MeshPoint> mps;
-                Array<double> params;
       // -------------------- DivideEdge -----------------
-                static constexpr size_t divide_edge_sections = 1000;
+      static constexpr size_t divide_edge_sections = 10000;
-                tdivide.Start();
       double hvalue[divide_edge_sections+1];
       hvalue[0] = 0;
 
@@ -309,7 +475,7 @@ namespace netgen
       }
       int nsubedges = max2(1, int(floor(hvalue[divide_edge_sections]+0.5)));
       tdivedgesections.Stop();
-                mps.SetSize(nsubedges-1);
+      points.SetSize(nsubedges-1);
       params.SetSize(nsubedges+1);
 
       int i = 1;
@@ -319,14 +485,14 @@ namespace netgen
           if (hvalue[i1]/hvalue[divide_edge_sections]*nsubedges >= i)
           {
               params[i] = (double(i1)/divide_edge_sections);
-                        mps[i-1] = MeshPoint(edge->GetPoint(params[i]));
+              points[i-1] = MeshPoint(edge->GetPoint(params[i]));
               i++;
           }
           i1++;
           if (i1 > divide_edge_sections)
           {
               nsubedges = i;
-                        mps.SetSize(nsubedges-1);
+              points.SetSize(nsubedges-1);
               params.SetSize(nsubedges+1);
               cout << "divide edge: local h too small" << endl;
           }
@@ -339,32 +505,141 @@ namespace netgen
       if(params[nsubedges] <= params[nsubedges-1])
       {
           cout << "CORRECTED" << endl;
-                    mps.SetSize (nsubedges-2);
+          points.SetSize (nsubedges-2);
           params.SetSize (nsubedges);
           params[nsubedges-1] = 1.;
       }
-                tdivide.Stop();
+  }
+
+  void NetgenGeometry :: FindEdges(Mesh& mesh,
+                                   const MeshingParameters& mparam) const
+  {
+    static Timer t1("MeshEdges"); RegionTimer regt(t1);
+    const char* savetask = multithread.task;
+    multithread.task = "Mesh Edges";
+
+    PointTree tree( bounding_box );
+
+    auto & identifications = mesh.GetIdentifications();
+
+    std::map<size_t, PointIndex> vert2meshpt;
+    for(auto & vert : vertices)
+      {
+        auto pi = mesh.AddPoint(vert->GetPoint());
+        tree.Insert(mesh[pi], pi);
+        vert2meshpt[vert->GetHash()] = pi;
+        mesh[pi].Singularity(vert->properties.hpref);
+
+        if(vert->properties.name)
+        {
+          Element0d el(pi, pi);
+          el.name = vert->properties.GetName();
+          mesh.SetCD3Name(pi, el.name);
+          mesh.pointelements.Append (el);
+        }
+      }
+
+    for(auto & vert : vertices)
+        for(auto & ident : vert->identifications)
+            identifications.Add(vert2meshpt[ident.from->GetHash()],
+                                vert2meshpt[ident.to->GetHash()],
+                                ident.name,
+                                ident.type);
+
+    size_t segnr = 0;
+    auto nedges = edges.Size();
+    Array<Array<PointIndex>> all_pnums(nedges);
+    Array<Array<double>> all_params(nedges);
+
+    for (auto edgenr : Range(edges))
+    {
+        auto edge = edges[edgenr].get();
+        PointIndex startp, endp;
+        // throws if points are not found
+        startp = vert2meshpt.at(edge->GetStartVertex().GetHash());
+        endp = vert2meshpt.at(edge->GetEndVertex().GetHash());
+
+        // ignore collapsed edges
+        if(startp == endp && edge->GetLength() < 1e-10 * bounding_box.Diam())
+            continue;
+
         // ----------- Add Points to mesh and create segments -----
-                Array<PointIndex> pnums(mps.Size() + 2);
+        auto & pnums = all_pnums[edgenr];
-                pnums[0] = startp;
+        auto & params = all_params[edgenr];
-                pnums[mps.Size()+1] = endp;
+        Array<Point<3>> edge_points;
+        Array<double> edge_params;
 
-                double eps = bounding_box.Diam() * 1e-8;
+        if(edge->primary == edge)
-
-                for(auto i : Range(mps))
         {
-                    bool exists = false;
+            // check if start and end vertex are identified (if so, we only insert one segement and do z-refinement later)
-                    for(auto pi : Range(mesh.Points()))
+            bool is_identified_edge = false;
+            auto v0 = vertices[edge->GetStartVertex().nr].get();
+            auto v1 = vertices[edge->GetEndVertex().nr].get();
+            for(auto & ident : v0->identifications)
             {
-                        if((mesh[pi] - mps[i]).Length() < eps)
+                auto other = ident.from == v0 ? ident.to : ident.from;
+                if(other->nr == v1->nr && ident.type == Identifications::CLOSESURFACES)
                 {
-                            exists = true;
+                    is_identified_edge = true;
-                            pnums[i+1] = pi;
                     break;
                 }
             }
-                    if(!exists)
+
-                      pnums[i+1] = mesh.AddPoint(mps[i]);
+            if(is_identified_edge)
+            {
+                params.SetSize(2);
+                params[0] = 0.;
+                params[1] = 1.;
+            }
+            else
+            {
+                DivideEdge(edge, mparam, mesh, edge_points, params);
+            }
+        }
+        else
+        {
+            auto nr_primary = edge->primary->nr;
+            auto & pnums_primary = all_pnums[nr_primary];
+            auto & params_primary = all_params[nr_primary];
+            auto trafo = edge->primary_to_me;
+
+            auto np = pnums_primary.Size();
+            edge_points.SetSize(np-2);
+            edge_params.SetSize(np-2);
+            for(auto i : Range(np-2))
+            {
+                edge_points[i] = trafo(mesh[pnums_primary[i+1]]);
+                EdgePointGeomInfo gi;
+                edge->ProjectPoint(edge_points[i], &gi);
+                edge_params[i] = gi.dist;
+            }
+
+            // reverse entries if we have decreasing parameters
+            if(edge_params.Size()>2 && edge_params[0] > edge_params.Last())
+                for(auto i : Range((np-2)/2))
+                {
+                    swap(edge_points[i], edge_points[np-3-i]);
+                    swap(edge_params[i], edge_params[np-3-i]);
+                }
+
+            params.SetSize(edge_params.Size()+2);
+            params[0] = 0.;
+            params.Last() = 1.;
+
+            for(auto i : Range(edge_params))
+                params[i+1] = edge_params[i];
+        }
+
+        pnums.SetSize(edge_points.Size() + 2);
+        pnums[0] = startp;
+        pnums.Last() = endp;
+
+
+        for(auto i : Range(edge_points))
+        {
+            auto pi = mesh.AddPoint(edge_points[i]);
+            tree.Insert(mesh[pi], pi);
+            pnums[i+1] = pi;
         }
 
         for(auto i : Range(pnums.Size()-1))
@@ -373,28 +648,33 @@ namespace netgen
             Segment seg;
             seg[0] = pnums[i];
             seg[1] = pnums[i+1];
-                    seg.edgenr = segnr;
+            seg.edgenr = edgenr+1;
+            seg.si = edgenr+1;
             seg.epgeominfo[0].dist = params[i];
             seg.epgeominfo[1].dist = params[i+1];
             seg.epgeominfo[0].edgenr = edgenr;
             seg.epgeominfo[1].edgenr = edgenr;
-                    seg.si = facenr+1;
+            seg.singedge_left = edge->properties.hpref;
-                    seg.surfnr1 = facenr+1;
+            seg.singedge_right = edge->properties.hpref;
-
+            seg.domin = edge->domin+1;
-                    // TODO: implement functionality to transfer edge parameter t to face parameters u,v
+            seg.domout = edge->domout+1;
-                    for(auto j : Range(2))
-                      face.CalcEdgePointGI(*edge, params[i+j],
-                                           seg.epgeominfo[j]);
-
-                    if(!oriented_edge.OrientedLikeGlobal())
-                      {
-                        swap (seg[0], seg[1]);
-                        swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
-                        swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
-                        swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
-                      }
             mesh.AddSegment(seg);
         }
+        mesh.SetCD2Name(edgenr+1, edge->properties.GetName());
+    }
+
+    for (auto & edge : edges)
+    {
+        // identify points on edge
+        for(auto & ident : edge->identifications)
+          if(ident.from == edge.get())
+          {
+            auto & pnums = all_pnums[edge->nr];
+            // start and end vertex are already identified
+            for(auto pi : pnums.Range(1, pnums.Size()-1))
+            {
+                auto pi_other = tree.Find(ident.trafo(mesh[pi]));
+                identifications.Add(pi, pi_other, ident.name, ident.type);
             }
           }
     }
@@ -402,15 +682,8 @@ namespace netgen
     multithread.task = savetask;
   }
 
-  void NetgenGeometry :: MeshSurface(Mesh& mesh,
+  bool NetgenGeometry :: MeshFace(Mesh& mesh, const MeshingParameters& mparam,
-                                     const MeshingParameters& mparam) const
+                     int k, FlatArray<int, PointIndex> glob2loc) const
-  {
-    static Timer t1("Surface Meshing"); RegionTimer regt(t1);
-    const char* savetask = multithread.task;
-    multithread.task = "Mesh Surface";
-
-    Array<int, PointIndex> glob2loc(mesh.GetNP());
-    for(auto k : Range(faces))
   {
     multithread.percent = 100. * k/faces.Size();
     const auto& face = *faces[k];
@@ -420,9 +693,8 @@ namespace netgen
     glob2loc = 0;
     int cntp = 0;
 
-        for(auto& seg : mesh.LineSegments())
+    auto segments = face.GetBoundary(mesh);
-          {
+    for(auto& seg : segments)
-            if(seg.si == k+1)
       {
         for(auto j : Range(2))
           {
@@ -435,10 +707,7 @@ namespace netgen
               }
           }
       }
-          }
+    for(auto & seg : segments)
-        for(auto & seg : mesh.LineSegments())
-          {
-            if(seg.si == k+1)
       {
         PointGeomInfo gi0, gi1;
         gi0.trignum = gi1.trignum = k+1;
@@ -450,7 +719,6 @@ namespace netgen
                                    glob2loc[seg[1]],
                                    gi0, gi1);
       }
-          }
 
     // TODO Set max area 2* area of face
 
@@ -464,10 +732,248 @@ namespace netgen
       {
         mesh.SurfaceElements()[i].SetIndex(k+1);
       }
+    return res != MESHING2_OK;
   }
+
+  void NetgenGeometry :: MeshSurface(Mesh& mesh,
+                                     const MeshingParameters& mparam) const
+  {
+    static Timer t1("Surface Meshing"); RegionTimer regt(t1);
+    const char* savetask = multithread.task;
+    multithread.task = "Mesh Surface";
+    mesh.ClearFaceDescriptors();
+
+    size_t n_failed_faces = 0;
+    Array<int, PointIndex> glob2loc(mesh.GetNP());
+    for(auto k : Range(faces))
+    {
+        auto & face = *faces[k];
+        FaceDescriptor fd(k+1, face.domin+1, face.domout+1, k+1);
+        mesh.AddFaceDescriptor(fd);
+        mesh.SetBCName(k, face.properties.GetName());
+        if(face.primary == &face)
+        {
+            // check if this face connects two identified closesurfaces
+            bool is_connecting_closesurfaces = false;
+            auto & idents = mesh.GetIdentifications();
+            std::set<int> relevant_edges;
+            auto segments = face.GetBoundary(mesh);
+            for(const auto &s : segments)
+                relevant_edges.insert(s.edgenr-1);
+
+            Array<bool, PointIndex> is_point_in_tree(mesh.Points().Size());
+            is_point_in_tree = false;
+            PointTree tree( bounding_box );
+            for(const auto &s : segments)
+                for(auto pi : s.PNums())
+                    if(!is_point_in_tree[pi])
+                    {
+                        tree.Insert(mesh[pi], pi);
+                        is_point_in_tree[pi] = true;
+                    }
+
+            Array<int> mapped_edges(edges.Size());
+            constexpr int UNINITIALIZED = -2;
+            constexpr int NOT_MAPPED = -1;
+            mapped_edges = UNINITIALIZED;
+
+            Transformation<3> trafo;
+
+            for(const auto &s : segments)
+            {
+                auto edgenr = s.edgenr-1;
+                auto & edge = *edges[edgenr];
+                ShapeIdentification *edge_mapping;
+
+                // have edgenr first time, search for closesurface identification
+                if(mapped_edges[edgenr] == UNINITIALIZED)
+                {
+                    mapped_edges[edgenr] = NOT_MAPPED;
+                    for(auto & edge_ident : edge.identifications)
+                    {
+                        if(edge_ident.type == Identifications::CLOSESURFACES && 
+                                edge_ident.from->nr == edgenr &&
+                                relevant_edges.count(edge_ident.to->nr) > 0
+                          )
+                        {
+                            trafo = edge_ident.trafo;
+                            mapped_edges[edgenr] = edge_ident.to->nr;
+                            is_connecting_closesurfaces = true;
+                            break;
+                        }
+                    }
+                }
+
+                // this edge has a closesurface mapping to another -> make connecting quad
+                if(mapped_edges[edgenr] != NOT_MAPPED)
+                {
+                    Element2d sel(4);
+                    sel[0] = s[0];
+                    sel[1] = s[1];
+                    sel[2] = tree.Find(trafo(mesh[s[1]]));
+                    sel[3] = tree.Find(trafo(mesh[s[0]]));
+                    for(auto i : Range(4))
+                        sel.GeomInfo()[i] = face.Project(mesh[sel[i]]);
+
+                    sel.SetIndex(face.nr+1);
+                    mesh.AddSurfaceElement(sel);
+                }
+            }
+
+            if(!is_connecting_closesurfaces)
+                if(MeshFace(mesh, mparam, k, glob2loc))
+                    n_failed_faces++;
+        }
+    }
+
+    if(n_failed_faces) 
+    {
+        cout << "WARNING! NOT ALL FACES HAVE BEEN MESHED" << endl;
+        cout << "SURFACE MESHING ERROR OCCURRED IN " << n_failed_faces << " FACES:" << endl;
+        return;
+    }
+
+    if (mparam.perfstepsend >= MESHCONST_OPTSURFACE)
+    {
+      mesh.CalcSurfacesOfNode();
+      OptimizeSurface(mesh, mparam);
+    }
+
+    bool have_identifications = false;
+    for(auto & face : faces)
+        if(face->primary != face.get())
+        {
+            have_identifications = true;
+            MapSurfaceMesh(mesh, *face);
+        }
+
+    // identify points on faces
+    if(have_identifications)
+    {
+        mesh.CalcSurfacesOfNode();
+        BitArray is_identified_face(faces.Size());
+        is_identified_face = false;
+        for(auto & face : faces)
+            for(auto & ident : face->identifications)
+            {
+                is_identified_face.SetBit(ident.from->nr);
+                is_identified_face.SetBit(ident.to->nr);
+            }
+
+        PointTree tree( bounding_box );
+        Array<int, PointIndex> pi_to_face(mesh.GetNP());
+        pi_to_face = -1;
+        Array<SurfaceElementIndex> si_of_face;
+        Array<Array<PointIndex>> pi_of_face(faces.Size());
+        for(auto & face : faces)
+            if(is_identified_face[face->nr])
+            {
+                mesh.GetSurfaceElementsOfFace(face->nr+1, si_of_face);
+                for(auto si : si_of_face)
+                    for(auto pi : mesh[si].PNums())
+                    {
+                        if(mesh[pi].Type() == SURFACEPOINT && pi_to_face[pi]==-1)
+                        {
+                            pi_to_face[pi] = face->nr;
+                            tree.Insert(mesh[pi], pi);
+                            pi_of_face[face->nr].Append(pi);
+                        }
+                    }
+            }
+
+        auto & mesh_ident = mesh.GetIdentifications();
+        for(auto & face : faces)
+            for(auto & ident : face->identifications)
+            {
+                if(ident.from == face.get())
+                    for(auto pi : pi_of_face[face->nr])
+                    {
+                        auto pi_other = tree.Find(ident.trafo(mesh[pi]));
+                        mesh_ident.Add(pi, pi_other, ident.name, ident.type);
+                    }
+            }
+    }
+
+    mesh.CalcSurfacesOfNode();
     multithread.task = savetask;
   }
 
+  void NetgenGeometry :: MapSurfaceMesh( Mesh & mesh, const GeometryFace & dst ) const
+  {
+    static Timer timer("MapSurfaceMesh");
+    RegionTimer rt(timer);
+
+    const auto & src = dynamic_cast<const GeometryFace&>(*dst.primary);
+    auto trafo = dst.primary_to_me;
+
+    PrintMessage(2, "Map face ", src.nr+1, " -> ", dst.nr+1);
+
+    // point map from src to dst
+    Array<PointIndex, PointIndex> pmap(mesh.Points().Size());
+    pmap = PointIndex::INVALID;
+
+    // first map points on edges (mapped points alread in mesh, use search tree)
+    Array<bool, PointIndex> is_point_in_tree(mesh.Points().Size());
+    is_point_in_tree = false;
+    PointTree tree( bounding_box );
+
+    for (Segment & seg : src.GetBoundary(mesh))
+        for(auto i : Range(2))
+          {
+            auto pi = seg[i];
+            if(!is_point_in_tree[pi])
+            {
+              tree.Insert(trafo(mesh[pi]), pi);
+              is_point_in_tree[pi] = true;
+            }
+          }
+
+    for (Segment & seg : dst.GetBoundary(mesh))
+        for(auto i : Range(2))
+          {
+            auto pi = seg[i];
+            if(pmap[pi].IsValid())
+              continue;
+
+            pmap[tree.Find(mesh[pi])] = pi;
+          }
+
+    xbool do_invert = maybe;
+
+    // now insert mapped surface elements
+    for(auto sei : mesh.SurfaceElements().Range())
+      {
+        auto sel = mesh[sei];
+        if(sel.GetIndex() != src.nr+1)
+          continue;
+
+        if(do_invert.IsMaybe())
+        {
+            auto n_src = src.GetNormal(mesh[sel[0]]);
+            auto n_dist = dst.GetNormal(mesh[sel[0]]);
+            Mat<3> normal_matrix;
+            CalcInverse(Trans(trafo.GetMatrix()), normal_matrix);
+            do_invert = n_src * (normal_matrix * n_dist) < 0.0;
+        }
+        auto sel_new = sel;
+        sel_new.SetIndex(dst.nr+1);
+        for(auto i : Range(sel.PNums()))
+          {
+            auto pi = sel[i];
+            if(!pmap[pi].IsValid())
+              {
+                pmap[pi] = mesh.AddPoint(trafo(mesh[pi]), 1, SURFACEPOINT);
+              }
+              sel_new[i] = pmap[pi];
+          }
+          if(do_invert.IsTrue())
+              sel_new.Invert();
+          for(auto i : Range(sel.PNums()))
+              dst.CalcPointGeomInfo(mesh[sel_new[i]], sel_new.GeomInfo()[i]);
+          mesh.AddSurfaceElement(sel_new);
+      }
+  }
+
   void NetgenGeometry :: OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const
   {
     const auto savetask = multithread.task;
@@ -505,6 +1011,13 @@ namespace netgen
     multithread.task = savetask;
   }
 
+  void NetgenGeometry :: FinalizeMesh(Mesh& mesh) const
+  {
+    for (int i = 0; i < mesh.GetNDomains(); i++)
+      if (auto name = solids[i]->properties.name)
+        mesh.SetMaterial (i+1, *name);
+  }
+  
   shared_ptr<NetgenGeometry> GeometryRegisterArray :: LoadFromMeshFile (istream & ist) const
   {
     if (!ist.good())
@@ -567,18 +1080,17 @@ namespace netgen
     if (mparam.perfstepsstart <= MESHCONST_MESHSURFACE)
       {
         MeshSurface(*mesh, mparam);
-        mesh->CalcSurfacesOfNode();
       }
 
-    if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHSURFACE)
-      return 0;
-
-    if (mparam.perfstepsstart <= MESHCONST_OPTSURFACE)
-      OptimizeSurface(*mesh, mparam);
-
     if (multithread.terminate || mparam.perfstepsend <= MESHCONST_OPTSURFACE)
       return 0;
 
+    if(dimension == 2)
+    {
+        FinalizeMesh(*mesh);
+        mesh->SetDimension(2);
+        return 0;
+    }
 
     if(mparam.perfstepsstart <= MESHCONST_MESHVOLUME)
       {

libsrc/meshing/basegeom.hpp

@@ -12,26 +12,82 @@ struct Tcl_Interp;
 
 namespace netgen
 {
-  class GeometryVertex
+  struct ShapeProperties
   {
-  public:
+    optional<string> name;
-    virtual ~GeometryVertex() {}
+    optional<Vec<4>> col;
-    virtual Point<3> GetPoint() const = 0;
+    double maxh = 1e99;
-    virtual size_t GetHash() const = 0;
+    double hpref = 0;  // number of hp refinement levels (will be multiplied by factor later)
+    void Merge(const ShapeProperties & prop2)
+    {
+      if (prop2.name) name = prop2.name;
+      if (prop2.col) col = prop2.col;
+      maxh = min2(maxh, prop2.maxh);
+      hpref = max2(hpref, prop2.hpref);
+    }
+
+    string GetName() const { return name ? *name : "default"; }
+    Vec<4> GetColor() { return col ? *col : Vec<4>{0., 1., 0., 1.}; }
+
+    void DoArchive(Archive& ar)
+    {
+        ar & name & col & maxh & hpref;
+    }
   };
 
-  class GeometryEdge
+  class GeometryShape;
+
+  struct ShapeIdentification
+  {
+    GeometryShape * from;
+    GeometryShape * to;
+    Transformation<3> trafo;
+    Identifications::ID_TYPE type;
+    string name = "";
+  };
+
+  class DLL_HEADER GeometryShape
   {
   public:
-    virtual ~GeometryEdge() {}
+    int nr = -1;
-    virtual const GeometryVertex& GetStartVertex() const = 0;
+    ShapeProperties properties;
-    virtual const GeometryVertex& GetEndVertex() const = 0;
+    Array<ShapeIdentification> identifications;
+    GeometryShape * primary;
+    Transformation<3> primary_to_me;
+
+    virtual ~GeometryShape() {}
+    virtual size_t GetHash() const = 0;
+    virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const;
+  };
+
+
+  class DLL_HEADER GeometryVertex : public GeometryShape
+  {
+  public:
+    virtual Point<3> GetPoint() const = 0;
+    virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const override;
+  };
+
+  class DLL_HEADER GeometryEdge : public GeometryShape
+  {
+  protected:
+      GeometryVertex *start, *end;
+  public:
+    int domin=-1, domout=-1;
+
+    GeometryEdge( GeometryVertex &start_, GeometryVertex &end_ )
+        : start(&start_), end(&end_)
+    {}
+
+    virtual const GeometryVertex& GetStartVertex() const { return *start; }
+    virtual const GeometryVertex& GetEndVertex() const { return *end; }
+    virtual GeometryVertex& GetStartVertex() { return *start; }
+    virtual GeometryVertex& GetEndVertex() { return *end; }
     virtual double GetLength() const = 0;
+    virtual Point<3> GetCenter() const = 0;
     virtual Point<3> GetPoint(double t) const = 0;
     // Calculate parameter step respecting edges sag value
     virtual double CalcStep(double t, double sag) const = 0;
-    virtual bool OrientedLikeGlobal() const = 0;
-    virtual size_t GetHash() const = 0;
     virtual void ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const = 0;
     virtual void PointBetween(const Point<3>& p1,
                               const Point<3>& p2,
@@ -46,15 +102,19 @@ namespace netgen
       ProjectPoint(newp, &newgi);
     }
     virtual Vec<3> GetTangent(double t) const = 0;
+    virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const override;
   };
 
-  class GeometryFace
+  class DLL_HEADER GeometryFace : public GeometryShape
   {
   public:
-    virtual ~GeometryFace() {}
+    Array<GeometryEdge*> edges;
+    int domin=-1, domout=-1;
+
+    virtual Point<3> GetCenter() const = 0;
     virtual size_t GetNBoundaries() const = 0;
-    virtual Array<unique_ptr<GeometryEdge>> GetBoundary(size_t index) const = 0;
+    virtual Array<Segment> GetBoundary(const Mesh& mesh) const = 0;
-    virtual string GetName() const { return "default"; }
+
     virtual PointGeomInfo Project(Point<3>& p) const = 0;
     // Project point using geo info. Fast if point is close to
     // parametrization in geo info.
@@ -93,6 +153,8 @@ namespace netgen
         newgi = Project(newp);
     }
 
+    virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const override;
+
   protected:
     void RestrictHTrig(Mesh& mesh,
                        const PointGeomInfo& gi0,
@@ -102,6 +164,9 @@ namespace netgen
                        int depth = 0, double h = 0.) const;
   };
 
+  class DLL_HEADER GeometrySolid : public GeometryShape
+  { };
+
   class DLL_HEADER NetgenGeometry
   {
     unique_ptr<Refinement> ref;
@@ -109,15 +174,31 @@ namespace netgen
     Array<unique_ptr<GeometryVertex>> vertices;
     Array<unique_ptr<GeometryEdge>> edges;
     Array<unique_ptr<GeometryFace>> faces;
+    Array<unique_ptr<GeometrySolid>> solids;
     Array<std::pair<Point<3>, double>> restricted_h;
     Box<3> bounding_box;
+    int dimension = 3;
+
+    std::map<size_t, GeometryVertex*> vertex_map;
+    std::map<size_t, GeometryEdge*> edge_map;
+    std::map<size_t, GeometryFace*> face_map;
+    std::map<size_t, GeometrySolid*> solid_map;
   public:
+
     NetgenGeometry()
     {
       ref = make_unique<Refinement>(*this);
     }
     virtual ~NetgenGeometry () { ; }
 
+    size_t GetNVertices() const { return vertices.Size(); }
+    size_t GetNEdges() const { return edges.Size(); }
+    size_t GetNFaces() const { return faces.Size(); }
+
+    const GeometryFace & GetFace(int i) const { return *faces[i]; }
+
+    void Clear();
+
     virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam);
 
     void RestrictH(const Point<3>& pnt, double maxh)
@@ -134,13 +215,17 @@ namespace netgen
   { throw NgException("DoArchive not implemented for " + Demangle(typeid(*this).name())); }
 
     virtual Mesh::GEOM_TYPE GetGeomType() const { return Mesh::NO_GEOM; }
+    virtual void ProcessIdentifications();
     virtual void Analyse(Mesh& mesh,
                          const MeshingParameters& mparam) const;
     virtual void FindEdges(Mesh& mesh, const MeshingParameters& mparam) const;
     virtual void MeshSurface(Mesh& mesh, const MeshingParameters& mparam) const;
+    virtual bool MeshFace(Mesh& mesh, const MeshingParameters& mparam,
+                     int nr, FlatArray<int, PointIndex> glob2loc) const;
+    virtual void MapSurfaceMesh( Mesh & mesh, const GeometryFace & dst ) const;
     virtual void OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const;
 
-    virtual void FinalizeMesh(Mesh& mesh) const {}
+    virtual void FinalizeMesh(Mesh& mesh) const;
 
     virtual PointGeomInfo ProjectPoint (int surfind, Point<3> & p) const
     {

libsrc/meshing/meshclass.cpp

@@ -665,6 +665,13 @@ namespace netgen
 
     outfile << "geomtype\n" << int(geomtype) << "\n";
 
+    outfile << "\n";
+    outfile << "# surfnr\tdomin\tdomout\ttlosurf\tbcprop\n";
+    outfile << "facedescriptors\n";
+    outfile << GetNFD() << "\n";
+    for(auto & fd : FaceDescriptors())
+        outfile << fd.SurfNr() << ' ' << fd.DomainIn() << ' ' << fd.DomainOut() << ' ' << fd.TLOSurface() << ' ' << fd.BCProperty() << '\n';
+
 
     outfile << "\n";
     outfile << "# surfnr    bcnr   domin  domout      np      p1      p2      p3"
@@ -1192,6 +1199,19 @@ namespace netgen
             geomtype = GEOM_TYPE(hi);
           }
 
+        if (strcmp (str, "facedescriptors") == 0)
+          {
+            int nfd;
+            infile >> nfd;
+            for(auto i : Range(nfd))
+            {
+                int surfnr, domin, domout, tlosurf, bcprop;
+                infile >> surfnr >> domin >> domout >> tlosurf >> bcprop;
+                auto faceind = AddFaceDescriptor (FaceDescriptor(surfnr, domin, domout, tlosurf));
+                GetFaceDescriptor(faceind).SetBCProperty(bcprop);
+            }
+          }
+
 
         if (strcmp (str, "surfaceelements") == 0 || strcmp (str, "surfaceelementsgi")==0 || strcmp (str, "surfaceelementsuv") == 0)
           {
@@ -6887,14 +6907,11 @@ namespace netgen
         int oldsize = bcnames.Size();
         bcnames.SetSize (bcnr+1);  // keeps contents
         for (int i = oldsize; i <= bcnr; i++)
-          bcnames[i] = nullptr;
+          bcnames[i] = new string("default");
       }
 
     if ( bcnames[bcnr] ) delete bcnames[bcnr];
-    if ( abcname != "default" )
     bcnames[bcnr] = new string ( abcname );
-    else
-      bcnames[bcnr] = nullptr;
 
     for (auto & fd : facedecoding)
       if (fd.BCProperty() <= bcnames.Size())

libsrc/meshing/meshfunc.cpp

@@ -156,8 +156,6 @@ namespace netgen
                   if(!pi0.IsValid() || !pi1.IsValid())
                       continue;
 
-                  if(pi1<pi0)
-                      Swap(pi0,pi1);
                   m_ident.Add(pi0, pi1, n);
               }
               m_ident.SetType( n, identifications.GetType(n) );
@@ -166,6 +164,59 @@ namespace netgen
       return ret;
   }
 
+  // Add between identified surface elements (only consider closesurface identifications)
+  void FillCloseSurface( MeshingData & md)
+  {
+      static Timer timer("FillCloseSurface"); RegionTimer rtimer(timer);
+
+      auto & mesh = md.mesh;
+      auto & identifications = mesh->GetIdentifications();
+      auto nmax = identifications.GetMaxNr();
+
+      bool have_closesurfaces = false;
+      for(auto i : Range(1,nmax+1))
+          if(identifications.GetType(i) == Identifications::CLOSESURFACES)
+              have_closesurfaces = true;
+      if(!have_closesurfaces)
+          return;
+
+      NgArray<int, PointIndex::BASE> map;
+      for(auto identnr : Range(1,nmax+1))
+      {
+          if(identifications.GetType(identnr) != Identifications::CLOSESURFACES)
+              continue;
+
+          identifications.GetMap(identnr, map);
+
+          for(auto & sel : mesh->SurfaceElements())
+          {
+              bool is_mapped = true;
+              for(auto pi : sel.PNums())
+                  if(!PointIndex(map[pi]).IsValid())
+                      is_mapped = false;
+
+              if(!is_mapped)
+                  continue;
+              
+              // in case we have symmetric mapping (used in csg), only map in one direction
+              if(map[map[sel[0]]] == sel[0] && map[sel[0]] < sel[0])
+                  continue;
+
+              // insert prism
+              auto np = sel.GetNP();
+              Element el(2*np);
+              for(auto i : Range(np))
+              {
+                  el[i] = sel[i];
+                  el[i+np] = map[sel[i]];
+              }
+
+              el.SetIndex(md.domain);
+              mesh->AddVolumeElement(el);
+          }
+      }
+  }
+
   void CloseOpenQuads( MeshingData & md)
   {
     auto & mesh = *md.mesh;
@@ -488,6 +539,9 @@ namespace netgen
            if (md[i].mesh->CheckOverlappingBoundary())
              throw NgException ("Stop meshing since boundary mesh is overlapping");
          
+         // TODO: FillCloseSurface is not working with CSG closesurfaces
+         if(md[i].mesh->GetGeometry()->GetGeomType() == Mesh::GEOM_OCC)
+           FillCloseSurface( md[i] );
          CloseOpenQuads( md[i] );
          MeshDomain(md[i]);
        });

libsrc/meshing/meshtype.cpp

@@ -2687,6 +2687,7 @@ namespace netgen
     identifiedpoints_nr.Set (tripl, 1);
 
     if (identnr > maxidentnr) maxidentnr = identnr;
+    names.SetSize(maxidentnr);
 
     if (identnr+1 > idpoints_table.Size())
       idpoints_table.ChangeSize (identnr+1);

libsrc/meshing/meshtype.hpp

@@ -1497,6 +1497,7 @@ namespace netgen
 
     /// number of identifications (or, actually used identifications ?)
     int maxidentnr;
+    Array<string> names;
 
   public:
     ///
@@ -1511,7 +1512,12 @@ namespace netgen
       identification nr identnr
     */
     DLL_HEADER void Add (PointIndex pi1, PointIndex pi2, int identnr);
-
+    void Add (PointIndex pi1, PointIndex pi2, string name, ID_TYPE type)
+    {
+        auto nr = GetNr(name);
+        Add(pi1, pi2, nr);
+        SetType(nr, type);
+    }
 
     int Get (PointIndex pi1, PointIndex pi2) const;
     int GetSymmetric (PointIndex pi1, PointIndex pi2) const;
@@ -1560,6 +1566,13 @@ namespace netgen
     ///
     int GetMaxNr () const { return maxidentnr; }  
 
+    int GetNr(string name)
+    {
+      if(!names.Contains(name))
+         names.Append(name);
+      return names.Pos(name)+1;
+    }
+
     /// remove secondorder
     void SetMaxPointNr (int maxpnum);
 

libsrc/meshing/python_mesh.cpp

@@ -117,6 +117,15 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
   m.def("_PushStatus", [](string s) { PushStatus(MyStr(s)); });
   m.def("_SetThreadPercentage", [](double percent) { SetThreadPercent(percent); });
 
+  py::enum_<Identifications::ID_TYPE>(m,"IdentificationType")
+    .value("UNDEFINED", Identifications::UNDEFINED)
+    .value("PERIODIC", Identifications::PERIODIC)
+    .value("CLOSESURFACES", Identifications::CLOSESURFACES)
+    .value("CLOSEEDGES", Identifications::CLOSEEDGES)
+    ;
+
+  py::implicitly_convertible<int, Identifications::ID_TYPE>();
+
   py::class_<NgMPI_Comm> (m, "MPI_Comm")
 #ifdef NG_MPI4PY
     .def(py::init([] (mpi4py_comm comm)
@@ -944,19 +953,19 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     .def ("GetCD3Name", &Mesh::GetCD3Name)
     .def ("SetCD3Name", &Mesh::SetCD3Name)
 
-    .def ("AddPointIdentification", [](Mesh & self, py::object pindex1, py::object pindex2, int identnr, int type)
+    .def ("AddPointIdentification", [](Mesh & self, py::object pindex1, py::object pindex2, int identnr, Identifications::ID_TYPE type)
                            {
 			     if(py::extract<PointIndex>(pindex1).check() && py::extract<PointIndex>(pindex2).check())
 			       {
 				 self.GetIdentifications().Add (py::extract<PointIndex>(pindex1)(), py::extract<PointIndex>(pindex2)(), identnr);
-				 self.GetIdentifications().SetType(identnr, Identifications::ID_TYPE(type)); // type = 2 ... periodic
+				 self.GetIdentifications().SetType(identnr, type); // type = 2 ... periodic
 			       }
                            },
           //py::default_call_policies(),
           py::arg("pid1"),
            py::arg("pid2"),
            py::arg("identnr"),
-           py::arg("type"))
+           py::arg("type")=Identifications::PERIODIC)
     .def("IdentifyPeriodicBoundaries", &Mesh::IdentifyPeriodicBoundaries,
          py::arg("face1"), py::arg("face2"), py::arg("mapping"), py::arg("point_tolerance") = -1.)
     .def("GetNrIdentifications", [](Mesh& self)

libsrc/occ/CMakeLists.txt

@@ -1,9 +1,12 @@
+if(USE_OCC)
+
 add_definitions(-DNGINTERFACE_EXPORTS)
 add_library(occ ${NG_LIB_TYPE}
         Partition_Inter2d.cxx Partition_Inter3d.cxx 
         Partition_Loop.cxx Partition_Loop2d.cxx Partition_Loop3d.cxx Partition_Spliter.cxx 
         occgenmesh.cpp occgeom.cpp occmeshsurf.cpp python_occ.cpp
         python_occ_basic.cpp python_occ_shapes.cpp
+        occ_face.cpp occ_edge.cpp occ_vertex.cpp occ_utils.cpp
         )
 if(USE_GUI)
         add_library(occvis ${NG_LIB_TYPE} vsocc.cpp)
@@ -14,11 +17,11 @@ target_link_libraries(occ PUBLIC ngcore PRIVATE "$<BUILD_INTERFACE:netgen_python
 
 if(NOT WIN32)
     target_link_libraries( occ PRIVATE ${OCC_LIBRARIES} )
-    if(USE_OCC AND APPLE)
+    if(APPLE)
         # Link AppKit in case OCE was built as static libraries
         find_library(AppKit AppKit)
         target_link_libraries( occ PRIVATE ${AppKit} )
-    endif(USE_OCC AND APPLE)
+    endif(APPLE)
     install( TARGETS occ ${NG_INSTALL_DIR})
     if (USE_GUI)
         target_link_libraries( occvis PUBLIC occ )
@@ -27,6 +30,9 @@ if(NOT WIN32)
 endif(NOT WIN32)
 
 install(FILES
-  occgeom.hpp occmeshsurf.hpp vsocc.hpp
+  occgeom.hpp occmeshsurf.hpp vsocc.hpp occ_utils.hpp
+  occ_vertex.hpp occ_edge.hpp occ_face.hpp occ_solid.hpp
   DESTINATION ${NG_INSTALL_DIR_INCLUDE}/occ COMPONENT netgen_devel
 )
+
+endif(USE_OCC)

libsrc/occ/occ_edge.cpp

@@ -0,0 +1,75 @@
+#include <BRepGProp.hxx>
+#include <BRep_Tool.hxx>
+#include <GeomAPI_ProjectPointOnCurve.hxx>
+
+#include "occ_edge.hpp"
+#include "occgeom.hpp"
+
+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);
+        BRepGProp::LinearProperties(edge, props);
+
+        auto verts = GetVertices(edge);
+        if(verts.size() != 2)
+            throw Exception("OCC edge does not have 2 vertices");
+
+        if(start != end)
+        {
+            // swap start/end if necessary
+            double d00 = Dist(GetPoint(0), start->GetPoint());
+            double d01 = Dist(GetPoint(0), end->GetPoint());
+            if(d01 < d00)
+                swap(start, end);
+        }
+    }
+
+    double OCCEdge::GetLength() const
+    {
+        return props.Mass();
+    }
+
+    Point<3> OCCEdge::GetCenter() const
+    {
+        return occ2ng( props.CentreOfMass() );
+    }
+
+    Point<3> OCCEdge::GetPoint(double t) const
+    {
+        return occ2ng( curve->Value(s0+t*(s1-s0)) );
+    }
+
+    double OCCEdge::CalcStep(double t, double sag) const
+    {
+        throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__));
+    }
+
+    size_t OCCEdge::GetHash() const
+    {
+        return reinterpret_cast<size_t>(tedge.get());
+    }
+
+    void OCCEdge::ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const
+    {
+        auto pnt = ng2occ(p);
+        GeomAPI_ProjectPointOnCurve proj(pnt, curve);
+        pnt = proj.NearestPoint();
+        if(gi)
+            gi->dist = (proj.LowerDistanceParameter() - s0)/(s1-s0);
+        p = occ2ng(pnt);
+    }
+
+    Vec<3> OCCEdge::GetTangent(double t) const
+    {
+        t = s0 + t*(s1-s0);
+        gp_Pnt p;
+        gp_Vec v;
+        curve->D1(t, p, v);
+        return occ2ng(v);
+    }
+}

libsrc/occ/occ_edge.hpp

@@ -0,0 +1,40 @@
+#ifndef FILE_OCC_EDGE_INCLUDED
+#define FILE_OCC_EDGE_INCLUDED
+
+#include <GProp_GProps.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Edge.hxx>
+#include <Geom_Curve.hxx>
+#include <BRep_TEdge.hxx>
+
+#include "occ_vertex.hpp"
+#include "meshing.hpp"
+
+namespace netgen
+{
+    class OCCEdge : public GeometryEdge 
+    {
+        public:
+        T_Shape tedge;
+        TopoDS_Edge edge;
+        Handle(Geom_Curve) curve;
+        double s0, s1;
+        GProp_GProps props;
+
+        public:
+        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;
+        Point<3> GetPoint(double t) const override;
+        double CalcStep(double t, double sag) const override;
+        size_t GetHash() const override;
+        void ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const override;
+        Vec<3> GetTangent(double t) const override;
+    };
+}
+
+#endif // FILE_OCCEDGE_INCLUDED

libsrc/occ/occ_face.cpp

@@ -0,0 +1,252 @@
+#include <BRepGProp.hxx>
+#include <BRep_Tool.hxx>
+#include <GeomAPI_ProjectPointOnCurve.hxx>
+
+#include "occ_edge.hpp"
+#include "occ_face.hpp"
+#include "occgeom.hpp"
+
+namespace netgen
+{
+    OCCFace::OCCFace(TopoDS_Shape dshape)
+        : tface(dshape.TShape()),
+          face(TopoDS::Face(dshape))
+    {
+        BRepGProp::LinearProperties(face, props);
+        bbox = ::netgen::GetBoundingBox(face);
+
+        surface = BRep_Tool::Surface(face);
+        shape_analysis = new ShapeAnalysis_Surface( surface );
+        tolerance = BRep_Tool::Tolerance( face );
+    }
+
+    size_t OCCFace::GetNBoundaries() const
+    {
+        return 0;
+    }
+
+    size_t OCCFace::GetHash() const
+    {
+        return reinterpret_cast<size_t>(tface.get());
+    }
+
+    Point<3> OCCFace::GetCenter() const
+    {
+        return occ2ng( props.CentreOfMass() );
+    }
+
+    Array<Segment> OCCFace::GetBoundary(const Mesh& mesh) const
+    {
+        auto & geom = dynamic_cast<OCCGeometry&>(*mesh.GetGeometry());
+
+        auto n_edges = geom.edge_map.size();
+        constexpr int UNUSED = 0;
+        constexpr int FORWARD = 1;
+        constexpr int REVERSED = 2;
+        constexpr int BOTH = 3;
+
+        Array<int> edge_orientation(n_edges);
+        edge_orientation = UNUSED;
+
+        Array<Handle(Geom2d_Curve)> curve_on_face[BOTH];
+        curve_on_face[FORWARD].SetSize(n_edges);
+        curve_on_face[REVERSED].SetSize(n_edges);
+
+        Array<TopoDS_Edge> edge_on_face[BOTH];
+        edge_on_face[FORWARD].SetSize(n_edges);
+        edge_on_face[REVERSED].SetSize(n_edges);
+
+        for(auto edge_ : GetEdges(face))
+        {
+            auto edge = TopoDS::Edge(edge_);
+            if(geom.edge_map.count(edge.TShape())==0)
+                continue;
+            auto edgenr = geom.edge_map[edge.TShape()];
+            auto & orientation = edge_orientation[edgenr];
+            double s0, s1;
+            auto cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1);
+            if(edge.Orientation() == TopAbs_FORWARD)
+            {
+                curve_on_face[FORWARD][edgenr] = cof;
+                orientation += FORWARD;
+                edge_on_face[FORWARD][edgenr] = edge;
+            }
+            if(edge.Orientation() == TopAbs_REVERSED)
+            {
+                curve_on_face[REVERSED][edgenr] = cof;
+                orientation += REVERSED;
+                edge_on_face[REVERSED][edgenr] = edge;
+            }
+
+            if(orientation > BOTH)
+                throw Exception("have edge more than twice in face " + ToString(nr) + " " + properties.GetName() + ", orientation: " + ToString(orientation));
+        }
+
+        Array<Segment> boundary;
+        for (auto seg : mesh.LineSegments())
+        {
+            auto edgenr = seg.epgeominfo[0].edgenr;
+            auto orientation = edge_orientation[edgenr];
+
+            if(orientation == UNUSED)
+                continue;
+
+            for(const auto ORIENTATION : {FORWARD, REVERSED})
+            {
+                if((orientation & ORIENTATION) == 0)
+                    continue;
+
+                // auto cof = curve_on_face[ORIENTATION][edgenr];
+                auto edge = edge_on_face[ORIENTATION][edgenr];
+                double s0, s1;
+                auto cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1);
+
+                double s[2] = { seg.epgeominfo[0].dist, seg.epgeominfo[1].dist };
+
+                // dist is in [0,1], map parametrization to [s0, s1]
+                s[0] = s0 + s[0]*(s1-s0);
+                s[1] = s0 + s[1]*(s1-s0);
+
+                // fixes normal-vector roundoff problem when endpoint is cone-tip
+                double delta = s[1]-s[0];
+                s[0] += 1e-10*delta;
+                s[1] -= 1e-10*delta;
+
+                for(auto i : Range(2))
+                {
+                    auto uv = cof->Value(s[i]);
+                    seg.epgeominfo[i].u = uv.X();
+                    seg.epgeominfo[i].v = uv.Y();
+                }
+
+                if(ORIENTATION == REVERSED)
+                {
+                    swap(seg[0], seg[1]);
+                    swap(seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
+                    swap(seg.epgeominfo[0].u, seg.epgeominfo[1].u);
+                    swap(seg.epgeominfo[0].v, seg.epgeominfo[1].v);
+                }
+
+                boundary.Append(seg);
+            }
+        }
+        return boundary;
+    }
+
+    PointGeomInfo OCCFace::Project(Point<3>& p) const
+    {
+        auto suval = shape_analysis->ValueOfUV(ng2occ(p), tolerance);
+        double u,v;
+        suval.Coord(u, v);
+        p = occ2ng(surface->Value( u, v ));
+
+        PointGeomInfo gi;
+        gi.trignum = nr+1;
+        gi.u = u;
+        gi.v = v;
+        return gi;
+    }
+
+    bool OCCFace::ProjectPointGI(Point<3>& p_, PointGeomInfo& gi) const
+    {
+        double u = gi.u;
+        double v = gi.v;
+        auto p = ng2occ(p_);
+        auto x = surface->Value (u,v);
+      
+        if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true;
+      
+        gp_Vec du, dv;
+        surface->D1(u,v,x,du,dv);
+      
+        int count = 0;
+        gp_Pnt xold;
+        gp_Vec n;
+        double det, lambda, mu;
+      
+        do {
+           count++;
+
+           n = du^dv;
+
+           det = Det3 (n.X(), du.X(), dv.X(),
+              n.Y(), du.Y(), dv.Y(),
+              n.Z(), du.Z(), dv.Z());
+
+           if (det < 1e-15) return false;
+
+           lambda = Det3 (n.X(), p.X()-x.X(), dv.X(),
+              n.Y(), p.Y()-x.Y(), dv.Y(),
+              n.Z(), p.Z()-x.Z(), dv.Z())/det;
+
+           mu     = Det3 (n.X(), du.X(), p.X()-x.X(),
+              n.Y(), du.Y(), p.Y()-x.Y(),
+              n.Z(), du.Z(), p.Z()-x.Z())/det;
+
+           u += lambda;
+           v += mu;
+
+           xold = x;
+           surface->D1(u,v,x,du,dv);
+
+        } while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50);
+
+        //    (*testout) << "FastProject count: " << count << endl;
+
+        if (count == 50) return false;
+
+        p_ = occ2ng(x);
+
+        return true;
+    }
+
+    Point<3> OCCFace::GetPoint(const PointGeomInfo& gi) const
+    {
+        return occ2ng(surface->Value( gi.u, gi.v ));
+    }
+
+    void OCCFace::CalcEdgePointGI(const GeometryEdge& edge,
+            double t,
+            EdgePointGeomInfo& egi) const
+    {
+        throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__));
+    }
+
+    Box<3> OCCFace::GetBoundingBox() const
+    {
+        return bbox;
+    }
+
+
+    double OCCFace::GetCurvature(const PointGeomInfo& gi) const
+    {
+        throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__));
+    }
+
+    void OCCFace::RestrictH(Mesh& mesh, const MeshingParameters& mparam) const
+    {
+        throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__));
+    }
+
+    Vec<3> OCCFace::GetNormal(const Point<3>& p, const PointGeomInfo* gi) const
+    {
+        PointGeomInfo gi_;
+        if(gi==nullptr)
+        {
+            auto p_ = p;
+            gi_ = Project(p_);
+            gi = &gi_;
+        }
+
+        gp_Pnt pnt;
+        gp_Vec du, dv;
+        surface->D1(gi->u,gi->v,pnt,du,dv);
+        auto n = Cross (occ2ng(du), occ2ng(dv));
+        n.Normalize();
+        if (face.Orientation() == TopAbs_REVERSED)
+            n *= -1;
+        return n;
+    }
+
+
+}

libsrc/occ/occ_face.hpp

@@ -0,0 +1,49 @@
+#ifndef FILE_OCC_FACE_INCLUDED
+#define FILE_OCC_FACE_INCLUDED
+
+#include <GProp_GProps.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Face.hxx>
+#include <ShapeAnalysis_Surface.hxx>
+
+#include "occ_vertex.hpp"
+#include "meshing.hpp"
+
+namespace netgen
+{
+    class OCCFace : public GeometryFace 
+    {
+        T_Shape tface;
+        TopoDS_Face face;
+        GProp_GProps props;
+        Box<3> bbox;
+
+        Handle( Geom_Surface ) surface;
+        Handle( ShapeAnalysis_Surface ) shape_analysis;
+        double tolerance;
+
+        public:
+        OCCFace(TopoDS_Shape dshape);
+
+        T_Shape TShape() { return tface; }
+
+        size_t GetHash() const override;
+        Point<3> GetCenter() const override;
+        virtual size_t GetNBoundaries() const override;
+        virtual Array<Segment> GetBoundary(const Mesh& mesh) const override;
+        virtual PointGeomInfo Project(Point<3>& p) const override;
+        virtual bool ProjectPointGI(Point<3>& p, PointGeomInfo& gi) const override;
+        virtual Point<3> GetPoint(const PointGeomInfo& gi) const override;
+        virtual void CalcEdgePointGI(const GeometryEdge& edge,
+                double t,
+                EdgePointGeomInfo& egi) const override;
+        virtual Box<3> GetBoundingBox() const override;
+
+        virtual double GetCurvature(const PointGeomInfo& gi) const override;
+
+        virtual void RestrictH(Mesh& mesh, const MeshingParameters& mparam) const override;
+        virtual Vec<3> GetNormal(const Point<3>& p, const PointGeomInfo* gi = nullptr) const override;
+    };
+}
+
+#endif // FILE_OCC_FACE_INCLUDED

libsrc/occ/occ_solid.hpp

@@ -0,0 +1,27 @@
+#ifndef FILE_OCC_SOLID_INCLUDED
+#define FILE_OCC_SOLID_INCLUDED
+
+#include <TopoDS.hxx>
+#include <TopoDS_Solid.hxx>
+
+#include "meshing.hpp"
+
+namespace netgen
+{
+    class OCCSolid : public GeometrySolid 
+    {
+        T_Shape tsolid;
+        TopoDS_Solid solid;
+
+        public:
+        OCCSolid(TopoDS_Shape dshape)
+            : tsolid(dshape.TShape()),
+              solid(TopoDS::Solid(dshape))
+        { }
+
+        T_Shape TShape() { return tsolid; }
+        size_t GetHash() const override { return reinterpret_cast<size_t>(tsolid.get()); }
+    };
+}
+
+#endif // FILE_OCC_SOLID_INCLUDED

libsrc/occ/occ_utils.cpp

@@ -0,0 +1,40 @@
+#include <Bnd_Box.hxx>
+#include <BRepBndLib.hxx>
+#include <BRep_TVertex.hxx>
+
+#include "occ_utils.hpp"
+
+namespace netgen
+{
+    Point<3> occ2ng (Handle(TopoDS_TShape) shape)
+    {
+        return occ2ng( Handle(BRep_TVertex)::DownCast(shape)->Pnt() );
+    }
+
+    Transformation<3> occ2ng (const gp_Trsf & occ_trafo)
+    {
+        Transformation<3> trafo;
+        auto v = occ_trafo.TranslationPart();
+        auto m = occ_trafo.VectorialPart();
+        auto & tv = trafo.GetVector();
+        auto & tm = trafo.GetMatrix();
+        for(auto i : Range(3))
+        {
+            tv[i] = v.Coord(i+1);
+            for(auto k : Range(3))
+                tm(i,k) = m(i+1,k+1);
+        }
+        return trafo;
+    }
+
+    Box<3> GetBoundingBox( const TopoDS_Shape & shape )
+    {
+        Bnd_Box bb;
+#if OCC_VERSION_HEX < 0x070000
+        BRepBndLib::Add (shape, bb);
+#else
+        BRepBndLib::Add (shape, bb, true);
+#endif
+        return {occ2ng(bb.CornerMin()), occ2ng(bb.CornerMax())};
+    }
+}

libsrc/occ/occ_utils.hpp

@@ -0,0 +1,283 @@
+#ifndef FILE_OCC_UTILS_INCLUDED
+#define FILE_OCC_UTILS_INCLUDED
+
+#include <BRepGProp.hxx>
+#include <BRep_Tool.hxx>
+#include <GProp_GProps.hxx>
+#include <Standard_Version.hxx>
+#include <TopExp_Explorer.hxx>
+#include <TopTools_IndexedMapOfShape.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Vertex.hxx>
+#include <gp_Trsf.hxx>
+
+#include "meshing.hpp"
+
+#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4
+#define OCC_HAVE_DUMP_JSON
+#endif
+
+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());
+    }
+
+    inline Point<2> occ2ng (const gp_Pnt2d & p)
+    {
+        return Point<2> (p.X(), p.Y());
+    }
+
+    inline Vec<3> occ2ng (const gp_Vec & v)
+    {
+        return Vec<3> (v.X(), v.Y(), v.Z());
+    }
+
+    DLL_HEADER Point<3> occ2ng (T_Shape shape);
+
+    inline Point<3> occ2ng (const TopoDS_Shape & s)
+    {
+        return occ2ng(s.TShape());
+    }
+
+    inline Point<3> occ2ng (const TopoDS_Vertex & v)
+    {
+        return occ2ng (BRep_Tool::Pnt (v));
+    }
+
+    DLL_HEADER Transformation<3> occ2ng (const gp_Trsf & t);
+
+    inline gp_Pnt ng2occ (const Point<3> & p)
+    {
+        return gp_Pnt(p(0), p(1), p(2));
+    }
+
+    DLL_HEADER Box<3> GetBoundingBox( const TopoDS_Shape & shape );
+
+    class OCCIdentification
+    {
+    public:
+      T_Shape from;
+      T_Shape to;
+      Transformation<3> trafo;
+      string name;
+      Identifications::ID_TYPE type;
+      bool opposite_direction;
+    };
+
+
+    class MyExplorer
+    {
+      class Iterator
+      {
+        TopExp_Explorer exp;
+      public:
+        Iterator (TopoDS_Shape ashape, TopAbs_ShapeEnum atoFind, TopAbs_ShapeEnum atoAvoid)
+          : exp(ashape, atoFind, atoAvoid) { }
+        auto operator*() { return exp.Current(); }
+        Iterator & operator++() { exp.Next(); return *this; }
+        bool operator!= (nullptr_t nu) { return exp.More(); }
+      };
+
+    public:
+      TopoDS_Shape shape;
+      TopAbs_ShapeEnum toFind;
+      TopAbs_ShapeEnum toAvoid;
+      MyExplorer (TopoDS_Shape ashape, TopAbs_ShapeEnum atoFind, TopAbs_ShapeEnum atoAvoid = TopAbs_SHAPE)
+        : shape(ashape), toFind(atoFind), toAvoid(atoAvoid) { ; }
+      Iterator begin() { return Iterator(shape, toFind, toAvoid); }
+      auto end() { return nullptr; }
+    };
+
+    inline auto Explore (TopoDS_Shape shape, TopAbs_ShapeEnum toFind, TopAbs_ShapeEnum toAvoid = TopAbs_SHAPE)
+    {
+      return MyExplorer (shape, toFind, toAvoid);
+    }
+
+
+    class IndexMapIterator
+    {
+      class Iterator
+      {
+        const TopTools_IndexedMapOfShape & indmap;
+        int i;
+      public:
+        Iterator (const TopTools_IndexedMapOfShape & aindmap, int ai)
+          : indmap(aindmap), i(ai) { ; }
+        auto operator*() { return tuple(i, indmap(i)); }
+        Iterator & operator++() { i++; return *this; }
+        bool operator!= (const Iterator & i2) { return i != i2.i; }
+      };
+
+    public:
+      const TopTools_IndexedMapOfShape & indmap;
+      IndexMapIterator (const TopTools_IndexedMapOfShape & aindmap) : indmap(aindmap) { }
+      Iterator begin() { return Iterator(indmap, 1); }
+      Iterator end() { return Iterator(indmap, indmap.Extent()+1); }
+    };
+
+    inline auto Enumerate (const TopTools_IndexedMapOfShape & indmap)
+    {
+      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:
+      DLL_HEADER TopoDS_Shape Max(gp_Vec dir);
+      DLL_HEADER TopoDS_Shape Nearest(gp_Pnt pnt);
+      DLL_HEADER ListOfShapes SubShapes(TopAbs_ShapeEnum type) const;
+
+      ListOfShapes Solids() const
+      {
+        return SubShapes(TopAbs_SOLID);
+      }
+      ListOfShapes Faces() const
+      {
+        return SubShapes(TopAbs_FACE);
+      }
+      ListOfShapes Wires() const
+      {
+        return SubShapes(TopAbs_WIRE);
+      }
+      ListOfShapes Edges() const
+      {
+        return SubShapes(TopAbs_EDGE);
+      }
+      ListOfShapes Vertices() const
+      {
+        return SubShapes(TopAbs_VERTEX);
+      }
+
+      ListOfShapes operator*(const ListOfShapes& other) const
+      {
+        ListOfShapes common;
+        for(const auto& shape : (*this))
+          for(const auto& shape_o : other)
+            if(shape.IsSame(shape_o))
+              common.push_back(shape);
+        return common;
+      }
+    };
+
+    inline ListOfShapes GetSolids(const TopoDS_Shape & shape)
+    {
+      ListOfShapes sub;
+      for (TopExp_Explorer e(shape, TopAbs_SOLID); e.More(); e.Next())
+        sub.push_back(e.Current());
+      return sub;
+    }
+
+    inline ListOfShapes GetFaces(const TopoDS_Shape & shape)
+    {
+      ListOfShapes sub;
+      for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next())
+        sub.push_back(e.Current());
+      return sub;
+    }
+
+    inline ListOfShapes GetWires(const TopoDS_Shape & shape)
+    {
+      ListOfShapes sub;
+      for (TopExp_Explorer e(shape, TopAbs_WIRE); e.More(); e.Next())
+        sub.push_back(e.Current());
+      return sub;
+    }
+
+    inline ListOfShapes GetEdges(const TopoDS_Shape & shape)
+    {
+      ListOfShapes sub;
+      for (TopExp_Explorer e(shape, TopAbs_EDGE); e.More(); e.Next())
+        sub.push_back(e.Current());
+      return sub;
+    }
+
+    inline ListOfShapes GetVertices(const TopoDS_Shape & shape)
+    {
+      ListOfShapes sub;
+      for (TopExp_Explorer e(shape, TopAbs_VERTEX); e.More(); e.Next())
+        sub.push_back(e.Current());
+      return sub;
+    }
+
+  class DirectionalInterval
+  {
+  public:
+    gp_Vec dir;
+    double minval = -1e99;
+    double maxval = 1e99;
+    bool openmin = false, openmax = false;
+
+    DirectionalInterval (gp_Vec adir) : dir(adir) { ; }
+    DirectionalInterval (const DirectionalInterval & i2)
+      : dir(i2.dir), minval(i2.minval), maxval(i2.maxval) { ; }
+
+    DirectionalInterval operator< (double val) const
+    {
+      DirectionalInterval i2 = *this;
+      i2.maxval = val;
+      return i2;
+    }
+
+    DirectionalInterval operator> (double val) const
+    {
+      DirectionalInterval i2 = *this;
+      i2.minval = val;
+      return i2;
+    }
+
+
+    DirectionalInterval Intersect (const DirectionalInterval & i2)
+    {
+      DirectionalInterval res = *this;
+      res.minval = max(res.minval, i2.minval);
+      res.maxval = min(res.maxval, i2.maxval);
+      return res;
+    }
+
+    bool Contains (gp_Pnt p, double eps = 1e-8)
+    {
+      // cout << "Contains point " << p.X() << "," << p.Y() << "," << p.Z() << " ? " << endl;
+      double val = dir.X()*p.X() + dir.Y()*p.Y() + dir.Z() * p.Z();
+      // cout << "minval = " << minval << ", val = " << val << " maxval = " << maxval << endl;
+      if (openmin) {
+        if (val < minval+eps) return false;
+      } else {
+        if (val < minval-eps) return false;
+      }
+      if (openmax) {
+        if (val > maxval-eps) return false;
+      } else {
+        if (val > maxval+eps) return false;
+      }
+      return true;
+    }
+  };
+
+
+  inline gp_Pnt Center (TopoDS_Shape shape)
+  {
+    GProp_GProps props;
+    switch (shape.ShapeType())
+      {
+      case TopAbs_FACE:
+        BRepGProp::SurfaceProperties (shape, props); break;
+      default:
+        BRepGProp::LinearProperties(shape, props);
+      }
+    return props.CentreOfMass();
+  }
+
+}
+#endif // FILE_OCC_UTILS_INCLUDED

libsrc/occ/occ_vertex.cpp

@@ -0,0 +1,25 @@
+#include <BRepGProp.hxx>
+#include <BRep_Tool.hxx>
+
+#include "occ_vertex.hpp"
+
+namespace netgen
+{
+
+    OCCVertex::OCCVertex( TopoDS_Shape s )
+        : vertex(TopoDS::Vertex(s)),
+          tvertex(s.TShape())
+    {
+        p = occ2ng(vertex);
+    }
+
+    Point<3> OCCVertex::GetPoint() const
+    {
+        return p;
+    }
+
+    size_t OCCVertex::GetHash() const
+    {
+        return reinterpret_cast<size_t>(tvertex.get());
+    }
+}

libsrc/occ/occ_vertex.hpp

@@ -0,0 +1,28 @@
+#ifndef FILE_OCC_VERTEX_INCLUDED
+#define FILE_OCC_VERTEX_INCLUDED
+
+#include <TopoDS.hxx>
+#include <BRep_TVertex.hxx>
+
+#include "meshing.hpp"
+#include "occ_utils.hpp"
+
+namespace netgen
+{
+    class OCCVertex : public GeometryVertex
+    {
+        TopoDS_Vertex vertex;
+        T_Shape tvertex;
+        Point<3> p;
+
+        public:
+        OCCVertex( ) = default;
+        OCCVertex( TopoDS_Shape s );
+        ~OCCVertex() {}
+        Point<3> GetPoint() const override;
+        size_t GetHash() const override;
+        T_Shape TShape() { return tvertex; }
+    };
+}
+
+#endif // FILE_OCC_VERTEX_INCLUDED

libsrc/occ/occgenmesh.cpp

@@ -1,14 +1,27 @@
 #ifdef OCCGEOMETRY
 
 #include <mystdlib.h>
-#include <occgeom.hpp>
 #include <meshing.hpp>
 
+#include "occgeom.hpp"
+#include "occmeshsurf.hpp"
+
+#include <BRepAdaptor_Curve.hxx>
+#include <BRepGProp.hxx>
+#include <BRepLProp_CLProps.hxx>
+#include <BRepLProp_SLProps.hxx>
+#include <BRepMesh_IncrementalMesh.hxx>
+#include <BRepTools.hxx>
+#include <GProp_GProps.hxx>
+#include <Quantity_Color.hxx>
+#include <ShapeAnalysis.hxx>
+#include <TopExp.hxx>
+#include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
+#include <TopoDS_Edge.hxx>
 
 namespace netgen
 {
 
-#include "occmeshsurf.hpp"
 
 #define TCL_OK 0
 #define TCL_ERROR 1
@@ -218,536 +231,27 @@ namespace netgen
       }
   }
 
-
+  bool OCCMeshFace (const OCCGeometry & geom, Mesh & mesh, FlatArray<int, PointIndex> glob2loc,
-
+                       const MeshingParameters & mparam, int nr, int projecttype, bool delete_on_failure)
-  void DivideEdge (TopoDS_Edge & edge, NgArray<MeshPoint> & ps,
-                   Array<double> & params, Mesh & mesh,
-                   const MeshingParameters & mparam)
-  {
-    double s0, s1;
-    int nsubedges = 1;
-    gp_Pnt pnt, oldpnt;
-    double svalue[DIVIDEEDGESECTIONS];
-
-    GProp_GProps system;
-    BRepGProp::LinearProperties(edge, system);
-    double L = system.Mass();
-
-    Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1);
-
-    double hvalue[DIVIDEEDGESECTIONS+1];
-    hvalue[0] = 0;
-    pnt = c->Value(s0);
-
-    int tmpVal = (int)(DIVIDEEDGESECTIONS);
-
-    for (int i = 1; i <= tmpVal; i++)
-      {
-        oldpnt = pnt;
-        pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0));
-        hvalue[i] = hvalue[i-1] +
-          1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*
-          pnt.Distance(oldpnt);
-
-        //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))
-        //	   <<  " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl;
-      }
-
-    //  nsubedges = int(ceil(hvalue[DIVIDEEDGESECTIONS]));
-    nsubedges = max (1, int(floor(hvalue[DIVIDEEDGESECTIONS]+0.5)));
-
-    ps.SetSize(nsubedges-1);
-    params.SetSize(nsubedges+1);
-
-    int i = 1;
-    int i1 = 0;
-    do
-      {
-        if (hvalue[i1]/hvalue[DIVIDEEDGESECTIONS]*nsubedges >= i)
-          {
-            params[i] = s0+(i1/double(DIVIDEEDGESECTIONS))*(s1-s0);
-            pnt = c->Value(params[i]);
-            ps[i-1] = MeshPoint (Point3d(pnt.X(), pnt.Y(), pnt.Z()));
-            i++;
-          }
-        i1++;
-        if (i1 > DIVIDEEDGESECTIONS)
-          {
-            nsubedges = i;
-            ps.SetSize(nsubedges-1);
-            params.SetSize(nsubedges+1);
-            cout << "divide edge: local h too small" << endl;
-          }
-      } while (i < nsubedges);
-
-    params[0] = s0;
-    params[nsubedges] = s1;
-
-    if (params[nsubedges] <= params[nsubedges-1])
-      {
-        cout << "CORRECTED" << endl;
-        ps.SetSize (nsubedges-2);
-        params.SetSize (nsubedges);
-        params[nsubedges] = s1;
-      }
-  }
-
-
-
-
-  void OCCFindEdges (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam)
-  {
-    static Timer t("OCCFindEdges"); RegionTimer r(t);
-    static Timer tsearch("OCCFindEdges - search point");
-    const char * savetask = multithread.task;
-    multithread.task = "Edge meshing";
-
-    (*testout) << "edge meshing" << endl;
-
-    int nvertices = geom.vmap.Extent();
-    int nedges = geom.emap.Extent();
-
-    Array<Array<PointIndex>> alledgepnums(nedges);
-    Array<Array<double>> alledgeparams(nedges);     
-    
-    (*testout) << "nvertices = " << nvertices << endl;
-    (*testout) << "nedges = " << nedges << endl;
-
-    double eps = 1e-6 * geom.GetBoundingBox().Diam();
-
-    tsearch.Start();
-    for (auto [i,vshape] : Enumerate(geom.vmap))
-      {
-        TopoDS_Vertex vertex = TopoDS::Vertex(vshape);
-        gp_Pnt pnt = BRep_Tool::Pnt (vertex); 
-
-        mesh.AddPoint (occ2ng(pnt));
-        
-        double hpref = OCCGeometry::global_shape_properties[vertex.TShape()].hpref;
-        mesh.Points().Last().Singularity(hpref);
-        
-        double maxh = OCCGeometry::global_shape_properties[vertex.TShape()].maxh;
-        mesh.RestrictLocalH (occ2ng(pnt), maxh);
-      }
-    tsearch.Stop();
-    (*testout) << "different vertices = " << mesh.GetNP() << endl;
-
-    // int first_ep = mesh.GetNP()+1;
-    // PointIndex first_ep = mesh.Points().End();
-    PointIndex first_ep = *mesh.Points().Range().end();
-    auto vertexrange = mesh.Points().Range();
-
-    NgArray<int> face2solid[2];
-    for (int i = 0; i < 2; i++)
-      {
-        face2solid[i].SetSize (geom.fmap.Extent());
-        face2solid[i] = 0;
-      }
-
-    /*
-    int solidnr = 0;
-    for (TopExp_Explorer exp0(geom.shape, TopAbs_SOLID); exp0.More(); exp0.Next())
-      {
-        solidnr++;
-        for (TopExp_Explorer exp1(exp0.Current(), TopAbs_FACE); exp1.More(); exp1.Next())
-          {
-            TopoDS_Face face = TopoDS::Face(exp1.Current());
-            int facenr = geom.fmap.FindIndex(face);
-            if(facenr < 1) continue;
-
-            if (face2solid[0][facenr-1] == 0)
-              face2solid[0][facenr-1] = solidnr;
-            else
-              face2solid[1][facenr-1] = solidnr;
-          }
-      }
-    */
-    int solidnr = 0;
-    for (auto solid : Explore(geom.shape, TopAbs_SOLID))
-      {
-        solidnr++;
-        for (auto face : Explore (solid, TopAbs_FACE))
-          if (geom.fmap.Contains(face))
-            {
-              int facenr = geom.fmap.FindIndex(face);
-              if (face2solid[0][facenr-1] == 0)
-                face2solid[0][facenr-1] = solidnr;
-              else
-                face2solid[1][facenr-1] = solidnr;
-            }
-      }
-    
-
-
-    /*
-    int total = 0;
-    for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++)
-      for (TopExp_Explorer exp2(geom.fmap(i3), TopAbs_WIRE); exp2.More(); exp2.Next())
-        for (TopExp_Explorer exp3(exp2.Current(), TopAbs_EDGE); exp3.More(); exp3.Next())
-          total++;
-    */
-    int total = 0;
-    for (auto [i3, face] : Enumerate(geom.fmap))
-      for (auto wire : Explore(face, TopAbs_WIRE))
-        for (auto edge : Explore(wire, TopAbs_EDGE))
-          total++;
-      
-    
-
-    int facenr = 0;
-    // int edgenr = mesh.GetNSeg();
-
-    (*testout) << "faces = " << geom.fmap.Extent() << endl;
-    int curr = 0;
-
-    /*
-    for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++)
-      {
-        TopoDS_Face face = TopoDS::Face(geom.fmap(i3));
-    */
-    for (auto [i3,faceshape] : Enumerate(geom.fmap))
-      {
-        TopoDS_Face face = TopoDS::Face(faceshape);
-        facenr = geom.fmap.FindIndex (face);       // sollte doch immer == i3 sein ??? JS
-        if (facenr != i3)
-          cout << "info: facenr != i3, no problem, but please report to developers" << endl;
-        
-        int solidnr0 = face2solid[0][i3-1];
-        int solidnr1 = face2solid[1][i3-1];
-
-        /* auskommentiert am 3.3.05 von robert
-           for (exp2.Init (geom.somap(solidnr0), TopAbs_FACE); exp2.More(); exp2.Next())
-           {
-           TopoDS_Face face2 = TopoDS::Face(exp2.Current());
-           if (geom.fmap.FindIndex(face2) == facenr)
-           {
-           //		      if (face.Orientation() != face2.Orientation()) swap (solidnr0, solidnr1);
-           }
-           }
-        */
-
-        mesh.AddFaceDescriptor (FaceDescriptor(facenr, solidnr0, solidnr1, 0));
-
-        Vec<4> col = OCCGeometry::global_shape_properties[face.TShape()].col.value_or(Vec<4>(0,1,0,1));
-        mesh.GetFaceDescriptor(facenr).SetSurfColour(col);
-
-        if(auto & opt_name = geom.fprops[facenr-1]->name)
-          mesh.GetFaceDescriptor(facenr).SetBCName(*opt_name);
-        else
-          mesh.GetFaceDescriptor(facenr).SetBCName("bc_"+ToString(facenr));
-        mesh.GetFaceDescriptor(facenr).SetBCProperty(facenr);
-        // ACHTUNG! STIMMT NICHT ALLGEMEIN (RG)
-        // kA was RG damit meinte
-
-
-        Handle(Geom_Surface) occface = BRep_Tool::Surface(face);
-
-        /*
-        for (TopExp_Explorer exp2 (face, TopAbs_WIRE); exp2.More(); exp2.Next())
-          {
-            TopoDS_Shape wire = exp2.Current();
-        */
-        for (auto wire : MyExplorer (face, TopAbs_WIRE))
-          {
-            // for (TopExp_Explorer exp3 (wire, TopAbs_EDGE); exp3.More(); exp3.Next())
-            for (auto edgeshape : MyExplorer (wire, TopAbs_EDGE))
-              {
-                TopoDS_Edge edge = TopoDS::Edge(edgeshape);
-                curr++;
-                (*testout) << "edge nr " << curr << endl;
-                multithread.percent = 100 * curr / double (total);
-                if (multithread.terminate) return;
-
-                // TopoDS_Edge edge = TopoDS::Edge (exp3.Current());
-                if (BRep_Tool::Degenerated(edge))
-                  {
-                    //(*testout) << "ignoring degenerated edge" << endl;
-                    continue;
-                  }
-
-                if (!geom.emap.Contains(edge)) continue;
-
-                if (geom.vmap.FindIndex(TopExp::FirstVertex (edge)) ==
-                    geom.vmap.FindIndex(TopExp::LastVertex (edge)))
-                  {
-                    GProp_GProps system;
-                    BRepGProp::LinearProperties(edge, system);
-
-                    if (system.Mass() < eps)
-                      {
-                        cout << "ignoring edge " << geom.emap.FindIndex (edge)
-                             << ". closed edge with length < " << eps << endl;
-                        continue;
-                      }
-                  }
-
-                double s0, s1;
-                Handle(Geom2d_Curve) cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1);
-
-                int geomedgenr = geom.emap.FindIndex(edge);
-                Array<PointIndex> pnums;
-                Array<double> params;
-                
-                
-                // check for identifications
-                bool copy_identified = false;
-                if (auto it = geom.identifications.find(edge.TShape()); it != geom.identifications.end())
-                  for (auto & ids : it->second)
-                    {
-                      cout << "edge has identification with trafo " << ids.name << ", inv = " << ids.inverse << endl;
-                      int otherind = geom.emap.FindIndex(ids.other);
-                      Array<Segment> othersegs;
-                      for (auto & seg : mesh.LineSegments())
-                        if (seg.edgenr == otherind)
-                          othersegs.Append (seg);
-
-                      if (othersegs.Size())
-                        {
-                          cout << "other has already segs" << endl;
-                          copy_identified = true;
-                                                              
-                          Array<PointIndex> pnums_other;
-                          pnums_other.Append (othersegs[0][0]);
-                          for (auto & seg : othersegs)
-                            pnums_other.Append (seg[1]);
-
-                          auto inv = ids.trafo.CalcInverse();
-                          // for (auto & pi : pnums)
-                          for (auto oi : Range(pnums_other))
-                            {
-                              PointIndex piother = pnums_other[pnums_other.Size()-oi-1];
-                              Point<3> pother = mesh[piother];
-                              Point<3> p = inv(pother);
-                              
-                              bool found = false;
-                              PointIndex pi;
-                              for (PointIndex piv : vertexrange)
-                                if (Dist2 (mesh[piv], p) < eps*eps)
-                                  {
-                                    pi = piv;
-                                    found = true;
-                                  }
-                              
-                              if (!found)
-                                pi = mesh.AddPoint (p);
-                              
-                              // params.Add ( find parameter p );
-                              double s0, s1;
-                              Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, s0, s1);
-                              
-                              GeomAPI_ProjectPointOnCurve proj(ng2occ(p), curve);
-                              params.Append (proj.LowerDistanceParameter());
-                              pnums.Append (pi);
-                              mesh.GetIdentifications().Add (pi, piother, geomedgenr);
-                              
-                            }
-                          mesh.GetIdentifications().SetType(geomedgenr,Identifications::PERIODIC);
-                          
-                          copy_identified = true;
-                          break;
-                        }
-                    }
-
-                
-                if (!copy_identified)
-                  {
-
-                    
-                    if (alledgepnums[geomedgenr-1].Size())
-                      {
-                        pnums = alledgepnums[geomedgenr-1];
-                        params = alledgeparams[geomedgenr-1];
-                      }
-                    else
-                      {
-                        NgArray <MeshPoint> mp;
-                        DivideEdge (edge, mp, params, mesh, mparam);
-                        
-                        pnums.SetSize(mp.Size()+2);
-                        if (!merge_solids)
-                          {
-                            pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1;
-                            pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1;
-                          }
-                        else
-                          {
-                            Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge)));
-                            Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge)));
-                            
-                            pnums[0] = PointIndex::INVALID;
-                            pnums.Last() = PointIndex::INVALID;
-                            for (PointIndex pi : vertexrange)
-                              {
-                                if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi;
-                                if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi;
-                              }
-                          }
-                        
-                        for (size_t i = 1; i <= mp.Size(); i++)
-                          {
-                            bool exists = false;
-                            tsearch.Start();
-
-                            /*
-                            // for (PointIndex j = first_ep; j < mesh.Points().End(); j++)
-                            for (PointIndex j = first_ep; j < *mesh.Points().Range().end(); j++)
-                              if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps)
-                                {
-                                  exists = true;
-                                  pnums[i] = j;
-                                  break;
-                                }
-                            */
-                            tsearch.Stop();
-                            
-                            if (!exists)
-                              pnums[i] = mesh.AddPoint (mp[i-1]);
-                          }
-                      }
-                    
-
-                    alledgepnums[geomedgenr-1] = pnums;
-                    alledgeparams[geomedgenr-1] = params;
-                  }
-
-                auto name = geom.eprops[geom.emap.FindIndex(edge)-1]->name.value_or("");
-                mesh.SetCD2Name(geomedgenr, name);
-                  
-                (*testout) << "NP = " << mesh.GetNP() << endl;
-                //(*testout) << pnums[pnums.Size()-1] << endl;
-
-                double hpref = OCCGeometry::global_shape_properties[edge.TShape()].hpref;
-                
-                // for (size_t i = 1; i <= mp.Size()+1; i++)
-                for (size_t i = 1; i < pnums.Size(); i++)
-                  {
-                    // edgenr++;
-                    Segment seg;
-
-                    seg[0] = pnums[i-1];
-                    seg[1] = pnums[i];
-                    // seg.edgenr = edgenr;
-                    seg.edgenr = geomedgenr;
-                    seg.si = facenr;
-                    seg.epgeominfo[0].dist = params[i-1];
-                    seg.epgeominfo[1].dist = params[i];
-                    seg.epgeominfo[0].edgenr = geomedgenr;
-                    seg.epgeominfo[1].edgenr = geomedgenr;
-
-                    double s0 = params[i-1];
-                    double s1 = params[i];
-                    double delta = s1-s0;
-                    s0 += 1e-10*delta;   // fixes normal-vector roundoff problem when endpoint is cone-tip
-                    s1 -= 1e-10*delta;  
-                    gp_Pnt2d p2d1 = cof->Value(s0);
-                    gp_Pnt2d p2d2 = cof->Value(s1);
-                    seg.epgeominfo[0].u = p2d1.X(); 
-                    seg.epgeominfo[0].v = p2d1.Y(); 
-                    seg.epgeominfo[1].u = p2d2.X(); 
-                    seg.epgeominfo[1].v = p2d2.Y(); 
-
-                    seg.singedge_left = hpref;
-                    seg.singedge_right = hpref;
-                    /*
-                      if (occface->IsUPeriodic())
-                      {
-                      cout << "U Periodic" << endl;
-                      if (fabs(seg.epgeominfo[1].u-seg.epgeominfo[0].u) >
-                      fabs(seg.epgeominfo[1].u-
-                      (seg.epgeominfo[0].u-occface->UPeriod())))
-                      seg.epgeominfo[0].u = p2d.X()+occface->UPeriod();
-
-                      if (fabs(seg.epgeominfo[1].u-seg.epgeominfo[0].u) >
-                      fabs(seg.epgeominfo[1].u-
-                      (seg.epgeominfo[0].u+occface->UPeriod())))
-                      seg.epgeominfo[0].u = p2d.X()-occface->UPeriod();
-                      }
-
-                      if (occface->IsVPeriodic())
-                      {
-                      cout << "V Periodic" << endl;
-                      if (fabs(seg.epgeominfo[1].v-seg.epgeominfo[0].v) >
-                      fabs(seg.epgeominfo[1].v-
-                      (seg.epgeominfo[0].v-occface->VPeriod())))
-                      seg.epgeominfo[0].v = p2d.Y()+occface->VPeriod();
-
-                      if (fabs(seg.epgeominfo[1].v-seg.epgeominfo[0].v) >
-                      fabs(seg.epgeominfo[1].v-
-                      (seg.epgeominfo[0].v+occface->VPeriod())))
-                      seg.epgeominfo[0].v = p2d.Y()-occface->VPeriod();
-                      }
-                    */
-
-                    if (edge.Orientation() == TopAbs_REVERSED)
-                      {
-                        swap (seg[0], seg[1]);
-                        swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
-                        swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
-                        swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
-                      }
-
-                    mesh.AddSegment (seg);
-
-                    //edgesegments[geomedgenr-1]->Append(mesh.GetNSeg());
-
-                  }
-              }
-          }
-      }
-
-    //	for(i=1; i<=mesh.GetNSeg(); i++)
-    //		(*testout) << "edge " << mesh.LineSegment(i).edgenr << " face " << mesh.LineSegment(i).si
-    //				<< " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl;
-    //	exit(10);
-
-    mesh.CalcSurfacesOfNode();
-    multithread.task = savetask;
-  }
-
-
-
-
-  void OCCMeshSurface (const OCCGeometry & geom, Mesh & mesh,
-                       const MeshingParameters & mparam)
-  {
-    static Timer t("OCCMeshSurface"); RegionTimer r(t);
-    
-    // int i, j, k;
-    // int changed;
-
-    const char * savetask = multithread.task;
-    multithread.task = "Surface meshing";
-
-    geom.facemeshstatus = 0;
-
-    int noldp = mesh.GetNP();
-
-    double starttime = GetTime();
-
-    NgArray<int, PointIndex::BASE> glob2loc(noldp);
-
-    //int projecttype = PARAMETERSPACE;
-
-    int projecttype = PARAMETERSPACE;
-
-    int notrys = 1;
-
-    int surfmesherror = 0;
-
-    for (int k = 1; k <= mesh.GetNFD(); k++)
   {
+    auto k = nr+1;
     if(1==0 && !geom.fvispar[k-1].IsDrawable())
       {
         (*testout) << "ignoring face " << k << endl;
         cout << "ignoring face " << k << endl;
-            continue;
+        return true;
       }
 
+    // if(master_faces[k]!=k)
+    //     continue;
+
     (*testout) << "mesh face " << k << endl;
     multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL);
     geom.facemeshstatus[k-1] = -1;
 
     FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
+    auto face = TopoDS::Face(geom.fmap(k));
+    auto fshape = face.TShape();
 
     int oldnf = mesh.GetNSE();
 
@@ -758,33 +262,33 @@ namespace netgen
     
     static Timer tinit("init");
     tinit.Start();
-        Meshing2OCCSurfaces meshing(geom, TopoDS::Face(geom.fmap(k)), bb, projecttype, mparam);
+    Meshing2OCCSurfaces meshing(geom, face, bb, projecttype, mparam);
     tinit.Stop();
 
 
     static Timer tprint("print");
     tprint.Start();
     if (meshing.GetProjectionType() == PLANESPACE)
-          PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (plane space projection)");
+      PrintMessage (2, "Face ", k, " / ", geom.GetNFaces(), " (plane space projection)");
     else
-          PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (parameter space projection)");
+      PrintMessage (2, "Face ", k, " / ", geom.GetNFaces(), " (parameter space projection)");
     tprint.Stop();
-        if (surfmesherror)
-          cout << "Surface meshing error occurred before (in " << surfmesherror << " faces)" << endl;
 
     //      Meshing2OCCSurfaces meshing(f2, bb);
-        meshing.SetStartTime (starttime);
+    // meshing.SetStartTime (starttime);
     //(*testout) << "Face " << k << endl << endl;
     
 
+    auto segments = geom.GetFace(k-1).GetBoundary(mesh);
+
     if (meshing.GetProjectionType() == PLANESPACE)
       {
         static Timer t("MeshSurface: Find edges and points - Physical"); RegionTimer r(t);
         int cntp = 0;
         glob2loc = 0;
 
-            for (Segment & seg : mesh.LineSegments())
+        for (Segment & seg : segments)
-              if (seg.si == k)
+          // if (seg.si == k)
             for (int j = 0; j < 2; j++)
               {
                 PointIndex pi = seg[j];
@@ -801,8 +305,9 @@ namespace netgen
           {
             Segment & seg = mesh.LineSegment(i);
         */
-            for (Segment & seg : mesh.LineSegments())                
+        // for (Segment & seg : mesh.LineSegments())                
-              if (seg.si == k)
+        for (Segment & seg : segments)
+          //if (seg.si == k)
             {
               PointGeomInfo gi0, gi1;
               gi0.trignum = gi1.trignum = k;
@@ -811,93 +316,58 @@ namespace netgen
               gi1.u = seg.epgeominfo[1].u;
               gi1.v = seg.epgeominfo[1].v;
               
+              //if(orientation & 1)
               meshing.AddBoundaryElement (glob2loc[seg[0]], glob2loc[seg[1]], gi0, gi1);
+
             }
       }
     else
       {
         static Timer t("MeshSurface: Find edges and points - Parameter"); RegionTimer r(t);
         
-            int cntp = 0;
+        Array<PointGeomInfo> gis(2*segments.Size());
-
-            /*
-            for (int i = 1; i <= mesh.GetNSeg(); i++)
-              if (mesh.LineSegment(i).si == k)
-                cntp+=2;
-            */
-            for (Segment & seg : mesh.LineSegments())
-              if (seg.si == k)
-                cntp += 2;
-
-            NgArray<PointGeomInfo> gis;
-
-            gis.SetAllocSize (cntp);
         gis.SetSize (0);
 
-            for (int i = 1; i <= mesh.GetNSeg(); i++)
+        Box<2> uv_box(Box<2>::EMPTY_BOX);
+        for(auto & seg : segments)
+            for(auto i : Range(2))
+                uv_box.Add( {seg.epgeominfo[i].u, seg.epgeominfo[i].v } );
+
+        BoxTree<2> uv_tree(uv_box);
+        Array<int> found_points;
+
+        for(auto & seg : segments)
         {
-                Segment & seg = mesh.LineSegment(i);
+            PointGeomInfo gi[2];
-                if (seg.si == k)
+            gi[0].trignum = gi[1].trignum = k;
-                  {
+            gi[0].u = seg.epgeominfo[0].u;
-                    PointGeomInfo gi0, gi1;
+            gi[0].v = seg.epgeominfo[0].v;
-                    gi0.trignum = gi1.trignum = k;
+            gi[1].u = seg.epgeominfo[1].u;
-                    gi0.u = seg.epgeominfo[0].u;
+            gi[1].v = seg.epgeominfo[1].v;
-                    gi0.v = seg.epgeominfo[0].v;
-                    gi1.u = seg.epgeominfo[1].u;
-                    gi1.v = seg.epgeominfo[1].v;
 
             int locpnum[2] = {0, 0};
 
             for (int j = 0; j < 2; j++)
             {
-                        PointGeomInfo gi = (j == 0) ? gi0 : gi1;
+                Point<2> uv = {gi[j].u, gi[j].v};
+                uv_tree.GetIntersecting(uv, uv, found_points);
 
-                        /*
+                if(found_points.Size())
-                        int l;
+                    locpnum[j] = found_points[0];
-                        for (l = 0; l < gis.Size() && locpnum[j] == 0; l++)
+                else
-                          {
-                            double dist = sqr (gis[l].u-gi.u)+sqr(gis[l].v-gi.v);
-
-                            if (dist < 1e-10)
-                              locpnum[j] = l+1;
-                          }
-
-                        if (locpnum[j] == 0)
                 {
                     PointIndex pi = seg[j];
                     meshing.AddPoint (mesh.Point(pi), pi);
 
-                            gis.SetSize (gis.Size()+1);
+                    gis.Append (gi[j]);
-                            gis[l] = gi;
-                            locpnum[j] = l+1;
-                          }
-                        */
-                        for (int l = 0; l < gis.Size(); l++)
-                          {
-                            double dist = sqr (gis[l].u-gi.u)+sqr(gis[l].v-gi.v);
-                            if (dist < 1e-10)
-                              {
-                                locpnum[j] = l+1;
-                                break;
-                              }
-                          }
-
-                        if (locpnum[j] == 0)
-                          {
-                            PointIndex pi = seg[j];
-                            meshing.AddPoint (mesh.Point(pi), pi);
-
-                            gis.Append (gi);
                     locpnum[j] = gis.Size();
+                    uv_tree.Insert(uv, locpnum[j]);
                 }
             }
 
-                    meshing.AddBoundaryElement (locpnum[0], locpnum[1], gi0, gi1);
+            meshing.AddBoundaryElement (locpnum[0], locpnum[1], gi[0], gi[1]);
         }
       }
-          }
-
-
 
 
     // Philippose - 15/01/2009
@@ -936,171 +406,25 @@ namespace netgen
         res = MESHING2_GIVEUP;
       }
 
-        projecttype = PARAMETERSPACE;
     static Timer t1("rest of loop"); RegionTimer reg1(t1);
       
-        if (res != MESHING2_OK)
+    bool meshing_failed = res != MESHING2_OK;
-          {
+    if(meshing_failed && delete_on_failure)
-            if (notrys == 1)
     {
         for (SurfaceElementIndex sei = noldsurfel; sei < mesh.GetNSE(); sei++)
             mesh.Delete(sei);
 
         mesh.Compress();
-
-                (*testout) << "retry Surface " << k << endl;
-
-                k--;
-                // projecttype*=-1;
-                projecttype = PLANESPACE;
-                notrys++;
-                continue;
     }
-            else
-              {
-                geom.facemeshstatus[k-1] = -1;
-                PrintError ("Problem in Surface mesh generation");
-                surfmesherror++;
-                //	      throw NgException ("Problem in Surface mesh generation");
-              }
-          }
-        else
-          {
-            geom.facemeshstatus[k-1] = 1;
-          }
-
-        notrys = 1;
 
     for (SurfaceElementIndex sei = oldnf; sei < mesh.GetNSE(); sei++)
       mesh[sei].SetIndex (k);
 
     auto n_illegal_trigs = mesh.FindIllegalTrigs();
     PrintMessage (3, n_illegal_trigs, " illegal triangles");
+    return meshing_failed;
   }
 
-    //      ofstream problemfile("occmesh.rep");
-
-    //      problemfile << "SURFACEMESHING" << endl << endl;
-
-    if (surfmesherror)
-      {
-        cout << "WARNING! NOT ALL FACES HAVE BEEN MESHED" << endl;
-        cout << "SURFACE MESHING ERROR OCCURRED IN " << surfmesherror << " FACES:" << endl;
-        for (int i = 1; i <= geom.fmap.Extent(); i++)
-          if (geom.facemeshstatus[i-1] == -1)
-            {
-              cout << "Face " << i << endl;
-              //               problemfile << "problem with face " << i << endl;
-              //               problemfile << "vertices: " << endl;
-              TopExp_Explorer exp0,exp1,exp2;
-              for ( exp0.Init(TopoDS::Face (geom.fmap(i)), TopAbs_WIRE); exp0.More(); exp0.Next() )
-                {
-                  TopoDS_Wire wire = TopoDS::Wire(exp0.Current());
-                  for ( exp1.Init(wire,TopAbs_EDGE); exp1.More(); exp1.Next() )
-                    {
-                      TopoDS_Edge edge = TopoDS::Edge(exp1.Current());
-                      for ( exp2.Init(edge,TopAbs_VERTEX); exp2.More(); exp2.Next() )
-                        {
-                          TopoDS_Vertex vertex = TopoDS::Vertex(exp2.Current());
-                          gp_Pnt point = BRep_Tool::Pnt(vertex);
-                          //                        problemfile << point.X() << " " << point.Y() << " " << point.Z() << endl;
-                        }
-                    }
-                }
-              //               problemfile << endl;
-
-            }
-        cout << endl << endl;
-        cout << "for more information open IGES/STEP Topology Explorer" << endl;
-        //            problemfile.close();
-        throw NgException ("Problem in Surface mesh generation");
-      }
-    else
-      {
-        //         problemfile << "OK" << endl << endl;
-        //         problemfile.close();
-      }
-    
-    for (int i = 0; i < mesh.GetNFD(); i++)
-      mesh.SetBCName (i, mesh.GetFaceDescriptor(i+1).GetBCName());
-    multithread.task = savetask;
-  }
-
-  void OCCOptimizeSurface(OCCGeometry & geom, Mesh & mesh,
-                          const MeshingParameters & mparam)
-  {
-    const char * savetask = multithread.task;
-    multithread.task = "Optimizing surface";
-
-    static Timer timer_opt2d("Optimization 2D");
-    timer_opt2d.Start();
-
-    for (int k = 1; k <= mesh.GetNFD(); k++)
-      {
-        //      if (k != 42) continue;
-        //      if (k != 36) continue;
-
-        //      (*testout) << "optimize face " << k << endl;
-        multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL);
-
-        FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
-
-        PrintMessage (1, "Optimize Surface ", k);
-        for (int i = 1; i <= mparam.optsteps2d; i++)
-          {
-            //	  (*testout) << "optstep " << i << endl;
-            if (multithread.terminate) return;
-
-            {
-              MeshOptimize2d meshopt(mesh);
-              meshopt.SetFaceIndex (k);
-              meshopt.SetImproveEdges (0);
-              meshopt.SetMetricWeight (mparam.elsizeweight);
-              meshopt.SetWriteStatus (0);
-              meshopt.EdgeSwapping (i > mparam.optsteps2d/2);
-            }
-
-            if (multithread.terminate) return;
-            {
-              MeshOptimize2d meshopt(mesh);
-              meshopt.SetFaceIndex (k);
-              meshopt.SetImproveEdges (0);
-              meshopt.SetMetricWeight (mparam.elsizeweight);
-              meshopt.SetWriteStatus (0);
-              meshopt.ImproveMesh (mparam);
-            }
-
-            {
-              MeshOptimize2d meshopt(mesh);
-              meshopt.SetFaceIndex (k);
-              meshopt.SetImproveEdges (0);
-              meshopt.SetMetricWeight (mparam.elsizeweight);
-              meshopt.SetWriteStatus (0);
-              meshopt.CombineImprove ();
-            }
-
-            if (multithread.terminate) return;
-            {
-              MeshOptimize2d meshopt(mesh);
-              meshopt.SetFaceIndex (k);
-              meshopt.SetImproveEdges (0);
-              meshopt.SetMetricWeight (mparam.elsizeweight);
-              meshopt.SetWriteStatus (0);
-              meshopt.ImproveMesh (mparam);
-            }
-          }
-
-      }
-
-
-    mesh.CalcSurfacesOfNode();
-    mesh.Compress();
-    timer_opt2d.Stop();
-
-    multithread.task = savetask;
-  }
-
-
 
   void OCCSetLocalMeshSize(const OCCGeometry & geom, Mesh & mesh,
                            const MeshingParameters & mparam, const OCCParameters& occparam)

libsrc/occ/occgeom.cpp

@@ -1,52 +1,69 @@
 
 #ifdef OCCGEOMETRY
 
-#include <mystdlib.h>
-#include <occgeom.hpp>
-#include <core/register_archive.hpp>
 #include <cstdio>
-#include "ShapeAnalysis_ShapeTolerance.hxx"
+#include <set>
-#include "ShapeAnalysis_ShapeContents.hxx"
-#include "ShapeAnalysis_CheckSmallFace.hxx"
-#include "ShapeAnalysis_DataMapOfShapeListOfReal.hxx"
-#include "ShapeAnalysis_Surface.hxx"
 
-#include "BRepCheck_Analyzer.hxx"
+#include <mystdlib.h>
-#include "BRepLib.hxx"
+#include <core/register_archive.hpp>
-#include "ShapeBuild_ReShape.hxx"
-#include "ShapeFix.hxx"
-#include "ShapeFix_FixSmallFace.hxx"
-#include "Partition_Spliter.hxx"
-#include "BRepAlgoAPI_Fuse.hxx"
-#include "Interface_InterfaceModel.hxx"
 
-#include "XSControl_WorkSession.hxx"
+#include "occ_vertex.hpp"
-#include "XSControl_TransferReader.hxx"
+#include "occ_edge.hpp"
-#include "StepRepr_RepresentationItem.hxx"
+#include "occ_face.hpp"
-#include "StepBasic_ProductDefinitionRelationship.hxx"
+#include "occ_solid.hpp"
-#include "Transfer_TransientProcess.hxx"
+#include "occgeom.hpp"
-#include "TransferBRep.hxx"
-#ifndef _Standard_Version_HeaderFile
-#include <Standard_Version.hxx>
-#endif
 
+#include <BOPAlgo_Builder.hxx>
+#include <BRepBndLib.hxx>
+#include <BRepBuilderAPI_MakeSolid.hxx>
+#include <BRepBuilderAPI_MakeVertex.hxx>
+#include <BRepCheck_Analyzer.hxx>
+#include <BRepExtrema_DistShapeShape.hxx>
+#include <BRepGProp.hxx>
+#include <BRepLib.hxx>
+#include <BRepMesh_IncrementalMesh.hxx>
+#include <BRepOffsetAPI_Sewing.hxx>
+#include <BRepTools.hxx>
+#include <IGESCAFControl_Reader.hxx>
+#include <IGESControl_Writer.hxx>
+#include <Interface_InterfaceModel.hxx>
+#include <Interface_Static.hxx>
+#include <Partition_Spliter.hxx>
+#include <STEPCAFControl_Writer.hxx>
 #include <STEPConstruct.hxx>
-#include <Transfer_FinderProcess.hxx>
+#include <STEPControl_Writer.hxx>
-#include <TDataStd_Name.hxx>
+#include <ShapeAnalysis_CheckSmallFace.hxx>
-#include <XCAFPrs.hxx>
+#include <ShapeAnalysis_DataMapOfShapeListOfReal.hxx>
+#include <ShapeAnalysis_ShapeContents.hxx>
+#include <ShapeBuild_ReShape.hxx>
+#include <ShapeFix_Face.hxx>
+#include <ShapeFix_FixSmallFace.hxx>
+#include <ShapeFix_Shape.hxx>
+#include <ShapeFix_Wire.hxx>
+#include <ShapeFix_Wireframe.hxx>
+#include <StepBasic_ProductDefinitionRelationship.hxx>
+#include <StepRepr_RepresentationItem.hxx>
+#include <StlAPI_Writer.hxx>
 #include <TopoDS_Shape.hxx>
+#include <TransferBRep.hxx>
+#include <Transfer_FinderProcess.hxx>
+#include <Transfer_TransientProcess.hxx>
+#include <XCAFApp_Application.hxx>
+#include <XCAFDoc_ColorTool.hxx>
+#include <XCAFDoc_DocumentTool.hxx>
+#include <XCAFDoc_ShapeTool.hxx>
+#include <XCAFPrs.hxx>
 #include <XCAFPrs_Style.hxx>
-#include <TopTools_ShapeMapHasher.hxx>
+#include <XSControl_TransferReader.hxx>
-#include <NCollection_IndexedDataMap.hxx>
+#include <XSControl_WorkSession.hxx>
-#include <StepShape_TopologicalRepresentationItem.hxx>
 
 #if OCC_VERSION_HEX < 0x070000
 // pass
 #elif OCC_VERSION_HEX < 0x070200
-   #include "StlTransfer.hxx"
+   #include <StlTransfer.hxx>
-   #include "TopoDS_Iterator.hxx"
+   #include <TopoDS_Iterator.hxx>
 #else
-   #include "TopoDS_Iterator.hxx"
+   #include <TopoDS_Iterator.hxx>
 #endif
 
 namespace netgen
@@ -57,6 +74,67 @@ namespace netgen
   std::map<Handle(TopoDS_TShape), ShapeProperties> OCCGeometry::global_shape_properties;
   std::map<Handle(TopoDS_TShape), std::vector<OCCIdentification>> OCCGeometry::identifications;
   
+  TopoDS_Shape ListOfShapes::Max(gp_Vec dir)
+  {
+    double maxval = -1e99;
+    TopoDS_Shape maxshape;
+    for (auto shape : *this)
+      {
+        GProp_GProps props;
+        gp_Pnt center;
+
+        switch (shape.ShapeType())
+          {
+          case TopAbs_VERTEX:
+            center = BRep_Tool::Pnt (TopoDS::Vertex(shape)); break;
+          case TopAbs_FACE:
+            BRepGProp::SurfaceProperties (shape, props);
+            center = props.CentreOfMass();
+            break;
+          default:
+            BRepGProp::LinearProperties(shape, props);
+            center = props.CentreOfMass();
+          }
+
+        double val = center.X()*dir.X() + center.Y()*dir.Y() + center.Z() * dir.Z();
+        if (val > maxval)
+          {
+            maxval = val;
+            maxshape = shape;
+          }
+      }
+    return maxshape;
+  }
+  TopoDS_Shape ListOfShapes::Nearest(gp_Pnt pnt)
+  {
+    double mindist = 1e99;
+    TopoDS_Shape nearestshape;
+    auto vertex = BRepBuilderAPI_MakeVertex (pnt).Vertex();
+    
+    for (auto shape : *this)
+      {
+        double dist = BRepExtrema_DistShapeShape(shape, vertex).Value();
+        if (dist < mindist)
+          {
+            nearestshape = shape;
+            mindist = dist;
+          }
+      }
+    return nearestshape;
+  }
+
+  ListOfShapes ListOfShapes::SubShapes(TopAbs_ShapeEnum type) const
+  {
+    std::set<TopoDS_Shape, ShapeLess> unique_shapes;
+    for(const auto& shape : *this)
+      for(TopExp_Explorer e(shape, type); e.More(); e.Next())
+        unique_shapes.insert(e.Current());
+    ListOfShapes sub;
+    for(const auto& shape : unique_shapes)
+      sub.push_back(shape);
+    return sub;
+  }
+
   OCCGeometry::OCCGeometry(const TopoDS_Shape& _shape, int aoccdim, bool copy)
   {
     if(copy)
@@ -64,14 +142,14 @@ namespace netgen
         auto filename = GetTempFilename();
         step_utils::WriteSTEP(_shape, filename.c_str());
         LoadOCCInto(this, filename.c_str());
-        occdim = aoccdim;
+        dimension = aoccdim;
         std::remove(filename.c_str());
       }
     else
       {
         shape = _shape;
         changed = 1;
-        occdim = aoccdim;
+        dimension = aoccdim;
         BuildFMap();
         CalcBoundingBox();
         PrintContents (this);
@@ -117,23 +195,23 @@ namespace netgen
     OCCSetLocalMeshSize(*this, mesh, mparam, occparam);
   }
 
-  void OCCGeometry :: FindEdges(Mesh& mesh,
+  bool OCCGeometry :: MeshFace(Mesh& mesh,
-                                const MeshingParameters& mparam) const
+                                  const MeshingParameters& mparam, int nr, FlatArray<int, PointIndex> glob2loc) const
   {
-    OCCFindEdges(*this, mesh, mparam);
+      bool failed = OCCMeshFace(*this, mesh, glob2loc, mparam, nr, PARAMETERSPACE, true);
-  }
+      if(failed)
+          failed = OCCMeshFace(*this, mesh, glob2loc, mparam, nr, PLANESPACE, false);
 
-  void OCCGeometry :: MeshSurface(Mesh& mesh,
+      if(failed)
-                                  const MeshingParameters& mparam) const
       {
-    OCCMeshSurface(*this, mesh, mparam);
+          facemeshstatus[nr] = -1;
+          PrintError ("Problem in Surface mesh generation");
       }
-
+      else
-  void OCCGeometry :: FinalizeMesh(Mesh& mesh) const
       {
-    for (int i = 0; i < mesh.GetNDomains(); i++)
+          facemeshstatus[nr] = 1;
-      if (auto name = sprops[i]->name)
+      }
-        mesh.SetMaterial (i+1, *name);
+      return failed;
   }
 
    void OCCGeometry :: PrintNrShapes ()
@@ -1029,31 +1107,125 @@ namespace netgen
 
       fsingular = esingular = vsingular = false;
 
+      NetgenGeometry::Clear();
+      edge_map.clear();
+      vertex_map.clear();
+      face_map.clear();
+      solid_map.clear();
 
-      sprops.SetSize(somap.Extent());
+      // Add shapes
-      for (TopExp_Explorer e(shape, TopAbs_SOLID); e.More(); e.Next())
+      for(auto i1 : Range(1, vmap.Extent()+1))
       {
-          auto s = e.Current();
+          auto v = vmap(i1);
-          sprops[somap.FindIndex(s)-1] = &global_shape_properties[s.TShape()];
+          auto tshape = v.TShape();
+          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;
+
+          if(global_shape_properties.count(tshape)>0)
+              occ_vertex->properties = global_shape_properties[tshape];
+          vertices.Append(std::move(occ_vertex));
       }
 
-      fprops.SetSize(fmap.Extent());
+      for(auto i1 : Range(1, emap.Extent()+1))
-      for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next())
       {
-          auto s = e.Current();
+          auto e = emap(i1);
-          fprops[fmap.FindIndex(s)-1] = &global_shape_properties[s.TShape()];
+          auto tshape = e.TShape();
+          auto edge = TopoDS::Edge(e);
+          if(edge_map.count(tshape)!=0)
+              continue;
+          edge_map[tshape] = 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()]] );
+          occ_edge->properties = global_shape_properties[tshape];
+          edges.Append(std::move(occ_edge));
       }
 
-      eprops.SetSize(emap.Extent());
+      for(auto i1 : Range(1, fmap.Extent()+1))
-      /*
-      for (TopExp_Explorer e(shape, TopAbs_EDGE); e.More(); e.Next())
       {
-          auto s = e.Current();
+          auto f = fmap(i1);
-          eprops[emap.FindIndex(s)-1] = &global_shape_properties[s.TShape()];
+          auto tshape = f.TShape();
+          if(face_map.count(tshape)==0)
+          {
+
+              auto k = faces.Size();
+              face_map[tshape] = k;
+              auto occ_face = make_unique<OCCFace>(f);
+
+              for(auto e : GetEdges(f))
+                  occ_face->edges.Append( edges[edge_map[e.TShape()]].get() );
+
+              if(global_shape_properties.count(tshape)>0)
+                  occ_face->properties = global_shape_properties[tshape];
+              faces.Append(std::move(occ_face));
+
+              if(dimension==2)
+                  for(auto e : GetEdges(f))
+                  {
+                      auto & edge = *edges[edge_map[e.TShape()]];
+                      if(e.Orientation() == TopAbs_REVERSED)
+                          edge.domout = k;
+                      else
+                          edge.domin = k;
                   }
-      */
+          }
-      for (auto [nr,s] : Enumerate(emap))
+      }
-        eprops[nr-1] = &global_shape_properties[s.TShape()];
+
+
+      for(auto i1 : Range(1, somap.Extent()+1))
+      {
+          auto s = somap(i1);
+          auto tshape = s.TShape();
+          int k;
+          if(solid_map.count(tshape)==0)
+          {
+              k = solids.Size();
+              solid_map[tshape] = k;
+              auto occ_solid = make_unique<OCCSolid>(s);
+              if(global_shape_properties.count(tshape)>0)
+                  occ_solid->properties = global_shape_properties[tshape];
+              solids.Append(std::move(occ_solid));
+          }
+
+          for(auto f : GetFaces(s))
+          {
+              auto face_nr = face_map[f.TShape()];
+              auto & face = faces[face_nr];
+              if(face->domin==-1)
+                  face->domin = k;
+              else
+                  face->domout = k;
+          }
+      }
+
+      // Add identifications
+      auto add_identifications = [&](auto & shapes, auto & shape_map)
+      {
+          for(auto &[tshape, nr] : shape_map)
+            if(identifications.count(tshape))
+                for(auto & ident : identifications[tshape])
+                {
+                    if(shape_map.count(ident.from)==0 || shape_map.count(ident.to)==0)
+                        continue;
+
+                    ShapeIdentification si{
+                        shapes[shape_map[ident.from]].get(),
+                        shapes[shape_map[ident.to]].get(),
+                        ident.trafo,
+                        ident.type,
+                        ident.name
+                    };
+                    shapes[nr]->identifications.Append(si);
+                }
+      };
+      add_identifications( vertices, vertex_map );
+      add_identifications( edges, edge_map );
+      add_identifications( faces, face_map );
+
+      bounding_box = ::netgen::GetBoundingBox( shape );
+      ProcessIdentifications();
    }
 
 
@@ -1156,286 +1328,11 @@ namespace netgen
 
    void OCCGeometry :: CalcBoundingBox ()
    {
-      Bnd_Box bb;
+      boundingbox = ::netgen::GetBoundingBox(shape);
-#if OCC_VERSION_HEX < 0x070000
+      (*testout) << "Bounding Box = [" << boundingbox.PMin() << " - " << boundingbox.PMax() << "]" << endl;
-      BRepBndLib::Add (shape, bb);
-#else
-      BRepBndLib::Add ((const TopoDS_Shape) shape, bb,(Standard_Boolean)true);
-#endif
-
-      double x1,y1,z1,x2,y2,z2;
-      bb.Get (x1,y1,z1,x2,y2,z2);
-      Point<3> p1 = Point<3> (x1,y1,z1);
-      Point<3> p2 = Point<3> (x2,y2,z2);
-
-      (*testout) << "Bounding Box = [" << p1 << " - " << p2 << "]" << endl;
-      boundingbox = Box<3> (p1,p2);
       SetCenter();
    }
 
-   PointGeomInfo OCCGeometry :: ProjectPoint(int surfi, Point<3> & p) const
-   {
-      static int cnt = 0;
-      if (++cnt % 1000 == 0) cout << "Project cnt = " << cnt << endl;
-
-      gp_Pnt pnt(p(0), p(1), p(2));
-
-      double u,v;
-      Handle( Geom_Surface ) thesurf = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
-      Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( thesurf );
-      gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfi)) ) );
-      suval.Coord( u, v);
-      pnt = thesurf->Value( u, v );
-
-      PointGeomInfo gi;
-      gi.trignum = surfi;
-      gi.u = u;
-      gi.v = v;
-      p = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
-      return gi;
-   }
-
-  bool OCCGeometry :: ProjectPointGI(int surfind, Point<3>& p, PointGeomInfo& gi) const
-  {
-    double u = gi.u;
-    double v = gi.v;
-
-    Point<3> hp = p;
-    if (FastProject (surfind, hp, u, v))
-      {
-	p = hp;
-	return 1;
-      }
-    ProjectPoint (surfind, p);
-    return CalcPointGeomInfo (surfind, gi, p);
-  }
-
-  void OCCGeometry :: ProjectPointEdge(int surfind, INDEX surfind2,
-                                       Point<3> & p, EdgePointGeomInfo* gi) const
-  {
-    TopExp_Explorer exp0, exp1;
-    bool done = false;
-    Handle(Geom_Curve) c;
-
-    for (exp0.Init(fmap(surfind), TopAbs_EDGE); !done && exp0.More(); exp0.Next())
-      for (exp1.Init(fmap(surfind2), TopAbs_EDGE); !done && exp1.More(); exp1.Next())
-	{
-	  if (TopoDS::Edge(exp0.Current()).IsSame(TopoDS::Edge(exp1.Current())))
-	    {
-	      done = true;
-	      double s0, s1;
-	      c = BRep_Tool::Curve(TopoDS::Edge(exp0.Current()), s0, s1);
-	    }
-	}
-
-    gp_Pnt pnt(p(0), p(1), p(2));
-    GeomAPI_ProjectPointOnCurve proj(pnt, c);
-    pnt = proj.NearestPoint();
-    p(0) = pnt.X();
-    p(1) = pnt.Y();
-    p(2) = pnt.Z();
-
-  }
-
-   bool OCCGeometry :: FastProject (int surfi, Point<3> & ap, double& u, double& v) const
-   {
-      gp_Pnt p(ap(0), ap(1), ap(2));
-
-      Handle(Geom_Surface) surface = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
-
-      gp_Pnt x = surface->Value (u,v);
-
-      if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true;
-
-      gp_Vec du, dv;
-
-      surface->D1(u,v,x,du,dv);
-
-      int count = 0;
-
-      gp_Pnt xold;
-      gp_Vec n;
-      double det, lambda, mu;
-
-      do {
-         count++;
-
-         n = du^dv;
-
-         det = Det3 (n.X(), du.X(), dv.X(),
-            n.Y(), du.Y(), dv.Y(),
-            n.Z(), du.Z(), dv.Z());
-
-         if (det < 1e-15) return false;
-
-         lambda = Det3 (n.X(), p.X()-x.X(), dv.X(),
-            n.Y(), p.Y()-x.Y(), dv.Y(),
-            n.Z(), p.Z()-x.Z(), dv.Z())/det;
-
-         mu     = Det3 (n.X(), du.X(), p.X()-x.X(),
-            n.Y(), du.Y(), p.Y()-x.Y(),
-            n.Z(), du.Z(), p.Z()-x.Z())/det;
-
-         u += lambda;
-         v += mu;
-
-         xold = x;
-         surface->D1(u,v,x,du,dv);
-
-      } while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50);
-
-      //    (*testout) << "FastProject count: " << count << endl;
-
-      if (count == 50) return false;
-
-      ap = Point<3> (x.X(), x.Y(), x.Z());
-
-      return true;
-   }
-
-  Vec<3> OCCGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* geominfo) const
-  {
-    if(geominfo)
-      {
-        gp_Pnt pnt;
-        gp_Vec du, dv;
-
-        Handle(Geom_Surface) occface;
-        occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind)));
-
-        occface->D1(geominfo->u,geominfo->v,pnt,du,dv);
-
-        auto n = Cross (Vec<3>(du.X(), du.Y(), du.Z()),
-                        Vec<3>(dv.X(), dv.Y(), dv.Z()));
-        n.Normalize();
-
-        if (fmap(surfind).Orientation() == TopAbs_REVERSED) n *= -1;
-        return n;
-      }
-    Standard_Real u,v;
-
-    gp_Pnt pnt(p(0), p(1), p(2));
-
-    Handle(Geom_Surface) occface;
-    occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind)));
-
-    /*
-    GeomAPI_ProjectPointOnSurf proj(pnt, occface);
-
-    if (proj.NbPoints() < 1)
-      {
-	cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
-	     << endl;
-	cout << p << endl;
-	return;
-      }
- 
-    proj.LowerDistanceParameters (u, v);
-    */
-    
-    Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
-    gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfind)) ) );
-    suval.Coord( u, v);
-    pnt = occface->Value( u, v );
-
-    gp_Vec du, dv;
-    occface->D1(u,v,pnt,du,dv);
-
-    /*
-      if (!occface->IsCNu (1) || !occface->IsCNv (1))
-      (*testout) << "SurfOpt: Differentiation FAIL" << endl;
-    */
-
-    auto n = Cross (Vec3d(du.X(), du.Y(), du.Z()),
-	       Vec3d(dv.X(), dv.Y(), dv.Z()));
-    n.Normalize();
-
-    if (fmap(surfind).Orientation() == TopAbs_REVERSED) n *= -1;
-    return n;
-  }
-
-  bool OCCGeometry :: CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p) const
-  {
-    Standard_Real u,v;
-
-    gp_Pnt pnt(p(0), p(1), p(2));
-
-    Handle(Geom_Surface) occface;
-    occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind)));
-
-    /*
-    GeomAPI_ProjectPointOnSurf proj(pnt, occface);
-
-    if (proj.NbPoints() < 1)
-      {
-	cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
-	     << endl;
-	cout << p << endl;
-	return 0;
-      }
- 
-    proj.LowerDistanceParameters (u, v);  
-    */
-
-    Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
-    gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfind)) ) );
-    suval.Coord( u, v);
-    //pnt = occface->Value( u, v );
-    
-
-    gi.u = u;
-    gi.v = v;
-    return true;
-  }
-
-  void OCCGeometry :: PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
-                                   int surfi, 
-                                   const PointGeomInfo & gi1, 
-                                   const PointGeomInfo & gi2,
-                                   Point<3> & newp, PointGeomInfo & newgi) const
-  {
-    Point<3> hnewp;
-    hnewp = p1+secpoint*(p2-p1);
-
-    if (surfi > 0)
-      {
-	double u = gi1.u+secpoint*(gi2.u-gi1.u);
-	double v = gi1.v+secpoint*(gi2.v-gi1.v);
-
-        auto savept = hnewp;
-	if (!FastProject(surfi, hnewp, u, v) || Dist(hnewp, savept) > Dist(p1,p2))
-	  {
-            //  cout << "Fast projection to surface fails! Using OCC projection" << endl;
-            hnewp = savept;
-	    ProjectPoint(surfi, hnewp);
-	  }
-	newgi.trignum = 1;
-        newgi.u = u;
-        newgi.v = v;
-      }
-    newp = hnewp;
-  }
-
-
-  void OCCGeometry :: PointBetweenEdge(const Point<3> & p1,
-                                       const Point<3> & p2, double secpoint,
-                                       int surfi1, int surfi2, 
-                                       const EdgePointGeomInfo & ap1, 
-                                       const EdgePointGeomInfo & ap2,
-                                       Point<3> & newp, EdgePointGeomInfo & newgi) const
-  {
-    double s0, s1;
-
-    newp = p1+secpoint*(p2-p1);
-    if(ap1.edgenr > emap.Size() || ap1.edgenr == 0)
-      return;
-    gp_Pnt pnt(newp(0), newp(1), newp(2));
-    GeomAPI_ProjectPointOnCurve proj(pnt, BRep_Tool::Curve(TopoDS::Edge(emap(ap1.edgenr)), s0, s1));
-    pnt = proj.NearestPoint();
-    newp = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
-    newgi = ap1;
-  };
-
 
 //    void OCCGeometry :: WriteOCC_STL(char * filename)
 //    {
@@ -1699,26 +1596,18 @@ namespace netgen
         auto occ_hash = key.HashCode(1<<31UL);
         return std::hash<decltype(occ_hash)>()(occ_hash);
     };
-    std::unordered_map<TopoDS_Shape, int, decltype(my_hash)> shape_map(10, my_hash);
-    Array<TopoDS_Shape> shape_list;
-
     std::map<Handle(TopoDS_TShape), int> tshape_map;
     Array<Handle(TopoDS_TShape)> tshape_list;
 
-    ar & occdim;
+    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)
-            {
-              shape_map[ds] = shape_list.Size();
-              shape_list.Append(ds);
-            }
           if(tshape_map.count(ts)==0)
             {
-              tshape_map[ts] = shape_list.Size();
+              tshape_map[ts] = tshape_list.Size();
               tshape_list.Append(ts);
             }
         }
@@ -1741,12 +1630,18 @@ namespace netgen
             for(auto i : Range(n_idents))
               {
                 auto & id = idents[i];
-                int shape_id;
+                int id_from, id_to;
                 if(ar.Output())
-                    shape_id = shape_map[id.other];
+                {
-                ar & shape_id & id.trafo & id.inverse & id.name;
+                    id_from = tshape_map[id.from];
+                    id_to = tshape_map[id.to];
+                }
+                ar & id_from & id_to & id.trafo & id.name;
                 if(ar.Input())
-                    id.other = shape_list[shape_id];
+                {
+                    id.from = tshape_list[id_from];
+                    id.to = tshape_list[id_to];
+                }
               }
           }
       }
@@ -2053,6 +1948,114 @@ namespace netgen
       return false;
    }
 
+  Point<3> GetCenter(const TopoDS_Shape & shape)
+  {
+      GProp_GProps props;
+      BRepGProp::LinearProperties(shape, props);
+      return occ2ng( props.CentreOfMass() );
+  }
+
+  bool IsMappedShape(const Transformation<3> & trafo, const TopoDS_Shape & me, const TopoDS_Shape & you)
+  {
+      if(me.ShapeType() != you.ShapeType()) return false;
+
+      Bnd_Box bbox;
+      BRepBndLib::Add(me, bbox);
+      BRepBndLib::Add(you, bbox);
+      BoxTree<3> tree( occ2ng(bbox.CornerMin()), occ2ng(bbox.CornerMax()) );
+
+      Point<3> c_me = GetCenter(me);
+      Point<3> c_you = GetCenter(you);
+      if(tree.GetTolerance() < Dist(trafo(c_me), c_you))
+          return false;
+
+      std::map<T_Shape, T_Shape> vmap;
+
+      auto verts_me = GetVertices(me);
+      auto verts_you = GetVertices(you);
+
+      if(verts_me.size() != verts_you.size())
+          return false;
+
+      for (auto i : Range(verts_me.size()))
+      {
+          auto s = verts_me[i].TShape();
+          if(vmap.count(s)>0)
+              continue;
+          auto p = trafo(occ2ng(s));
+          tree.Insert( p, i );
+          vmap[s] = nullptr;
+      }
+          
+      for (auto vert : verts_you)
+      {
+          auto s = vert.TShape();
+          auto p = occ2ng(s);
+          bool vert_mapped = false;
+          tree.GetFirstIntersecting( p, p, [&](size_t i ) {
+                  vmap[verts_me[i].TShape()] = s;
+                  vert_mapped = true;
+                  return true;
+          });
+          if(!vert_mapped)
+              return false;
+      }
+      return true;
+  }
+
+  void Identify(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional<gp_Trsf> opt_trafo) 
+  {
+    gp_Trsf trafo;
+    if(opt_trafo)
+    {
+        trafo = *opt_trafo;
+    }
+    else
+    {
+        auto v = GetCenter(you) - GetCenter(me);
+        trafo.SetTranslation(gp_Vec(v[0], v[1], v[2]));
+    }
+
+    ListOfShapes list_me, list_you;
+    list_me.push_back(me);
+    list_you.push_back(you);
+    Identify(list_me, list_you, name, type, trafo);
+  }
+
+  void Identify(const ListOfShapes & me, const ListOfShapes & you, string name, Identifications::ID_TYPE type, gp_Trsf occ_trafo) 
+  {
+    Transformation<3> trafo = occ2ng(occ_trafo);
+
+    ListOfShapes id_me;
+    ListOfShapes id_you;
+
+    if(auto faces_me = me.Faces(); faces_me.size()>0)
+    {
+        id_me = faces_me;
+        id_you = you.Faces();
+    }
+    else if(auto edges_me = me.Edges(); edges_me.size()>0)
+    {
+        id_me = edges_me;
+        id_you = you.Edges();
+    }
+    else
+    {
+        id_me = me.Vertices();
+        id_you = you.Vertices();
+    }
+
+    for(auto shape_me : id_me)
+        for(auto shape_you : id_you)
+        {
+            if(!IsMappedShape(trafo, shape_me, shape_you))
+                continue;
+
+            OCCGeometry::identifications[shape_me.TShape()].push_back
+                (OCCIdentification { shape_me.TShape(), shape_you.TShape(), trafo, name, type });
+        }
+  }
+
   void OCCParameters :: Print(ostream & ost) const
    {
       ost << "OCC Parameters:" << endl
@@ -2063,6 +2066,18 @@ namespace netgen
   DLL_HEADER extern OCCParameters occparam;
   OCCParameters occparam;
 
+
+
+
+
+
+
+
+
+
+
+
+
   // int OCCGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam)
   //  {
   //    return OCCGenerateMesh (*this, mesh, mparam, occparam);
@@ -2208,8 +2223,7 @@ namespace netgen
           for(auto & ident : identifications)
           {
               Array<Handle(StepRepr_RepresentationItem)> items;
-              items.Append(STEPConstruct::FindEntity(finder, ident.other));
+              // items.Append(STEPConstruct::FindEntity(finder, ident.other)); // TODO!
-              items.Append(MakeInt(static_cast<int>(ident.inverse)));
               auto & m = ident.trafo.GetMatrix();
               for(auto i : Range(9))
                   items.Append(MakeReal(m(i)));
@@ -2239,8 +2253,7 @@ namespace netgen
               auto id_item = Handle(StepRepr_CompoundRepresentationItem)::DownCast(idents->ItemElementValue(i));
               OCCIdentification ident;
               ident.name = id_item->Name()->ToCString();
-              ident.other = TransferBRep::ShapeResult(transProc->Find(id_item->ItemElementValue(1)));
+              // ident.other = TransferBRep::ShapeResult(transProc->Find(id_item->ItemElementValue(1))); /TODO!
-              ident.inverse = static_cast<bool>(ReadInt(id_item->ItemElementValue(2)));
 
               auto & m = ident.trafo.GetMatrix();
               for(auto i : Range(9))

libsrc/occ/occgeom.hpp

@@ -9,100 +9,30 @@
 
 #ifdef OCCGEOMETRY
 
+#include <set>
+
 #include <meshing.hpp>
+#include "occ_utils.hpp"
+#include "occmeshsurf.hpp"
 
-#include <Standard_Version.hxx>
+#include <Quantity_ColorRGBA.hxx>
-#include "BRep_Tool.hxx"
+#include <STEPCAFControl_Reader.hxx>
-#include "Geom_Curve.hxx"
-#include "Geom2d_Curve.hxx"
-#include "Geom_Surface.hxx"
-#include "GeomAPI_ProjectPointOnSurf.hxx"
-#include "GeomAPI_ProjectPointOnCurve.hxx"
-#include "BRepTools.hxx"
-#include "TopExp.hxx"
-#include "BRepBuilderAPI_MakeVertex.hxx"
-#include "BRepBuilderAPI_MakeShell.hxx"
-#include "BRepBuilderAPI_MakeSolid.hxx"
-#include "BRepOffsetAPI_Sewing.hxx"
-#include "BRepLProp_SLProps.hxx"
-#include "BRepAdaptor_Surface.hxx"
-#include "Poly_Triangulation.hxx"
-#include "Poly_Array1OfTriangle.hxx"
-#include "TColgp_Array1OfPnt2d.hxx"
-#include "Poly_Triangle.hxx"
-#include "GProp_GProps.hxx"
-#include "BRepGProp.hxx"
-#include "gp_Pnt.hxx"
-#include "TopoDS.hxx"
-#include "TopoDS_Solid.hxx"
-#include "TopExp_Explorer.hxx"
-#include "TopTools_ListIteratorOfListOfShape.hxx"
-#include "TopoDS_Wire.hxx"
-#include "BRepTools_WireExplorer.hxx"
-#include "TopTools_IndexedMapOfShape.hxx"
-#include "BRepLProp_CLProps.hxx"
-#include "BRepAdaptor_Curve.hxx"
-#include "TopoDS_Shape.hxx"
-#include "TopoDS_Face.hxx"
-#include "IGESToBRep_Reader.hxx"
-#include "Interface_Static.hxx"
-#include "GeomAPI_ExtremaCurveCurve.hxx"
-#include "Standard_ErrorHandler.hxx"
-#include "Standard_Failure.hxx"
-#include "ShapeUpgrade_ShellSewing.hxx"
-#include "ShapeFix_Shape.hxx"
-#include "ShapeFix_Wireframe.hxx"
-#include "BRepMesh_IncrementalMesh.hxx"
-#include "BRepBndLib.hxx"
-#include "Bnd_Box.hxx"
-#include "ShapeAnalysis.hxx"
-#include "ShapeBuild_ReShape.hxx"
-#include "BOPAlgo_Builder.hxx"
-
-// Philippose - 29/01/2009
-// OpenCascade XDE Support
-// Include support for OpenCascade XDE Features
-#include "TDocStd_Document.hxx"
-#include "Quantity_Color.hxx"
-#include "XCAFApp_Application.hxx"
-#include "XCAFDoc_ShapeTool.hxx"
-#include "XCAFDoc_Color.hxx"
-#include "XCAFDoc_ColorTool.hxx"
-#include "XCAFDoc_ColorType.hxx"
-#include "XCAFDoc_LayerTool.hxx"
-#include "XCAFDoc_DimTolTool.hxx"
-#include "XCAFDoc_MaterialTool.hxx"
-#include "XCAFDoc_DocumentTool.hxx"
-#include "TDF_Label.hxx"
-#include "TDF_LabelSequence.hxx"
-#include "STEPCAFControl_Reader.hxx"
-#include "STEPCAFControl_Writer.hxx"
-#include "IGESCAFControl_Reader.hxx"
-#include "IGESCAFControl_Writer.hxx"
-
-#include "IGESControl_Reader.hxx"
-#include "STEPControl_Reader.hxx"
-#include "IGESControl_Writer.hxx"
-#include "STEPControl_Writer.hxx"
-
-#include <StepRepr_ValueRepresentationItem.hxx>
-#include <StepRepr_IntegerRepresentationItem.hxx>
-#include <StepRepr_CompoundRepresentationItem.hxx>
 #include <StepBasic_MeasureValueMember.hxx>
-
+#include <StepRepr_CompoundRepresentationItem.hxx>
-#include "StlAPI_Writer.hxx"
+#include <StepRepr_IntegerRepresentationItem.hxx>
-#include "STEPControl_StepModelType.hxx"
+#include <StepRepr_ValueRepresentationItem.hxx>
+#include <TCollection_HAsciiString.hxx>
+#include <TDocStd_Document.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Shape.hxx>
+#include <Transfer_FinderProcess.hxx>
 
 #if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4
 #define OCC_HAVE_HISTORY
 #endif
 
-
-
-
 namespace netgen
 {
-#include "occmeshsurf.hpp"
 
   // extern DLL_HEADER MeshingParameters mparam;
 
@@ -116,28 +46,7 @@ namespace netgen
 #define OCCGEOMETRYVISUALIZATIONFULLCHANGE 1  // Compute transformation matrices and redraw
 #define OCCGEOMETRYVISUALIZATIONHALFCHANGE 2  // Redraw
 
-
+  bool IsMappedShape(const Transformation<3> & trafo, const TopoDS_Shape & me, const TopoDS_Shape & you);
-  inline Point<3> occ2ng (const gp_Pnt & p)
-  {
-    return Point<3> (p.X(), p.Y(), p.Z());
-  }
-
-  inline Point<2> occ2ng (const gp_Pnt2d & p)
-  {
-    return Point<2> (p.X(), p.Y());
-  }
-  
-  inline Vec<3> occ2ng (const gp_Vec & v)
-  {
-    return Vec<3> (v.X(), v.Y(), v.Z());
-  }
-
-  inline gp_Pnt ng2occ (const Point<3> & p)
-  {
-    return gp_Pnt(p(0), p(1), p(2));
-  }
-
-  
 
   class EntityVisualizationCode
   {
@@ -219,108 +128,20 @@ namespace netgen
   };
 
 
-  class ShapeProperties
-  {
-  public:
-    optional<string> name;
-    optional<Vec<4>> col;
-    double maxh = 1e99;
-    double hpref = 0;  // number of hp refinement levels (will be multiplied by factor later)
-    void Merge(const ShapeProperties & prop2)
-    {
-      if (prop2.name) name = prop2.name;
-      if (prop2.col) col = prop2.col;
-      maxh = min2(maxh, prop2.maxh);
-    }
-
-    void DoArchive(Archive& ar)
-    {
-        ar & name & col & maxh & hpref;
-    }
-  };
-
-
-  class OCCIdentification
-  {
-  public:
-    TopoDS_Shape other;
-    Transformation<3> trafo;
-    bool inverse;
-    string name;
-  };
-
-
-  class MyExplorer
-  {
-    class Iterator
-    {
-      TopExp_Explorer exp;
-    public:
-      Iterator (TopoDS_Shape ashape, TopAbs_ShapeEnum atoFind, TopAbs_ShapeEnum atoAvoid)
-        : exp(ashape, atoFind, atoAvoid) { }
-      auto operator*() { return exp.Current(); }
-      Iterator & operator++() { exp.Next(); return *this; }
-      bool operator!= (nullptr_t nu) { return exp.More(); }
-    };
-
-  public:
-    TopoDS_Shape shape;
-    TopAbs_ShapeEnum toFind;
-    TopAbs_ShapeEnum toAvoid;
-    MyExplorer (TopoDS_Shape ashape, TopAbs_ShapeEnum atoFind, TopAbs_ShapeEnum atoAvoid = TopAbs_SHAPE)
-      : shape(ashape), toFind(atoFind), toAvoid(atoAvoid) { ; }
-    Iterator begin() { return Iterator(shape, toFind, toAvoid); }
-    auto end() { return nullptr; }
-  };
-
-  inline auto Explore (TopoDS_Shape shape, TopAbs_ShapeEnum toFind, TopAbs_ShapeEnum toAvoid = TopAbs_SHAPE)
-  {
-    return MyExplorer (shape, toFind, toAvoid);
-  }
-
-  
-  class IndexMapIterator
-  {
-    class Iterator
-    {
-      const TopTools_IndexedMapOfShape & indmap;
-      int i;
-    public:
-      Iterator (const TopTools_IndexedMapOfShape & aindmap, int ai)
-        : indmap(aindmap), i(ai) { ; }
-      auto operator*() { return tuple(i, indmap(i)); }
-      Iterator & operator++() { i++; return *this; }
-      bool operator!= (const Iterator & i2) { return i != i2.i; }
-    };
-
-  public:
-    const TopTools_IndexedMapOfShape & indmap;
-    IndexMapIterator (const TopTools_IndexedMapOfShape & aindmap) : indmap(aindmap) { }
-    Iterator begin() { return Iterator(indmap, 1); }
-    Iterator end() { return Iterator(indmap, indmap.Extent()+1); }
-  };
-  
-  inline auto Enumerate (const TopTools_IndexedMapOfShape & indmap)
-  {
-    return IndexMapIterator(indmap);
-  }
-
-  
   class DLL_HEADER OCCGeometry : public NetgenGeometry
   {
     Point<3> center;
     OCCParameters occparam;
   public:
-    static std::map<Handle(TopoDS_TShape), ShapeProperties> global_shape_properties;
+    static std::map<T_Shape, ShapeProperties> global_shape_properties;
-    static std::map<Handle(TopoDS_TShape), std::vector<OCCIdentification>> identifications;
+    static std::map<T_Shape, std::vector<OCCIdentification>> identifications;
 
     TopoDS_Shape shape;
-    TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap;
+    TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap; // legacy maps
     NgArray<bool> fsingular, esingular, vsingular;
     Box<3> boundingbox;
 
-    // should we use 1-based arrays (JS->MH) ? 
+    std::map<T_Shape, int> edge_map, vertex_map, face_map, solid_map;
-    Array<ShapeProperties*> fprops, eprops, sprops; // pointers to the gobal property map
 
     mutable int changed;
     mutable NgArray<int> facemeshstatus;
@@ -350,8 +171,6 @@ namespace netgen
     bool makesolids;
     bool splitpartitions;
 
-    int occdim = 3; // meshing is always done 3D, changed to 2D later of occdim=2
-     
     OCCGeometry()
     {
       somap.Clear();
@@ -372,36 +191,15 @@ namespace netgen
 
     void Analyse(Mesh& mesh,
                  const MeshingParameters& mparam) const override;
-    void FindEdges(Mesh& mesh,
+    bool MeshFace(Mesh& mesh, const MeshingParameters& mparam,
-                   const MeshingParameters& mparam) const override;
+                     int nr, FlatArray<int, PointIndex> glob2loc) const override;
-    void MeshSurface(Mesh& mesh,
+    // void OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const override {}
-                     const MeshingParameters& mparam) const override;
- 
-    void FinalizeMesh(Mesh& mesh) const override;
  
     void Save (string filename) const override;
     void SaveToMeshFile (ostream & /* ost */) const override;
      
     void DoArchive(Archive& ar) override;
 
-    PointGeomInfo ProjectPoint(int surfind, Point<3> & p) const override;
-    void ProjectPointEdge (int surfind, int surfind2, Point<3> & p,
-                           EdgePointGeomInfo* gi = nullptr) const override;
-    bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override;
-    Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const override;
-    bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const override;
-
-    void PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
-                          int surfi1, int surfi2, 
-                          const EdgePointGeomInfo & ap1, 
-                          const EdgePointGeomInfo & ap2,
-                          Point<3> & newp, EdgePointGeomInfo & newgi) const override;
-    void PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
-                      int surfi, 
-                      const PointGeomInfo & gi1, 
-                      const PointGeomInfo & gi2,
-                      Point<3> & newp, PointGeomInfo & newgi) const override;
-
     void BuildFMap();
 
     auto GetShape() const { return shape; }
@@ -548,9 +346,12 @@ namespace netgen
     //      void WriteOCC_STL(char * filename);
 
   private:
-    bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
+    //bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
   };
 
+  void Identify(const ListOfShapes & me, const ListOfShapes & you, string name, Identifications::ID_TYPE type, gp_Trsf occ_trafo);
+  void Identify(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional<gp_Trsf> opt_trafo);
+   
 
   void PrintContents (OCCGeometry * geom);
 
@@ -564,12 +365,107 @@ namespace netgen
   DLL_HEADER extern void OCCSetLocalMeshSize(const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam,
                                              const OCCParameters& occparam);
 
-  DLL_HEADER extern void OCCMeshSurface (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
+  DLL_HEADER extern bool OCCMeshFace (const OCCGeometry & geom, Mesh & mesh, FlatArray<int, PointIndex> glob2loc,
+                       const MeshingParameters & mparam, int nr, int projecttype, bool delete_on_failure);
 
-  DLL_HEADER extern void OCCOptimizeSurface (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
 
-  DLL_HEADER extern void OCCFindEdges (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam);
+  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<T_Shape, bool> tshape_handled;
+    Transformation<3> trafo_inv;
+    if(trafo)
+        trafo_inv = trafo->CalcInverse();
   
+    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();
+          mod_map[tshape].insert(tshape);
+          tshape_handled[tshape] = 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();
+  
+          for (auto mods : builder.Modified(e.Current()))
+              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();
+  
+          if(tshape_handled[tshape])
+              continue;
+          tshape_handled[tshape] = true;
+  
+          if(OCCGeometry::identifications.count(tshape)==0)
+              continue;
+  
+          auto tshape_mapped = mod_map[tshape];
+  
+          for(auto ident : OCCGeometry::identifications[tshape])
+          {
+              // nothing happened
+              if(mod_map[ident.to].size()==1 && mod_map[ident.from].size() ==1)
+                  continue;
+  
+              auto from = ident.from;
+              auto to = ident.to;
+  
+              for(auto from_mapped : mod_map[from])
+                  for(auto to_mapped : mod_map[to])
+                  {
+                      if(from==from_mapped && to==to_mapped)
+                          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)
+                      {
+                          Transformation<3> trafo_temp;
+                          trafo_temp.Combine(ident.trafo, trafo_inv);
+                          trafo_mapped.Combine(*trafo, trafo_temp);
+                      }
+  
+                      if(!IsMappedShape(trafo_mapped, s_from, s_to))
+                          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;
+                      OCCGeometry::identifications[id_owner].push_back(id_new);
+                  }
+          }
+      }
+  }
+  
+  template <class TBuilder>
+  void PropagateProperties (TBuilder & builder, TopoDS_Shape shape, std::optional<Transformation<3>> trafo = nullopt)
+  {
+    bool have_identifications = false;
+  
+    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 & prop = OCCGeometry::global_shape_properties[tshape];
+          for (auto mods : builder.Modified(e.Current()))
+            OCCGeometry::global_shape_properties[mods.TShape()].Merge(prop);
+          have_identifications |= OCCGeometry::identifications.count(tshape) > 0;
+        }
+    if(have_identifications)
+        PropagateIdentifications(builder, shape, trafo);
+  }
 
   namespace step_utils
   {

libsrc/occ/occmeshsurf.cpp

@@ -2,16 +2,16 @@
 
 #include <mystdlib.h>
 
-#include <occgeom.hpp>
 #include <meshing.hpp>
+#include "occgeom.hpp"
+
 #include <GeomLProp_SLProps.hxx>
 #include <ShapeAnalysis_Surface.hxx>
 
+#include "occmeshsurf.hpp"
 
 namespace netgen
 {
-#include "occmeshsurf.hpp"
-
 
   bool glob_testout(false);
 

libsrc/occ/occmeshsurf.hpp

@@ -6,9 +6,15 @@
 #include "occgeom.hpp"
 #include "mydefs.hpp"
 
+#include <TopoDS_Face.hxx>
+#include <Geom_Surface.hxx>
+#include <ShapeAnalysis.hxx>
+
 #define PARAMETERSPACE -1
 #define PLANESPACE     1
 
+namespace netgen
+{
 class OCCGeometry;
 
 class SingularMatrixException
@@ -144,8 +150,8 @@ protected:
 
 class OCCGeometry;
 
-#endif
+} // namespace netgen
-
+
-
+#endif
 
 #endif

libsrc/occ/occpkg.cpp

@@ -15,6 +15,9 @@
 
 #include "vsocc.hpp"
 
+#include <TopoDS_Edge.hxx>
+#include <IGESControl_Writer.hxx>
+
 // __declspec(dllimport) void AutoColourBcProps(Mesh & mesh, const char *bccolourfile);
 // __declspec(dllimport) void GetFaceColours(Mesh & mesh, NgArray<Vec3d> & face_colours);
 // __declspec(dllimport) bool ColourMatch(Vec3d col1, Vec3d col2, double eps = 2.5e-05);

libsrc/occ/python_occ.cpp

@@ -1,58 +1,25 @@
 #ifdef NG_PYTHON
 #ifdef OCCGEOMETRY
 
-#include <../general/ngpython.hpp>
-#include <core/python_ngcore.hpp>
-#include "../meshing/python_mesh.hpp"
 #include <memory>
 
+#include <general/ngpython.hpp>
+#include <core/python_ngcore.hpp>
+#include <meshing/python_mesh.hpp>
 #include <meshing.hpp>
-#include <occgeom.hpp>
 
-#include <gp_Ax1.hxx>
+#include "occgeom.hpp"
-#include <gp_Ax2.hxx>
-#include <gp_Ax2d.hxx>
-#include <gp_Trsf.hxx>
-#include <BRepPrimAPI_MakeSphere.hxx>
-#include <BRepPrimAPI_MakeCylinder.hxx>
-#include <BRepPrimAPI_MakeRevol.hxx>
-#include <BRepPrimAPI_MakeBox.hxx>
-#include <BRepPrimAPI_MakePrism.hxx>
-#include <BRepOffsetAPI_MakePipe.hxx>
-#include <BRepAlgoAPI_Cut.hxx>
-#include <BRepAlgoAPI_Common.hxx>
-#include <BRepAlgoAPI_Fuse.hxx>
-// #include <XCAFDoc_VisMaterialTool.hxx>
-#include <TDF_Attribute.hxx>
-#include <Standard_GUID.hxx>
-#include <Geom_TrimmedCurve.hxx>
-#include <Geom_Plane.hxx>
-#include <GC_MakeSegment.hxx>
-#include <GC_MakeCircle.hxx>
-#include <GC_MakeArcOfCircle.hxx>
-#include <GC_MakePlane.hxx>
-#include <BRepBuilderAPI_MakeEdge.hxx>
-#include <BRepBuilderAPI_MakeWire.hxx>
-#include <BRepBuilderAPI_Transform.hxx>
-#include <BRepBuilderAPI_MakeFace.hxx>
-#include <BRepFilletAPI_MakeFillet.hxx>
-#include <BRepOffsetAPI_ThruSections.hxx>
-
-#include <BRepGProp.hxx>
-#include <BRepOffsetAPI_MakeThickSolid.hxx>
-#include <BRepLib.hxx>
-
-#include <Geom2d_Curve.hxx>
-#include <Geom2d_Ellipse.hxx>
-#include <Geom2d_TrimmedCurve.hxx>
-#include <GCE2d_MakeSegment.hxx>
-#include <GCE2d_MakeCircle.hxx>
 
+#include <BOPAlgo_Builder.hxx>
+#include <BRepLProp_SLProps.hxx>
 #include <Message.hxx>
-
+#include <Standard_GUID.hxx>
-#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4
+#include <Standard_Version.hxx>
-#define OCC_HAVE_DUMP_JSON
+#include <TDF_Attribute.hxx>
-#endif
+#include <XCAFApp_Application.hxx>
+#include <XCAFDoc_DocumentTool.hxx>
+#include <XCAFDoc_MaterialTool.hxx>
+#include <XCAFDoc_ShapeTool.hxx>
 
 using namespace netgen;
 
@@ -199,6 +166,10 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                             {
                               self.SetFaceMaxH(fnr, meshsize);
                             }, "Set maximum meshsize for face fnr. Face numbers are 0 based.")
+    .def("Draw", [](shared_ptr<OCCGeometry> geo)
+    {
+      ng_geometry = geo;
+    })
     .def("_visualizationData", [] (shared_ptr<OCCGeometry> occ_geo)
          {
            std::vector<float> vertices;
@@ -294,12 +265,10 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                            geo->SetOCCParameters(occparam);
                            auto mesh = make_shared<Mesh>();
                            mesh->SetGeometry(geo);
+                           SetGlobalMesh(mesh);
                            auto result = geo->GenerateMesh(mesh, mp);
                            if(result != 0)
                              throw Exception("Meshing failed!");
-                           if (geo->occdim==2)
-                             mesh->SetDimension(2);
-                           SetGlobalMesh(mesh);
                            ng_geometry = geo;
                            return mesh;
                          }, py::arg("mp") = nullptr,

libsrc/occ/python_occ.hpp

@@ -1,68 +0,0 @@
-class DirectionalInterval
-{
-public:
-  gp_Vec dir;
-  double minval = -1e99;
-  double maxval = 1e99;
-  bool openmin = false, openmax = false;
-  
-  DirectionalInterval (gp_Vec adir) : dir(adir) { ; }
-  DirectionalInterval (const DirectionalInterval & i2)
-    : dir(i2.dir), minval(i2.minval), maxval(i2.maxval) { ; }
-  
-  DirectionalInterval operator< (double val) const
-  {
-    DirectionalInterval i2 = *this;
-    i2.maxval = val;
-    return i2;
-  }
-
-  DirectionalInterval operator> (double val) const
-  {
-    DirectionalInterval i2 = *this;
-    i2.minval = val;
-    return i2;
-  }
-
-
-  DirectionalInterval Intersect (const DirectionalInterval & i2)
-  {
-    DirectionalInterval res = *this;
-    res.minval = max(res.minval, i2.minval);
-    res.maxval = min(res.maxval, i2.maxval);
-    return res;
-  }
-  
-  bool Contains (gp_Pnt p, double eps = 1e-8)
-  {
-    // cout << "Contains point " << p.X() << "," << p.Y() << "," << p.Z() << " ? " << endl;
-    double val = dir.X()*p.X() + dir.Y()*p.Y() + dir.Z() * p.Z();
-    // cout << "minval = " << minval << ", val = " << val << " maxval = " << maxval << endl;
-    if (openmin) {
-      if (val < minval+eps) return false;
-    } else {
-      if (val < minval-eps) return false;
-    }
-    if (openmax) {
-      if (val > maxval-eps) return false;
-    } else {
-      if (val > maxval+eps) return false;
-    }
-    return true;
-  }
-};
-
-
-inline gp_Pnt Center (TopoDS_Shape shape)
-{
-  GProp_GProps props;
-  switch (shape.ShapeType())
-    {
-    case TopAbs_FACE:
-      BRepGProp::SurfaceProperties (shape, props); break;
-    default:
-      BRepGProp::LinearProperties(shape, props);
-    }
-  return props.CentreOfMass();
-}
-

libsrc/occ/python_occ_basic.cpp

@@ -1,59 +1,21 @@
 #ifdef NG_PYTHON
 #ifdef OCCGEOMETRY
 
-#include <../general/ngpython.hpp>
+#include <general/ngpython.hpp>
 #include <core/python_ngcore.hpp>
-#include "../meshing/python_mesh.hpp"
+#include <meshing/python_mesh.hpp>
-
 #include <meshing.hpp>
-#include <occgeom.hpp>
 
+#include "occgeom.hpp"
+
+#include <BRepBuilderAPI_Transform.hxx>
 #include <gp_Ax1.hxx>
 #include <gp_Ax2.hxx>
 #include <gp_Ax2d.hxx>
+#include <gp_Ax3.hxx>
 #include <gp_Trsf.hxx>
-#include <BRepPrimAPI_MakeSphere.hxx>
-#include <BRepPrimAPI_MakeCylinder.hxx>
-#include <BRepPrimAPI_MakeBox.hxx>
-#include <BRepPrimAPI_MakePrism.hxx>
-#include <BRepOffsetAPI_MakePipe.hxx>
-#include <BRepAlgoAPI_Cut.hxx>
-#include <BRepAlgoAPI_Common.hxx>
-#include <BRepAlgoAPI_Fuse.hxx>
-// #include <XCAFDoc_VisMaterialTool.hxx>
-#include <TDF_Attribute.hxx>
-#include <Standard_GUID.hxx>
-#include <Geom_TrimmedCurve.hxx>
-#include <GC_MakeSegment.hxx>
-#include <GC_MakeCircle.hxx>
-#include <GC_MakeArcOfCircle.hxx>
-#include <BRepBuilderAPI_MakeEdge.hxx>
-#include <BRepBuilderAPI_MakeWire.hxx>
-#include <BRepBuilderAPI_Transform.hxx>
-#include <BRepBuilderAPI_MakeFace.hxx>
-#include <BRepFilletAPI_MakeFillet.hxx>
-#include <BRepOffsetAPI_ThruSections.hxx>
-
-#include <BRepGProp.hxx>
-#include <BRepOffsetAPI_MakeThickSolid.hxx>
-#include <BRepLib.hxx>
-
-#include <Geom2d_Curve.hxx>
-#include <Geom2d_Ellipse.hxx>
-#include <Geom2d_TrimmedCurve.hxx>
-#include <GCE2d_MakeSegment.hxx>
-#include <GCE2d_MakeCircle.hxx>
-
-#include <python_occ.hpp>
-
-
-#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4
-#define OCC_HAVE_DUMP_JSON
-#endif
-
-
-
 
+using namespace netgen;
 
 DLL_HEADER void ExportNgOCCBasic(py::module &m) 
 {
@@ -87,6 +49,16 @@ DLL_HEADER void ExportNgOCCBasic(py::module &m)
     .def("__sub__", [](gp_Pnt p1, gp_Pnt p2) { return gp_Vec(p2, p1); }) 
     .def("__add__", [](gp_Pnt p, gp_Vec v) { return p.Translated(v); }) // gp_Pnt(p.X()+v.X(), p.Y()+v.Y(), p.Z()+v.Z()); })
     .def("__sub__", [](gp_Pnt p, gp_Vec v) { return p.Translated(-v); }) // gp_Pnt(p.X()-v.X(), p.Y()-v.Y(), p.Z()-v.Z()); })
+    .def("__getitem__", [](const gp_Pnt& p, int index)
+    {
+      if(index == 0)
+        return p.X();
+      if(index == 1)
+        return p.Y();
+      if(index == 2)
+        return p.Z();
+      throw std::out_of_range("Point index must be in range [0,3)!");
+    })
     ;
   
   py::class_<gp_Vec>(m, "gp_Vec", "3d OCC vector")
@@ -297,6 +269,7 @@ DLL_HEADER void ExportNgOCCBasic(py::module &m)
   py::class_<gp_Trsf>(m, "gp_Trsf")
     .def(py::init<>())    
     .def("SetMirror", [] (gp_Trsf & trafo, const gp_Ax1 & ax) { trafo.SetMirror(ax); return trafo; })
+    .def("Inverted", &gp_Trsf::Inverted)
     .def_static("Translation", [] (const gp_Vec & v) { gp_Trsf trafo; trafo.SetTranslation(v); return trafo; })
     .def_static("Scale", [] (const gp_Pnt & p, double s) { gp_Trsf trafo; trafo.SetScale(p,s); return trafo; })    
     .def_static("Mirror", [] (const gp_Ax1 & ax) { gp_Trsf trafo; trafo.SetMirror(ax); return trafo; })
@@ -308,7 +281,9 @@ DLL_HEADER void ExportNgOCCBasic(py::module &m)
                 { gp_Trsf trafo; trafo.SetTransformation(from, to); return trafo; })
     .def(py::self * py::self)
     .def("__call__", [] (gp_Trsf & trafo, const TopoDS_Shape & shape) {
-        return BRepBuilderAPI_Transform(shape, trafo).Shape();
+        BRepBuilderAPI_Transform builder(shape, trafo, true);
+        PropagateProperties(builder, shape, occ2ng(trafo));
+        return builder.Shape();
       })
     .def("__str__", [](gp_Trsf & trafo)
     {

libsrc/occ/python_occ_shapes.cpp

@@ -3,185 +3,76 @@
 
 #include <regex>
 
-#include <../general/ngpython.hpp>
+#include <general/ngpython.hpp>
 #include <core/python_ngcore.hpp>
-#include "../meshing/python_mesh.hpp"
+#include <meshing/python_mesh.hpp>
-
 #include <meshing.hpp>
-#include <occgeom.hpp>
 
-#include <gp_Ax1.hxx>
+#include "occgeom.hpp"
-#include <gp_Ax2.hxx>
+
-#include <gp_Ax2d.hxx>
+#include <BOPAlgo_Builder.hxx>
-#include <gp_Trsf.hxx>
+#include <BOPTools_AlgoTools.hxx>
-#include <BRepPrimAPI_MakeSphere.hxx>
-#include <BRepPrimAPI_MakeCylinder.hxx>
-#include <BRepPrimAPI_MakeRevol.hxx>
-#include <BRepPrimAPI_MakeBox.hxx>
-#include <BRepPrimAPI_MakePrism.hxx>
-#include <BRepPrimAPI_MakeHalfSpace.hxx>
-#include <BRepOffsetAPI_MakePipe.hxx>
-#include <BRepOffsetAPI_MakePipeShell.hxx>
-#include <BRepAlgoAPI_Cut.hxx>
 #include <BRepAlgoAPI_Common.hxx>
+#include <BRepAlgoAPI_Cut.hxx>
 #include <BRepAlgoAPI_Fuse.hxx>
-// #include <XCAFDoc_VisMaterialTool.hxx>
+#include <BRepAlgo_NormalProjection.hxx>
-#include <TDF_Attribute.hxx>
-#include <TDataStd_Real.hxx>
-#include <TDataStd_Name.hxx>
-#include <Standard_GUID.hxx>
-#include <Geom_TrimmedCurve.hxx>
-#include <Geom_Plane.hxx>
-#include <Geom_BSplineCurve.hxx>
-#include <Geom2d_BSplineCurve.hxx>
-#include <Geom_BezierCurve.hxx>
-#include <GeomAPI_PointsToBSpline.hxx>
-#include <Geom2dAPI_PointsToBSpline.hxx>
-#include <GeomAbs_Shape.hxx>
-#include <Geom_BSplineSurface.hxx>
-#include <GeomAPI_PointsToBSplineSurface.hxx>
-#include <GeomAPI_Interpolate.hxx>
-#include <Geom2dAPI_Interpolate.hxx>
-#include <GC_MakeSegment.hxx>
-#include <GC_MakeCircle.hxx>
-#include <GC_MakeArcOfCircle.hxx>
 #include <BRepBuilderAPI_MakeEdge.hxx>
-#include <BRepBuilderAPI_MakeEdge2d.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepBuilderAPI_MakeVertex.hxx>
 #include <BRepBuilderAPI_MakeWire.hxx>
 #include <BRepBuilderAPI_Transform.hxx>
-#include <BRepBuilderAPI_MakeFace.hxx>
-#include <BRepFilletAPI_MakeFillet.hxx>
-#include <BRepFilletAPI_MakeChamfer.hxx>
-#include <BRepOffsetAPI_ThruSections.hxx>
-#include <BRepOffsetAPI_MakeOffset.hxx>
 #include <BRepExtrema_DistShapeShape.hxx>
-
+#include <BRepFilletAPI_MakeChamfer.hxx>
+#include <BRepFilletAPI_MakeFillet.hxx>
 #include <BRepGProp.hxx>
-#include <BRepOffsetAPI_MakeThickSolid.hxx>
+#include <BRepLProp_SLProps.hxx>
 #include <BRepLib.hxx>
-
+#include <BRepMesh_IncrementalMesh.hxx>
+#include <BRepOffsetAPI_MakeOffset.hxx>
+#include <BRepOffsetAPI_MakePipe.hxx>
+#include <BRepOffsetAPI_MakePipeShell.hxx>
+#include <BRepOffsetAPI_MakeThickSolid.hxx>
+#include <BRepOffsetAPI_ThruSections.hxx>
+#include <BRepPrimAPI_MakeBox.hxx>
+#include <BRepPrimAPI_MakeCylinder.hxx>
+#include <BRepPrimAPI_MakeHalfSpace.hxx>
+#include <BRepPrimAPI_MakePrism.hxx>
+#include <BRepPrimAPI_MakeRevol.hxx>
+#include <BRepPrimAPI_MakeSphere.hxx>
+#include <BRepTools.hxx>
+#include <GCE2d_MakeArcOfCircle.hxx>
+#include <GCE2d_MakeCircle.hxx>
+#include <GCE2d_MakeSegment.hxx>
+#include <GC_MakeArcOfCircle.hxx>
+#include <GC_MakeCircle.hxx>
+#include <GC_MakeSegment.hxx>
+#include <GProp_GProps.hxx>
+#include <Geom2d_BSplineCurve.hxx>
 #include <Geom2d_Curve.hxx>
 #include <Geom2d_Ellipse.hxx>
 #include <Geom2d_TrimmedCurve.hxx>
-#include <GCE2d_MakeSegment.hxx>
+#include <Geom2dAPI_Interpolate.hxx>
-#include <GCE2d_MakeCircle.hxx>
+#include <Geom2dAPI_PointsToBSpline.hxx>
-#include <GCE2d_MakeArcOfCircle.hxx>
+#include <GeomAPI_Interpolate.hxx>
-#include <ShapeUpgrade_UnifySameDomain.hxx>
+#include <GeomAPI_PointsToBSpline.hxx>
+#include <GeomAPI_PointsToBSplineSurface.hxx>
 #include <GeomLProp_SLProps.hxx>
-
+#include <Geom_BezierCurve.hxx>
-#include <BOPTools_AlgoTools.hxx>
+#include <Geom_BSplineCurve.hxx>
+#include <Geom_BSplineSurface.hxx>
+#include <Geom_Plane.hxx>
+#include <Geom_TrimmedCurve.hxx>
 #include <IntTools_Context.hxx>
-#include <STEPControl_Writer.hxx>
 #include <ShapeAnalysis_FreeBounds.hxx>
-
+#include <ShapeUpgrade_UnifySameDomain.hxx>
-
+#include <gp_Ax1.hxx>
-#include <python_occ.hpp>
+#include <gp_Ax2.hxx>
-
+#include <gp_Ax2d.hxx>
-#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4
+#include <gp_Pln.hxx>
-#define OCC_HAVE_DUMP_JSON
+#include <gp_Trsf.hxx>
-#endif
 
 using namespace netgen;
 
-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:
-  TopoDS_Shape Max(gp_Vec dir)
-  {
-    double maxval = -1e99;
-    TopoDS_Shape maxshape;
-    for (auto shape : *this)
-      {
-        GProp_GProps props;
-        gp_Pnt center;
-        
-        switch (shape.ShapeType())
-          {
-          case TopAbs_VERTEX:
-            center = BRep_Tool::Pnt (TopoDS::Vertex(shape)); break;
-          case TopAbs_FACE:
-            BRepGProp::SurfaceProperties (shape, props);
-            center = props.CentreOfMass();
-            break;
-          default:
-            BRepGProp::LinearProperties(shape, props);
-            center = props.CentreOfMass();
-          }
-        
-        double val = center.X()*dir.X() + center.Y()*dir.Y() + center.Z() * dir.Z();
-        if (val > maxval)
-          {
-            maxval = val;
-            maxshape = shape;
-          }
-      }
-    return maxshape;
-  }
-
-  TopoDS_Shape Nearest(gp_Pnt pnt)
-  {
-    double mindist = 1e99;
-    TopoDS_Shape nearestshape;
-    auto vertex = BRepBuilderAPI_MakeVertex (pnt).Vertex();
-    
-    for (auto shape : *this)
-      {
-        double dist = BRepExtrema_DistShapeShape(shape, vertex).Value();
-        // cout << "dist = " << dist << endl;
-        if (dist < mindist)
-          nearestshape = shape;
-      }
-    return nearestshape;
-  }
-  
-  ListOfShapes SubShapes(TopAbs_ShapeEnum type) const
-  {
-    std::set<TopoDS_Shape, ShapeLess> unique_shapes;
-    for(const auto& shape : *this)
-      for(TopExp_Explorer e(shape, type); e.More(); e.Next())
-        unique_shapes.insert(e.Current());
-    ListOfShapes sub;
-    for(const auto& shape : unique_shapes)
-      sub.push_back(shape);
-    return sub;
-  }
-  ListOfShapes Solids() const
-  {
-    return SubShapes(TopAbs_SOLID);
-  }
-  ListOfShapes Faces() const
-  {
-    return SubShapes(TopAbs_FACE);
-  }
-  ListOfShapes Edges() const
-  {
-    return SubShapes(TopAbs_EDGE);
-  }
-  ListOfShapes Vertices() const
-  {
-    return SubShapes(TopAbs_VERTEX);
-  }
-
-  ListOfShapes operator*(const ListOfShapes& other) const
-  {
-    ListOfShapes common;
-    for(const auto& shape : (*this))
-      for(const auto& shape_o : other)
-        if(shape.IsSame(shape_o))
-          common.push_back(shape);
-    return common;
-  }
-};
-
-
 void ExtractEdgeData( const TopoDS_Edge & edge, int index, std::vector<double> * p, Box<3> & box )
 {
     if (BRep_Tool::Degenerated(edge)) return;
@@ -312,23 +203,6 @@ py::object CastShape(const TopoDS_Shape & s)
 };
 
 
-template <class TBuilder>
-void PropagateProperties (TBuilder & builder, TopoDS_Shape shape)
-{
-  // #ifdef OCC_HAVE_HISTORY  
-  // Handle(BRepTools_History) history = builder.History();
-  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()];
-        // for (auto mods : history->Modified(e.Current()))
-        for (auto mods : builder.Modified(e.Current()))
-          OCCGeometry::global_shape_properties[mods.TShape()].Merge(prop);
-      }
-  // #endif  
-}
-
-
 class WorkPlane : public enable_shared_from_this<WorkPlane>
 {
   gp_Ax3 axes;
@@ -819,36 +693,21 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
            return sub;
          }, py::arg("type"), "returns list of sub-shapes of type 'type'")
     
-    .def_property_readonly("solids", [] (const TopoDS_Shape & shape)
+    .def_property_readonly("solids", GetSolids,
+            "returns all sub-shapes of type 'SOLID'")
+    .def_property_readonly("faces", GetFaces,
+            "returns all sub-shapes of type 'FACE'")
+    .def_property_readonly("edges", GetEdges,
+            "returns all sub-shapes of type 'EDGE'")
+    .def_property_readonly("wires", GetWires,
+                           "returns all sub-shapes of type 'WIRE'")
+    .def_property_readonly("vertices", GetVertices,
+            "returns all sub-shapes of type 'VERTEX'")
+    .def_property_readonly("bounding_box", [] ( const TopoDS_Shape &shape )
             {
-      ListOfShapes solids;
+               auto box = GetBoundingBox(shape);
-      for(TopExp_Explorer e(shape, TopAbs_SOLID); e.More(); e.Next())
+               return py::make_tuple( ng2occ(box.PMin()), ng2occ(box.PMax()) );
-        solids.push_back(e.Current());
+            }, "returns bounding box (pmin, pmax)")
-      return solids;
-    }, "returns all sub-shapes of type 'SOLID'")
-    .def_property_readonly("faces", [] (const TopoDS_Shape & shape)
-         {
-           ListOfShapes sub;
-           for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next())
-             sub.push_back(e.Current());
-           return sub;
-         }, "returns all sub-shapes of type 'FACE'")
-    
-    .def_property_readonly("edges", [] (const TopoDS_Shape & shape)
-         {
-           ListOfShapes sub;
-           for (TopExp_Explorer e(shape, TopAbs_EDGE); e.More(); e.Next())
-             sub.push_back(e.Current());
-           return sub;
-         }, "returns all sub-shapes of type 'EDGE'")
-    
-    .def_property_readonly("vertices", [] (const TopoDS_Shape & shape)
-         {
-           ListOfShapes sub;
-           for (TopExp_Explorer e(shape, TopAbs_VERTEX); e.More(); e.Next())
-             sub.push_back(e.Current());
-           return sub;
-         }, "returns all sub-shapes of type 'VERTEX'")
 
     .def("Properties", [] (const TopoDS_Shape & shape)
          {
@@ -896,8 +755,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
            // version 1: Transoformation
            gp_Trsf trafo;
            trafo.SetTranslation(v);
-           BRepBuilderAPI_Transform builder(shape, trafo);
+           BRepBuilderAPI_Transform builder(shape, trafo, true);
-           PropagateProperties(builder, shape);
+           PropagateProperties(builder, shape, occ2ng(trafo));
            return builder.Shape();
            // version 2: change location
            // ...
@@ -908,8 +767,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
          {
            gp_Trsf trafo;
            trafo.SetRotation(ax, ang*M_PI/180);            
-           BRepBuilderAPI_Transform builder(shape, trafo);
+           BRepBuilderAPI_Transform builder(shape, trafo, true);
-           PropagateProperties(builder, shape);
+           PropagateProperties(builder, shape, occ2ng(trafo));
            return builder.Shape();
          }, py::arg("axis"), py::arg("ang"),
          "copy shape, and rotet copy by 'ang' degrees around 'axis'")
@@ -918,8 +777,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
          {
            gp_Trsf trafo;
            trafo.SetMirror(ax.Ax2());
-           BRepBuilderAPI_Transform builder(shape, trafo);
+           BRepBuilderAPI_Transform builder(shape, trafo, true);
-           PropagateProperties(builder, shape);
+           PropagateProperties(builder, shape, occ2ng(trafo));
            return builder.Shape();
          }, py::arg("axes"),
          "copy shape, and mirror over plane defined by 'axes'")
@@ -928,8 +787,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
          {
            gp_Trsf trafo;
            trafo.SetMirror(ax);
-           BRepBuilderAPI_Transform builder(shape, trafo);
+           BRepBuilderAPI_Transform builder(shape, trafo, true);
-           PropagateProperties(builder, shape);
+           PropagateProperties(builder, shape, occ2ng(trafo));
            return builder.Shape();
          }, py::arg("axes"),
          "copy shape, and mirror around axis 'axis'")
@@ -938,8 +797,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
          {
            gp_Trsf trafo;
            trafo.SetScale(p, s);
-           BRepBuilderAPI_Transform builder(shape, trafo);
+           BRepBuilderAPI_Transform builder(shape, trafo, true);
-           PropagateProperties(builder, shape);
+           PropagateProperties(builder, shape, occ2ng(trafo));
            return builder.Shape();
          }, py::arg("p"), py::arg("s"),
          "copy shape, and scale copy by factor 's'")
@@ -967,7 +826,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
           return *name;
         else
           return string();
-      }, [](const TopoDS_Shape & self, string name) {
+      }, [](const TopoDS_Shape & self, optional<string> name) {
         OCCGeometry::global_shape_properties[self.TShape()].name = name;            
       }, "'name' of shape")
     
@@ -1138,16 +997,24 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
     .def("Reversed", [](const TopoDS_Shape & shape) {
         return CastShape(shape.Reversed()); })
 
-    .def("Extrude", [](const TopoDS_Shape & shape, double h) {
+    .def("Extrude", [](const TopoDS_Shape & shape, double h,
+                       optional<gp_Vec> dir) {
         for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next())
           {
             Handle(Geom_Surface) surf = BRep_Tool::Surface (TopoDS::Face(e.Current()));
+            gp_Vec edir;
+            if(dir.has_value())
+              edir = *dir;
+            else
+              {
                 gp_Vec du, dv;
                 gp_Pnt p;
                 surf->D1 (0,0,p,du,dv);
-            BRepPrimAPI_MakePrism builder(shape, h*du^dv);
+                edir = du^dv;
+              }
+            BRepPrimAPI_MakePrism builder(shape, h*edir, true);
 
-            for (auto typ : { TopAbs_EDGE, TopAbs_VERTEX })
+            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()];
@@ -1158,7 +1025,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
             return builder.Shape();
           }
         throw Exception("no face found for extrusion");
-      }, py::arg("h"), "extrude shape to thickness 'h', shape must contain a plane surface")
+    }, py::arg("h"), py::arg("dir")=nullopt, "extrude shape to thickness 'h', shape must contain a plane surface, optionally give an extrusion direction")
     
     .def("Extrude", [] (const TopoDS_Shape & face, gp_Vec vec) {
         return BRepPrimAPI_MakePrism (face, vec).Shape();
@@ -1168,7 +1035,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
       // for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next())
         {
           // return BRepPrimAPI_MakeRevol (shape, A, D*M_PI/180).Shape();
-          BRepPrimAPI_MakeRevol builder(shape, A, D*M_PI/180);
+          BRepPrimAPI_MakeRevol builder(shape, A, D*M_PI/180, true);
             
           for (auto typ : { TopAbs_EDGE, TopAbs_VERTEX })
             for (TopExp_Explorer e(shape, typ); e.More(); e.Next())
@@ -1229,30 +1096,17 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
            BRepMesh_IncrementalMesh (shape, deflection, true);
          })
 
-    .def("Identify", [](const TopoDS_Shape & me, const TopoDS_Shape & you, string name) {
-        // only edges supported, by now
-        auto me_edge = TopoDS::Edge(me);
-        auto you_edge = TopoDS::Edge(you);
 
-        GProp_GProps props;
+    .def("Identify", py::overload_cast<const TopoDS_Shape &, const TopoDS_Shape &, string, Identifications::ID_TYPE, std::optional<gp_Trsf>>(&Identify),
-        BRepGProp::LinearProperties(me, props);
+            py::arg("other"), py::arg("name"),
-        gp_Pnt cme = props.CentreOfMass();
+            py::arg("type")=Identifications::PERIODIC, py::arg("trafo")=nullopt,
-        BRepGProp::LinearProperties(you, props);
+            "Identify shapes for periodic meshing")
-        gp_Pnt cyou = props.CentreOfMass();
-
-        double s0, s1;
-        auto curve_me = BRep_Tool::Curve(me_edge, s0, s1);
-        auto vme = occ2ng(curve_me->Value(s1))-occ2ng(curve_me->Value(s0));
-        auto curve_you = BRep_Tool::Curve(you_edge, s0, s1);
-        auto vyou = occ2ng(curve_you->Value(s1))-occ2ng(curve_you->Value(s0));
-        
-        bool inv = vme*vyou < 0;
-        OCCGeometry::identifications[me.TShape()].push_back
-                              (OCCIdentification { you, Transformation<3>(occ2ng(cyou) - occ2ng(cme)), inv, name });
-        OCCGeometry::identifications[you.TShape()].push_back
-                              (OCCIdentification { me, Transformation<3>(occ2ng(cme) - occ2ng(cyou)), inv, name });
-      }, py::arg("other"), py::arg("name"), "Identify shapes for periodic meshing")
 
+    .def("Distance", [](const TopoDS_Shape& self,
+                        const TopoDS_Shape& other)
+    {
+      return BRepExtrema_DistShapeShape(self, other).Value();
+    })
     
     .def("Triangulation", [](const TopoDS_Shape & shape)
          {
@@ -1582,6 +1436,14 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
         auto ax = gp_Ax3(p, du^dv, du);
         return make_shared<WorkPlane> (ax);
       })
+    .def("ProjectWire", [](const TopoDS_Face& face,
+                           const TopoDS_Wire& wire)
+    {
+      BRepAlgo_NormalProjection builder(face);
+      builder.Add(wire);
+      builder.Build();
+      return builder.Projection();
+    })
     ;
   py::class_<TopoDS_Solid, TopoDS_Shape> (m, "Solid");
   
@@ -1685,6 +1547,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
          })
     .def_property_readonly("solids", &ListOfShapes::Solids)
     .def_property_readonly("faces", &ListOfShapes::Faces)
+    .def_property_readonly("wires", &ListOfShapes::Wires)
     .def_property_readonly("edges", &ListOfShapes::Edges)
     .def_property_readonly("vertices", &ListOfShapes::Vertices)
     .def(py::self * py::self)
@@ -1731,12 +1594,15 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
     .def("Nearest", [] (ListOfShapes & shapes, gp_Pnt pnt) 
          { return CastShape(shapes.Nearest(pnt)); },
          py::arg("p"), "returns shape nearest to point 'p'")
+    .def("Nearest", [] (ListOfShapes & shapes, gp_Pnt2d pnt) 
+         { return CastShape(shapes.Nearest( { pnt.X(), pnt.Y(), 0 })); },
+         py::arg("p"), "returns shape nearest to point 'p'")
     
     .def_property("name", [](ListOfShapes& shapes)
     {
       throw Exception("Cannot get property of ListOfShapes, get the property from individual shapes!");
     },
-      [](ListOfShapes& shapes, std::string name)
+      [](ListOfShapes& shapes, optional<std::string> name)
       {
         for(auto& shape : shapes)
           {
@@ -1764,6 +1630,23 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
             OCCGeometry::global_shape_properties[shape.TShape()].maxh = maxh;
           }
       }, "set maxh for all elements of list")
+    .def_property("hpref", [](ListOfShapes& shapes)
+    {
+      throw Exception("Cannot get property of ListOfShapes, get the property from individual shapes!");
+    },
+      [](ListOfShapes& shapes, double hpref)
+      {
+        for(auto& shape : shapes)
+          {
+            auto& val = OCCGeometry::global_shape_properties[shape.TShape()].hpref;
+            val = max2(hpref, val);
+          }
+      }, "set hpref for all elements of list")
+    
+    .def("Identify", py::overload_cast<const ListOfShapes&, const ListOfShapes&, string, Identifications::ID_TYPE, gp_Trsf>(&Identify),
+            py::arg("other"), py::arg("name"),
+            py::arg("type")=Identifications::PERIODIC, py::arg("trafo"),
+            "Identify shapes for periodic meshing")
 
     ;
          
@@ -1873,13 +1756,13 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
     "create box with opposite points 'p1' and 'p2'");
 
   m.def("Prism", [] (const TopoDS_Shape & face, gp_Vec vec) {
-      return BRepPrimAPI_MakePrism (face, vec).Shape();
+      return BRepPrimAPI_MakePrism (face, vec, true).Shape();
     }, py::arg("face"), py::arg("v"),
     "extrude face along the vector 'v'");
 
   m.def("Revolve", [] (const TopoDS_Shape & face,const gp_Ax1 &A, const double D) {
-      //comvert angle from deg to rad
+      //convert angle from deg to rad
-      return BRepPrimAPI_MakeRevol (face, A, D*M_PI/180).Shape();
+      return BRepPrimAPI_MakeRevol (face, A, D*M_PI/180, true).Shape();
     });
 
   m.def("Pipe", [] (const TopoDS_Wire & spine, const TopoDS_Shape & profile,

libsrc/occ/vsocc.cpp

@@ -8,19 +8,18 @@
 
 #include <occgeom.hpp>
 
-#include "TopoDS_Shape.hxx"
+#include <BRepAdaptor_Surface.hxx>
-#include "TopoDS_Vertex.hxx"
+#include <BRepBndLib.hxx>
-#include "TopExp_Explorer.hxx"
+#include <BRepLProp_SLProps.hxx>
-#include "BRep_Tool.hxx"
+#include <BRep_Tool.hxx>
-#include "TopoDS.hxx"
+#include <Bnd_Box.hxx>
-#include "gp_Pnt.hxx"
+#include <Geom_Curve.hxx>
-#include "Geom_Curve.hxx"
+#include <Poly_PolygonOnTriangulation.hxx>
-#include "Poly_Triangulation.hxx"
+#include <Poly_Triangle.hxx>
-#include "Poly_Array1OfTriangle.hxx"
+#include <Poly_Triangulation.hxx>
-#include "TColgp_Array1OfPnt2d.hxx"
+#include <TopoDS.hxx>
-#include "Poly_Triangle.hxx"
+#include <TopoDS_Edge.hxx>
-#include "Poly_Polygon3D.hxx"
+#include <gp_Pnt.hxx>
-#include "Poly_PolygonOnTriangulation.hxx"
 
 #include <visual.hpp>
 
@@ -87,7 +86,7 @@ namespace netgen
       // Added clipping planes to Geometry view
       SetClippingPlane();
 
-      GLfloat matcoledge[] = {  0, 0, 1, 1};
+      GLfloat matcoledge[] = {  0, 0, 0, 1};
       GLfloat matcolhiedge[] = {  1, 0, 0, 1};
 
       glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcoledge);

nglib/CMakeLists.txt

@@ -9,7 +9,6 @@ if(WIN32)
           $<TARGET_OBJECTS:csg>
 
           $<TARGET_OBJECTS:visual>
-          $<TARGET_OBJECTS:occ>
        )
        if(USE_GUI)
            set(nglib_objects ${nglib_objects}
@@ -18,6 +17,9 @@ if(WIN32)
                $<TARGET_OBJECTS:csgvis>
            )
        endif(USE_GUI)
+       if(USE_OCC)
+           set(nglib_objects ${nglib_objects} $<TARGET_OBJECTS:occ>)
+       endif(USE_OCC)
 endif(WIN32)
 
 add_library(nglib SHARED nglib.cpp ${nglib_objects})

nglib/nglib.cpp

@@ -880,7 +880,7 @@ namespace nglib
 
       mp->Transfer_Parameters();
 
-      OCCFindEdges(*occgeom, *me, mparam);
+      occgeom->FindEdges(*me, mparam);
 
       if((me->GetNP()) && (me->GetNFD()))
       {
@@ -910,11 +910,6 @@ namespace nglib
       // parameters structure
       mp->Transfer_Parameters();
 
-
-      // Only go into surface meshing if the face descriptors have already been added
-      if(!me->GetNFD())
-         return NG_ERROR;
-
       numpoints = me->GetNP();
 
       // Initially set up only for surface meshing without any optimisation
@@ -926,8 +921,8 @@ namespace nglib
          perfstepsend = MESHCONST_OPTSURFACE;
       }
 
-      OCCMeshSurface(*occgeom, *me, mparam);
+      occgeom->MeshSurface(*me, mparam);
-      OCCOptimizeSurface(*occgeom, *me, mparam);
+      occgeom->OptimizeSurface(*me, mparam);
 
       me->CalcSurfacesOfNode();
       

python/occ.py

@@ -9,6 +9,10 @@ For more detailed documentation consult the OCCT docs, a good starting point is
 https://dev.opencascade.org/doc/refman/html/class_b_rep_builder_a_p_i___make_shape.html
 """
 
+from .config import USE_OCC
+if not USE_OCC:
+    raise ImportError("Netgen was not built with Opencascade support")
+
 from .libngpy._NgOCC import *
 from .meshing import meshsize
 

tests/dockerfile

@@ -11,12 +11,15 @@ RUN apt-get update && apt-get -y install \
         libcgns-dev \
         libglu1-mesa-dev \
         libhdf5-dev \
+        libocct-ocaf-dev \
+        libocct-visualization-dev \
         libocct-data-exchange-dev \
         libocct-draw-dev \
         libpython3-dev \
         libtbb-dev \
         libxi-dev \
         libxmu-dev \
+        occt-misc \
         python3 \
         python3-distutils \
         python3-numpy \

tests/dockerfile_mpi

@@ -1,6 +1,27 @@
 FROM ubuntu:20.04
 ENV DEBIAN_FRONTEND=noninteractive
 MAINTAINER Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>
-RUN apt-get update && apt-get -y install python3 libpython3-dev python3-pip libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-numpy python3-tk python3-mpi4py clang-tidy python3-distutils clang libopenmpi-dev openmpi-bin gfortran
+RUN apt-get update && apt-get -y install \
+        ccache \
+        clang \
+        clang-tidy \
+        cmake \
+        g++ \
+        gfortran \
+        git \
+        libglu1-mesa-dev \
+        libopenmpi-dev \
+        libpython3-dev \
+        libxmu-dev \
+        openmpi-bin \
+        python3 \
+        python3-distutils \
+        python3-mpi4py \
+        python3-numpy \
+        python3-pip \
+        python3-tk \
+        tcl-dev \
+        tk-dev
+
 RUN python3 -m pip install pytest-mpi pytest-check pytest
 ADD . /root/src/netgen

tests/pytest/compare_results.py

@@ -73,6 +73,13 @@ for bad1,bad2, f1, f2 in zip(data['badness'], data2['badness'], data['file'], da
     if bad2>0 and bad2<0.9*bad1:
         print(f"file {f1} got better: {bad1} -> {bad2}")
 
+for bad1,bad2, f1, f2 in zip(data['#trigs'], data2['#trigs'], data['file'], data2['file']):
+    assert f1==f2
+    if bad2>0 and bad2>1.1*bad1:
+        print(f"file {f1} got worse:  {bad1} -> {bad2}")
+    if bad2>0 and bad2<0.9*bad1:
+        print(f"file {f1} got better: {bad1} -> {bad2}")
+
 n = len(data)+1
 fig,ax = plt.subplots(figsize=(10,7))
 for i,d in enumerate(['min trig angle','min tet angle','max trig angle','max tet angle']):

tests/pytest/results.json

@@ -1590,8 +1590,8 @@
                 0.0,
                 0.0
             ],
-            "ne1d": 81,
+            "ne1d": 80,
-            "ne2d": 436,
+            "ne2d": 412,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -1606,7 +1606,7 @@
                 0.0
             ],
             "ne1d": 86,
-            "ne2d": 452,
+            "ne2d": 440,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -1621,7 +1621,7 @@
                 0.0
             ],
             "ne1d": 86,
-            "ne2d": 450,
+            "ne2d": 464,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -1635,8 +1635,8 @@
                 0.0,
                 0.0
             ],
-            "ne1d": 81,
+            "ne1d": 80,
-            "ne2d": 436,
+            "ne2d": 412,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -1650,8 +1650,8 @@
                 0.0,
                 0.0
             ],
-            "ne1d": 82,
+            "ne1d": 83,
-            "ne2d": 441,
+            "ne2d": 453,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -1665,8 +1665,8 @@
                 0.0,
                 0.0
             ],
-            "ne1d": 86,
+            "ne1d": 84,
-            "ne2d": 472,
+            "ne2d": 462,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -2790,8 +2790,38 @@
                 0.0,
                 0.0
             ],
-            "ne1d": 31,
+            "ne1d": 27,
-            "ne2d": 97,
+            "ne2d": 87,
+            "ne3d": 0,
+            "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
+            "total_badness": 0.0
+        },
+        {
+            "angles_tet": [
+                0.0,
+                0.0
+            ],
+            "angles_trig": [
+                0.0,
+                0.0
+            ],
+            "ne1d": 28,
+            "ne2d": 80,
+            "ne3d": 0,
+            "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
+            "total_badness": 0.0
+        },
+        {
+            "angles_tet": [
+                0.0,
+                0.0
+            ],
+            "angles_trig": [
+                0.0,
+                0.0
+            ],
+            "ne1d": 26,
+            "ne2d": 78,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -2806,7 +2836,7 @@
                 0.0
             ],
             "ne1d": 27,
-            "ne2d": 71,
+            "ne2d": 87,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -2820,8 +2850,8 @@
                 0.0,
                 0.0
             ],
-            "ne1d": 25,
+            "ne1d": 36,
-            "ne2d": 85,
+            "ne2d": 150,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -2835,38 +2865,8 @@
                 0.0,
                 0.0
             ],
-            "ne1d": 31,
+            "ne1d": 42,
-            "ne2d": 97,
+            "ne2d": 246,
-            "ne3d": 0,
-            "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
-            "total_badness": 0.0
-        },
-        {
-            "angles_tet": [
-                0.0,
-                0.0
-            ],
-            "angles_trig": [
-                0.0,
-                0.0
-            ],
-            "ne1d": 41,
-            "ne2d": 173,
-            "ne3d": 0,
-            "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
-            "total_badness": 0.0
-        },
-        {
-            "angles_tet": [
-                0.0,
-                0.0
-            ],
-            "angles_trig": [
-                0.0,
-                0.0
-            ],
-            "ne1d": 30,
-            "ne2d": 108,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -2883,7 +2883,7 @@
                 0.0
             ],
             "ne1d": 32,
-            "ne2d": 148,
+            "ne2d": 142,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -2913,7 +2913,7 @@
                 0.0
             ],
             "ne1d": 32,
-            "ne2d": 134,
+            "ne2d": 126,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -2928,7 +2928,7 @@
                 0.0
             ],
             "ne1d": 32,
-            "ne2d": 148,
+            "ne2d": 142,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -2942,8 +2942,8 @@
                 0.0,
                 0.0
             ],
-            "ne1d": 43,
+            "ne1d": 42,
-            "ne2d": 300,
+            "ne2d": 326,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -2957,8 +2957,8 @@
                 0.0,
                 0.0
             ],
-            "ne1d": 79,
+            "ne1d": 76,
-            "ne2d": 821,
+            "ne2d": 811,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -2975,7 +2975,7 @@
                 0.0
             ],
             "ne1d": 32,
-            "ne2d": 124,
+            "ne2d": 118,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -3005,7 +3005,7 @@
                 0.0
             ],
             "ne1d": 32,
-            "ne2d": 110,
+            "ne2d": 102,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -3020,7 +3020,7 @@
                 0.0
             ],
             "ne1d": 32,
-            "ne2d": 124,
+            "ne2d": 118,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -3034,8 +3034,8 @@
                 0.0,
                 0.0
             ],
-            "ne1d": 43,
+            "ne1d": 42,
-            "ne2d": 249,
+            "ne2d": 274,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0
@@ -3049,8 +3049,8 @@
                 0.0,
                 0.0
             ],
-            "ne1d": 79,
+            "ne1d": 76,
-            "ne2d": 687,
+            "ne2d": 672,
             "ne3d": 0,
             "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]",
             "total_badness": 0.0

tests/pytest/test_occ.py

@@ -0,0 +1,43 @@
+import pytest
+import math
+from pytest import approx
+from pytest_check import check
+
+def check_volume(shape, volume, dim=3):
+    from netgen.occ import OCCGeometry
+    geo = OCCGeometry(shape, dim=dim)
+
+    m = geo.GenerateMesh()
+    assert len(m.Elements2D()) > 0
+    assert len(m.Elements1D()) > 0
+    if dim==3:
+        assert len(m.Elements3D()) > 0
+
+    ngs = pytest.importorskip("ngsolve")
+    mesh = ngs.Mesh(m)
+    mesh.Curve(5)
+    with check: assert ngs.Integrate(1.0, mesh) == approx(volume)
+
+def test_rect_with_two_holes():
+    occ = pytest.importorskip("netgen.occ")
+
+    face = occ.WorkPlane().Rectangle(7,4) \
+    .Circle(2,2,1).Reverse() \
+    .Circle(5,2,1).Reverse().Face()
+    check_volume(face, 7*4-2*math.pi, 2)
+
+def test_unit_square():
+    occ = pytest.importorskip("netgen.occ")
+    check_volume(occ.unit_square.shape, 1, dim=2)
+
+def test_box_and_cyl():
+    occ = pytest.importorskip("netgen.occ")
+    box = occ.Box(occ.Pnt(0,0,0), occ.Pnt(1,1,1))
+    check_volume(box, 1)
+    r = 0.3
+    h = 0.5
+    vcyl = r*r*math.pi*h
+    cyl = occ.Cylinder(occ.Pnt(1,0.5,0.5), occ.X, r=r, h=h)
+    check_volume(cyl, vcyl)
+    fused = box+cyl
+    check_volume(fused, 1+vcyl)

tests/pytest/test_occ_identifications.py

@@ -0,0 +1,76 @@
+import pytest
+
+from netgen.meshing import IdentificationType
+idtype = IdentificationType.CLOSESURFACES
+
+def test_two_boxes():
+    occ = pytest.importorskip("netgen.occ")
+    inner = occ.Box((0,0,0), (1,1,1))
+    trafo = occ.gp_Trsf().Scale(inner.center, 1.1)
+    outer = trafo(inner)
+
+    inner.Identify(outer, "", idtype, trafo)
+    shape = occ.Glue([outer-inner, inner])
+
+    geo = occ.OCCGeometry(shape)
+    mesh = geo.GenerateMesh(maxh=0.3)
+    have_prisms = False
+
+    for el in mesh.Elements3D():
+        if len(el.vertices)==6:
+            have_prisms = True
+            break
+
+    assert have_prisms
+
+def test_two_circles():
+    occ = pytest.importorskip("netgen.occ")
+    circ1 = occ.WorkPlane().Circle(1).Face()
+    trafo = occ.gp_Trsf().Scale(circ1.center, 1.1)
+
+    circ2 = trafo(circ1)
+    circ1.edges[0].Identify(circ2.edges[0], "", idtype, trafo)
+    circ2 -= circ1
+    shape = occ.Glue([circ1, circ2])
+
+    geo = occ.OCCGeometry(shape, 2)
+    mesh = geo.GenerateMesh(maxh=0.2)
+    have_quads = False
+
+    for el in mesh.Elements2D():
+        if len(el.vertices)==4:
+            have_quads = True
+            break
+
+    assert have_quads
+
+def test_cut_identified_face():
+    occ = pytest.importorskip("netgen.occ")
+    from netgen.occ import Z, Box, Cylinder, Glue, OCCGeometry
+    box = Box((-1,-1,0), (1,1,1))
+    cyl = Cylinder( (0,0,0), Z, 0.5, 1 )
+
+    box.faces.Min(Z).Identify(box.faces.Max(Z), "", idtype)
+    shape = Glue([cyl, box])
+    geo = OCCGeometry(shape)
+    mesh = geo.GenerateMesh(maxh=0.5)
+
+    for el in mesh.Elements3D():
+        assert len(el.vertices)==6
+
+def test_identify_multiple_faces():
+    occ = pytest.importorskip("netgen.occ")
+    from netgen.occ import Z, Box, Cylinder, Glue, OCCGeometry, gp_Trsf
+    box = Box((-1,-1,0), (1,1,1))
+    cyl = Cylinder( (0,0,0), Z, 0.5, 1 )
+
+    shape = Glue([box, cyl])
+    bot_faces = shape.faces[Z < 0.1]
+    top_faces = shape.faces[Z > 0.1]
+    bot_faces.Identify(top_faces, "", idtype, gp_Trsf.Translation((0,0,1)))
+
+    geo = OCCGeometry(shape)
+    mesh = geo.GenerateMesh(maxh=0.3)
+
+    for el in mesh.Elements3D():
+        assert len(el.vertices)==6