Merge remote-tracking branch 'origin/master' into boundarylayer_fixes

CMakeLists.txt

@@ -87,12 +87,15 @@ set(NG_INSTALL_SUFFIX netgen CACHE STRING "Suffix appended to install directorie
 if(USE_PYTHON)
   if(${CMAKE_VERSION} VERSION_GREATER_EQUAL 3.18)
       find_package(Python3 REQUIRED COMPONENTS Development.Module)
+      if(NOT EMSCRIPTEN)
           find_package(Python3 COMPONENTS Interpreter Development.Embed)
+      endif()
   else()
     find_package(Python3 REQUIRED COMPONENTS Interpreter Development)
   endif()
 
   if(NOT CMAKE_CROSSCOMPILING)
+      find_package(Python3 REQUIRED COMPONENTS Interpreter)
       execute_process(COMMAND ${Python3_EXECUTABLE} -c "import os.path, sysconfig;print(os.path.relpath(sysconfig.get_path('platlib'), sysconfig.get_path('data')))" OUTPUT_VARIABLE PYTHON_PACKAGES_INSTALL_DIR OUTPUT_STRIP_TRAILING_WHITESPACE)
       file(TO_CMAKE_PATH ${PYTHON_PACKAGES_INSTALL_DIR} PYTHON_PACKAGES_INSTALL_DIR)
   endif(NOT CMAKE_CROSSCOMPILING)

cmake/SuperBuild.cmake

@@ -184,7 +184,9 @@ if (USE_PYTHON)
     endif( PYBIND_INCLUDE_DIR )
     if(${CMAKE_VERSION} VERSION_GREATER_EQUAL 3.18)
       find_package(Python3 COMPONENTS Interpreter Development.Module)
+      if(NOT EMSCRIPTEN)
           find_package(Python3 COMPONENTS Interpreter Development.Embed)
+      endif()
     else()
       find_package(Python3 REQUIRED COMPONENTS Interpreter Development)
     endif()

libsrc/meshing/basegeom.cpp

@@ -22,7 +22,7 @@ namespace netgen
           ArrayMem<int, 1> points;
           tree.GetIntersecting(p, p, points);
           if(points.Size()==0)
-              throw Exception("cannot find mapped point");
+              throw Exception("cannot find mapped point " + ToString(p));
           return points[0];
       }
 
@@ -313,23 +313,14 @@ namespace netgen
               {
                   bool need_inverse = ident.from == s.get();
                   auto other = need_inverse ? ident.to : ident.from;
-                  if(other->nr <= s->primary->nr)
+                  if(other->primary->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 = other->nr != s->primary->nr;
-                  }
-                  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 = other->primary->nr != s->primary->nr;
+                      changed = true;
                   }
               }
             }
@@ -568,12 +559,13 @@ namespace netgen
 
     auto & identifications = mesh.GetIdentifications();
 
-    std::map<size_t, PointIndex> vert2meshpt;
+    Array<PointIndex> vert2meshpt(vertices.Size());
+    vert2meshpt = PointIndex::INVALID;
     for(auto & vert : vertices)
       {
         auto pi = mesh.AddPoint(vert->GetPoint(), vert->properties.layer);
         tree.Insert(mesh[pi], pi);
-        vert2meshpt[vert->GetHash()] = pi;
+        vert2meshpt[vert->nr] = pi;
         mesh[pi].Singularity(vert->properties.hpref);
         mesh[pi].SetType(FIXEDPOINT);
 
@@ -585,8 +577,8 @@ namespace netgen
 
     for(auto & vert : vertices)
         for(auto & ident : vert->identifications)
-            identifications.Add(vert2meshpt[ident.from->GetHash()],
+            identifications.Add(vert2meshpt[ident.from->nr],
-                                vert2meshpt[ident.to->GetHash()],
+                                vert2meshpt[ident.to->nr],
                                 ident.name,
                                 ident.type);
 
@@ -600,8 +592,8 @@ namespace netgen
         auto edge = edges[edgenr].get();
         PointIndex startp, endp;
         // throws if points are not found
-        startp = vert2meshpt.at(edge->GetStartVertex().GetHash());
+        startp = vert2meshpt[edge->GetStartVertex().nr];
-        endp = vert2meshpt.at(edge->GetEndVertex().GetHash());
+        endp = vert2meshpt[edge->GetEndVertex().nr];
 
         // ignore collapsed edges
         if(startp == endp && edge->GetLength() < 1e-10 * bounding_box.Diam())
@@ -944,11 +936,12 @@ namespace netgen
     }
 
     bool have_identifications = false;
+    std::map<tuple<PointIndex, int>, PointIndex> mapto;
     for(auto & face : faces)
         if(face->primary != face.get())
         {
             have_identifications = true;
-            MapSurfaceMesh(mesh, *face);
+            MapSurfaceMesh(mesh, *face, mapto);
         }
 
     // identify points on faces
@@ -992,7 +985,8 @@ namespace netgen
                 if(ident.from == face.get())
                     for(auto pi : pi_of_face[face->nr])
                     {
-                        auto pi_other = tree.Find(ident.trafo(mesh[pi]));
+                        auto pi_primary = ident.from->primary->nr == ident.from->nr ? pi : mapto[{pi, ident.to->primary->nr}];
+                        auto pi_other = ident.to->primary->nr == ident.to->nr ? pi_primary : mapto[{pi_primary, ident.to->nr}];
                         mesh_ident.Add(pi, pi_other, ident.name, ident.type);
                     }
             }
@@ -1002,7 +996,7 @@ namespace netgen
     multithread.task = savetask;
   }
 
-  void NetgenGeometry :: MapSurfaceMesh( Mesh & mesh, const GeometryFace & dst ) const
+  void NetgenGeometry :: MapSurfaceMesh( Mesh & mesh, const GeometryFace & dst, std::map<tuple<PointIndex, int>, PointIndex> & mapto ) const
   {
     static Timer timer("MapSurfaceMesh");
     RegionTimer rt(timer);
@@ -1083,6 +1077,8 @@ namespace netgen
                 pmap[pi] = mesh.AddPoint(trafo(mesh[pi]), 1, SURFACEPOINT);
               }
               sel_new[i] = pmap[pi];
+              mapto[{pi, dst.nr}] = pmap[pi];
+              mapto[{pmap[pi], src.nr}] = pi;
           }
           if(do_invert.IsTrue())
               sel_new.Invert();

libsrc/meshing/basegeom.hpp

@@ -68,7 +68,6 @@ namespace netgen
     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;
   };
 
@@ -243,7 +242,7 @@ namespace netgen
     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 MapSurfaceMesh( Mesh & mesh, const GeometryFace & dst, std::map<tuple<PointIndex, int>, PointIndex> & mapto) const;
     virtual void OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const;
 
     virtual void FinalizeMesh(Mesh& mesh) const;
@@ -320,13 +319,6 @@ namespace netgen
       throw Exception("Base geometry get tangent called");
     }
 
-    virtual size_t GetEdgeIndex(const GeometryEdge& edge) const
-    {
-      for(auto i : Range(edges))
-        if(edge.GetHash() == edges[i]->GetHash())
-          return i;
-      throw Exception("Couldn't find edge index");
-    }
     virtual void Save (const filesystem::path & filename) const;
     virtual void SaveToMeshFile (ostream & /* ost */) const { ; }
   };

libsrc/meshing/improve3.hpp

@@ -12,8 +12,8 @@ extern double CalcTotalBad (const Mesh::T_POINTS & points,
 ///
 class MeshOptimize3d
 {
-  const MeshingParameters & mp;
   Mesh & mesh;
+  const MeshingParameters & mp;
   OPTIMIZEGOAL goal = OPT_QUALITY;
   double min_badness = 0;
 

libsrc/occ/occ_edge.cpp

@@ -53,11 +53,6 @@ namespace netgen
         throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__));
     }
 
-    size_t OCCEdge::GetHash() const
-    {
-      return edge.HashCode(std::numeric_limits<Standard_Integer>::max());
-    }
-
     void OCCEdge::ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const
     {
         auto pnt = ng2occ(p);

libsrc/occ/occ_edge.hpp

@@ -36,7 +36,6 @@ namespace netgen
         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;
         bool IsDegenerated(double) const override {

libsrc/occ/occ_face.cpp

@@ -30,11 +30,6 @@ namespace netgen
         return 0;
     }
 
-    size_t OCCFace::GetHash() const
-    {
-      return face.HashCode(std::numeric_limits<Standard_Integer>::max());
-    }
-
     Point<3> OCCFace::GetCenter() const
     {
         return occ2ng( props.CentreOfMass() );

libsrc/occ/occ_face.hpp

@@ -31,7 +31,6 @@ namespace netgen
 
         const TopoDS_Face Shape() const { return face; }
 
-        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;

libsrc/occ/occ_solid.hpp

@@ -16,8 +16,6 @@ namespace netgen
         OCCSolid(TopoDS_Shape dshape)
             : solid(TopoDS::Solid(dshape))
         { }
-
-        size_t GetHash() const override { return solid.HashCode(std::numeric_limits<Standard_Integer>::max()); }
     };
 }
 

libsrc/occ/occ_vertex.cpp

@@ -16,9 +16,4 @@ namespace netgen
     {
         return p;
     }
-
-    size_t OCCVertex::GetHash() const
-    {
-      return vertex.HashCode(std::numeric_limits<Standard_Integer>::max());
-    }
 }

libsrc/occ/occ_vertex.hpp

@@ -24,7 +24,6 @@ namespace netgen
         OCCVertex( TopoDS_Shape s );
         ~OCCVertex() {}
         Point<3> GetPoint() const override;
-        size_t GetHash() const override;
     };
 }
 

libsrc/occ/occgeom.cpp

@@ -1710,13 +1710,6 @@ namespace netgen
         BRepTools::Read(shape, ss, builder);
       }
 
-    /*
-    // enumerate shapes and archive only integers
-    auto my_hash = [](const TopoDS_Shape & key) {
-        auto occ_hash = key.HashCode(1<<31UL);
-        return std::hash<decltype(occ_hash)>()(occ_hash);
-    };
-    */
     TopTools_IndexedMapOfShape shape_map;
     Array<TopoDS_Shape> shape_list;
 
@@ -1784,7 +1777,7 @@ namespace netgen
    "Face", "Wire", "Edge", "Vertex"};
 
   const char * orientationstring[] =
-     {"+", "-"};
+    {"+", "-", "i", "e"};
 
 
 
@@ -1843,6 +1836,16 @@ namespace netgen
          }
 
          str << "{" << shapename[l] << " " << count2;
+         if(HaveProperties(e.Current()))
+           {
+             const auto& props = GetProperties(e.Current());
+             if(props.name || props.maxh < 1e99)
+               str << " - ";
+             if(props.name)
+               str << props.GetName();
+             if(props.maxh < 1e99)
+               str << " maxh(" << props.maxh << ")";
+           }
 
          if (l <= TopAbs_EDGE)
          {
@@ -2338,20 +2341,24 @@ namespace netgen
 
           for(auto & ident : identifications)
           {
+              const auto& to = STEPConstruct::FindEntity(finder, ident.from == shape ? ident.to : ident.from);
+              if(to.IsNull())
+                  continue;
               Array<Handle(StepRepr_RepresentationItem)> items;
-              // items.Append(STEPConstruct::FindEntity(finder, ident.other)); // TODO!
+              items.Append(MakeReal(ident.from == shape ? 1 : 0));
+              items.Append(to);
               auto & m = ident.trafo.GetMatrix();
               for(auto i : Range(9))
                   items.Append(MakeReal(m(i)));
               auto & v = ident.trafo.GetVector();
               for(auto i : Range(3))
                   items.Append(MakeReal(v(i)));
-              for(auto & item : items.Range(1,items.Size()))
+              items.Append(MakeInt(ident.type));
+              for(auto & item : items.Range(0, items.Size()))
                   model->AddEntity(item);
               ident_items.Append(MakeCompound(items, ident.name));
           }
-
+          for(auto & item : ident_items.Range(1, ident_items.Size()))
-          for(auto & item : ident_items.Range(1,ident_items.Size()))
             model->AddEntity(item);
           auto comp = MakeCompound(ident_items, "netgen_geometry_identification");
           model->AddEntity(comp);
@@ -2369,7 +2376,18 @@ 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))); /TODO!
+              auto is_from = ReadReal(id_item->ItemElementValue(1));
+              if(is_from)
+                {
+                  ident.from = shape_origin;
+                  ident.to = TransferBRep::ShapeResult(transProc->Find(id_item->ItemElementValue(2)));
+                }
+              else
+                {
+                  ident.from = TransferBRep::ShapeResult(
+                      transProc->Find(id_item->ItemElementValue(2)));
+                  ident.to = shape_origin;
+                }
 
               auto & m = ident.trafo.GetMatrix();
               for(auto i : Range(9))
@@ -2377,7 +2395,7 @@ namespace netgen
               auto & v = ident.trafo.GetVector();
               for(auto i : Range(3))
                   v(i) = ReadReal(id_item->ItemElementValue(12+i));
-
+              ident.type = Identifications::ID_TYPE(ReadInt(id_item->ItemElementValue(15)));
               result.push_back(ident);
           }
           OCCGeometry::GetIdentifications(shape_origin) = result;

python/__init__.py

@@ -41,6 +41,7 @@ if sys.platform.startswith('win'):
 del sys
 del os
 
+from pyngcore import Timer
 from . import libngpy
 
 from netgen.libngpy._meshing import _Redraw
@@ -48,7 +49,6 @@ from netgen.libngpy._meshing import _Redraw
 def Redraw(*args, **kwargs):
     return _Redraw(*args, **kwargs)
 
-from pyngcore import Timer
 def TimeFunction(func, name=None):
     name = name or func.__qualname__
     timer = Timer(name)