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

CMakeLists.txt

@@ -375,28 +375,34 @@ if (USE_OCC)
       TKGeomAlgo
       TKGeomBase
       TKHLR
-      TKIGES
       TKLCAF
       TKMath
       TKMesh
       TKOffset
       TKPrim
-      TKSTEP
-      TKSTEP209
-      TKSTEPAttr
-      TKSTEPBase
-      TKSTL
       TKService
       TKShHealing
       TKTopAlgo
       TKV3d
       TKVCAF
       TKXCAF
-      TKXDEIGES
-      TKXDESTEP
       TKXSBase
       TKernel
     )
+    if(${OpenCASCADE_MAJOR_VERSION}.${OpenCASCADE_MINOR_VERSION} VERSION_GREATER_EQUAL 7.8)
+      list(APPEND OCC_LIBRARIES TKDEIGES TKDESTEP TKDESTL)
+    else()
+      list(APPEND OCC_LIBRARIES
+        TKIGES
+        TKSTEP
+        TKSTL
+        TKXDEIGES
+        TKXDESTEP
+        TKSTEP209
+        TKSTEPAttr
+        TKSTEPBase
+        )
+    endif()
     if(UNIX AND NOT APPLE)
       list(PREPEND OCC_LIBRARIES -Wl,--start-group)
       list(APPEND OCC_LIBRARIES -Wl,--end-group)

cmake/SuperBuild.cmake

@@ -91,6 +91,8 @@ if(BUILD_OCC)
   ExternalProject_Add(project_occ
     URL https://github.com/Open-Cascade-SAS/OCCT/archive/refs/tags/V7_6_3.zip
     URL_MD5 2426e373903faabbd4f96a01a934b66d
+    # URL https://github.com/Open-Cascade-SAS/OCCT/archive/refs/tags/V7_8_0.zip
+    # URL_MD5 f4432df8e42cb6178ea09a7448427f6c
     DOWNLOAD_DIR ${CMAKE_CURRENT_SOURCE_DIR}/external_dependencies
     ${SUBPROJECT_ARGS}
     CMAKE_ARGS

libsrc/core/CMakeLists.txt

@@ -53,7 +53,7 @@ if(USE_PYTHON)
 endif(USE_PYTHON)
 
 if(WIN32)
-  target_compile_options(ngcore PUBLIC /bigobj /MP /W1 /wd4068)
+  target_compile_options(ngcore PUBLIC /bigobj $<BUILD_INTERFACE:/MP;/W1;/wd4068>)
   get_WIN32_WINNT(ver)
   target_compile_definitions(ngcore PUBLIC _WIN32_WINNT=${ver} WNT WNT_WINDOW NOMINMAX MSVC_EXPRESS _CRT_SECURE_NO_WARNINGS HAVE_STRUCT_TIMESPEC WIN32)
   target_link_options(ngcore PUBLIC /ignore:4273 /ignore:4217 /ignore:4049)

libsrc/general/table.hpp

@@ -169,6 +169,7 @@ public:
   /// Inserts element acont into row i. BASE-based. Does not test if already used, assumes to have enough memory
   inline void AddSave (int i, const T & acont)
     {
+      NETGEN_CHECK_RANGE(i, BASE, data.Size()+BASE);
       ((T*)data[i-BASE].col)[data[i-BASE].size] = acont;
       data[i-BASE].size++;
     }

libsrc/meshing/basegeom.cpp

@@ -476,6 +476,20 @@ namespace netgen
     static Timer tdivide("Divide Edges");
     RegionTimer rt(tdivide);
     // -------------------- DivideEdge -----------------
+    if(properties.partition)
+      {
+        points.SetSize(properties.partition->Size());
+        params.SetSize(properties.partition->Size()+2);
+        params[0] = 0.0;
+        params.Last() = 1.0;
+        for(auto i : Range(properties.partition->Size()))
+          {
+            params[i+1] = (*properties.partition)[i];
+            points[i] = GetPoint(params[i+1]);
+          }
+        return;
+      }
+
     tdivedgesections.Start();
     auto layer = properties.layer;
     double safety = 0.5*(1.-mparam.grading);

libsrc/meshing/basegeom.hpp

@@ -27,10 +27,12 @@ namespace netgen
     double hpref = 0;  // number of hp refinement levels (will be multiplied by factor later)
     int layer = 1;
     optional<bool> quad_dominated;
+    optional<Array<double>> partition;
     void Merge(const ShapeProperties & prop2)
     {
       if (!name && prop2.name) name = prop2.name;
       if (!col && prop2.col) col = prop2.col;
+      if (!partition && prop2.partition) partition = prop2.partition;
       maxh = min2(maxh, prop2.maxh);
       hpref = max2(hpref, prop2.hpref);
       if(!quad_dominated.has_value()) quad_dominated = prop2.quad_dominated;

libsrc/meshing/meshclass.cpp

@@ -1,5 +1,6 @@
 #include <mystdlib.h>
 #include <atomic>
+#include <regex>
 #include <set>
 #include "meshing.hpp"
 #include "../general/gzstream.h"
@@ -318,27 +319,35 @@ namespace netgen
     hglob = mesh2.hglob;
     hmin = mesh2.hmin;
     maxhdomain = mesh2.maxhdomain;
+    pointelements = mesh2.pointelements;
 
+    // Remap string* values to new mesh
+    std::map<const string*, string*> names_map;
+    for (auto fi : Range(facedecoding))
+      names_map[&mesh2.facedecoding[fi].bcname] = &facedecoding[fi].bcname;
+
+    auto get_name = [&](const string *old_name) -> string* {
+      if (!old_name) return nullptr;
+      if (names_map.count(old_name)) return names_map[old_name];
+      return new string(*old_name);
+    };
 
     materials.SetSize( mesh2.materials.Size() );
     for ( int i = 0; i < mesh2.materials.Size(); i++ )
-      if ( mesh2.materials[i] ) materials[i] = new string ( *mesh2.materials[i] );
+    {
+      const string * old_name = mesh2.materials[i];
+      if ( old_name ) materials[i] = dimension == 2 ? get_name(old_name) : new string ( *old_name );
       else materials[i] = 0;
+    }
 
-    std::map<const string*, string*> bcmap;
     bcnames.SetSize( mesh2.bcnames.Size() );
     for ( int i = 0; i < mesh2.bcnames.Size(); i++ )
     {
-      if ( mesh2.bcnames[i] ) bcnames[i] = new string ( *mesh2.bcnames[i] );
+      const string * old_name = mesh2.bcnames[i];
+      if ( old_name ) bcnames[i] = dimension == 3 ? get_name(old_name) : new string ( *old_name );
       else bcnames[i] = 0;
-      bcmap[mesh2.bcnames[i]] = bcnames[i];
     }
 
-    // Remap string* members in FaceDescriptor to new mesh
-    for (auto & f : facedecoding)
-      f.SetBCName( bcmap[&f.GetBCName()] );
-
-
     cd2names.SetSize(mesh2.cd2names.Size());
     for (int i=0; i < mesh2.cd2names.Size(); i++)
       if (mesh2.cd2names[i]) cd2names[i] = new string(*mesh2.cd2names[i]);
@@ -7201,6 +7210,83 @@ namespace netgen
     UpdateTopology();
   }
 
+  shared_ptr<Mesh> Mesh :: GetSubMesh(string domains, string faces) const
+  {
+    // Copy the mesh into a new one, then delete unwanted elements
+    // Unused points are deleted by the Compress() function at the end
+    auto mesh_ptr = make_unique<Mesh>();
+    auto & mesh = *mesh_ptr;
+    mesh = (*this);
+
+    auto ndomains = GetNDomains();
+    auto nfaces = GetNFD();
+
+    BitArray keep_point(GetNP()+1);
+    BitArray keep_face(nfaces+1);
+    BitArray keep_domain(ndomains+1);
+    keep_point.Clear();
+    keep_face.Clear();
+    keep_domain.Clear();
+
+    regex regex_faces(faces);
+    regex regex_domains(domains);
+
+    for(auto dom : Range(ndomains))
+    {
+      if(regex_match(mesh.GetMaterial(dom), regex_domains))
+        keep_domain.SetBit(dom);
+    }
+
+    for(auto fi : Range(nfaces))
+    {
+      auto & fd = mesh.FaceDescriptors()[fi];
+      if (keep_domain[fd.DomainIn()] || keep_domain[fd.DomainOut()] || regex_match(fd.GetBCName(), regex_faces))
+        keep_face.SetBit(fd.BCProperty());
+    }
+
+    auto filter_elements = [&mesh, &keep_point](auto & elements, auto & keep_region)
+    {
+      for(auto & el : elements)
+      {
+        if(keep_region[el.GetIndex()])
+          for (auto pi : el.PNums())
+            keep_point.SetBit(pi);
+        else
+          el.Delete();
+      }
+    };
+
+    filter_elements(mesh.VolumeElements(), keep_domain);
+    filter_elements(mesh.SurfaceElements(), keep_face);
+
+    // Keep line segments only if all points are kept
+    // Check them in reverse order because they are deleted from the end
+    auto nsegments = mesh.LineSegments().Size();
+    for(auto i : Range(nsegments))
+    {
+      SegmentIndex segi = nsegments-i-1;
+      auto seg = mesh[segi];
+      bool keep = true;
+      for(auto pi : seg.PNums())
+        keep &= keep_point[pi];
+
+      if(!keep)
+        mesh.LineSegments().DeleteElement(segi);
+    }
+
+    // Check in reverse order because they are deleted from the end
+    auto npointelements = mesh.pointelements.Size();
+    for(auto i : Range(npointelements))
+    {
+      auto pel = mesh.pointelements[npointelements-i-1];
+      if(!keep_point[pel.pnum])
+        mesh.pointelements.DeleteElement(npointelements-i-1);
+    }
+
+    mesh.Compress();
+    return mesh_ptr;
+  }
+
   void Mesh :: SetMaterial (int domnr, const string & mat)
   {
     if (domnr > materials.Size())

libsrc/meshing/meshclass.hpp

@@ -700,6 +700,7 @@ namespace netgen
     void SetCommunicator(NgMPI_Comm acomm);
     
     DLL_HEADER void SplitFacesByAdjacentDomains();
+    DLL_HEADER shared_ptr<Mesh> GetSubMesh(string domains="", string faces="") const;
 
     ///
     DLL_HEADER void SetMaterial (int domnr, const string & mat);

libsrc/meshing/meshfunc.cpp

@@ -688,13 +688,16 @@ namespace netgen
 
     // optimize only bad elements first
     optmesh.SetMinBadness(1000.);
+    bool do_split = mp.optimize3d.find('d') != string::npos;
+    bool do_swap = mp.optimize3d.find('s') != string::npos;
+    bool do_swap2 = mp.optimize3d.find('t') != string::npos;
     for(auto i : Range(mp.optsteps3d))
       {
         auto [total_badness, max_badness, bad_els] = optmesh.UpdateBadness();
         if(bad_els==0) break;
-        optmesh.SplitImprove();
-        optmesh.SwapImprove();
-        optmesh.SwapImprove2();
+        if(do_split) optmesh.SplitImprove();
+        if(do_swap) optmesh.SwapImprove();
+        if(do_swap2) optmesh.SwapImprove2();
       }
 
     // Now optimize all elements

libsrc/meshing/python_mesh.cpp

@@ -1251,6 +1251,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     .def ("GetCD3Name", &Mesh::GetCD3Name)
     .def ("SetCD3Name", &Mesh::SetCD3Name)
     .def ("SplitFacesByAdjacentDomains", &Mesh::SplitFacesByAdjacentDomains)
+    .def ("GetSubMesh", &Mesh::GetSubMesh, py::arg("domains")="", py::arg("faces")="")
     .def("GetIdentifications", [](Mesh & self) -> py::list
          {
            py::list points;

libsrc/occ/Partition_Loop3d.hxx

@@ -29,7 +29,6 @@
 #if OCC_VERSION_HEX < 0x070000
 #else
    #include <TopTools_ShapeMapHasher.hxx>
-   #include <TopTools_OrientedShapeMapHasher.hxx>
    #include <TopTools_MapOfOrientedShape.hxx>
 #endif
 

libsrc/occ/python_occ_shapes.cpp

@@ -683,6 +683,10 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
     .export_values()
     ;
   
+  m.def("ResetGlobalShapeProperties", [] () {
+    OCCGeometry::global_shape_properties.clear();
+    OCCGeometry::global_shape_property_indices.Clear();
+  });
   
   py::class_<TopoDS_Shape> (m, "TopoDS_Shape")
     .def("__str__", [] (const TopoDS_Shape & shape)
@@ -1398,7 +1402,20 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
                              return tuple(s0, s1);
                            },
                            "parameter interval of curve")
-    
+    .def_property("partition",
+       [](TopoDS_Shape & self) -> optional<Array<double>>
+       {
+         if (OCCGeometry::HaveProperties(self))
+           return OCCGeometry::GetProperties(self).partition;
+         return nullopt;
+       },
+       [](TopoDS_Shape &self, py::array_t<double> val)
+       {
+         Array<double> partition(val.size());
+         for(auto i : Range(partition))
+           partition[i] = val.at(i);
+         OCCGeometry::GetProperties(self).partition = std::move(partition);
+       })
     
     .def("Split", [](const TopoDS_Edge& self, py::args args)
     {

nglib/CMakeLists.txt

@@ -6,7 +6,7 @@ if(USE_OCC)
 endif(USE_OCC)
 
 if(EMSCRIPTEN)
-  target_compile_options(nglib PUBLIC $<TARGET_PROPERTY:ngcore,INTERFACE_COMPILE_OPTIONS>)
+  target_compile_options(nglib PUBLIC $<BUILD_INTERFACE:$<TARGET_PROPERTY:ngcore,INTERFACE_COMPILE_OPTIONS>>)
   target_compile_definitions(nglib PUBLIC $<TARGET_PROPERTY:ngcore,INTERFACE_COMPILE_DEFINITIONS>)
   target_include_directories(nglib PUBLIC $<TARGET_PROPERTY:ngcore,INTERFACE_INCLUDE_DIRECTORIES>)
 else(EMSCRIPTEN)

python/CMakeLists.txt

@@ -42,8 +42,12 @@ else()
     message(WARNING "pybind11-stubgen version is too old, if you want to create stub files for better autocompletion support upgrade it with pip.")
   else()
   message("-- Found pybind11-stubgen version: ${stubgen_version}")
-  install(CODE "execute_process(COMMAND ${Python3_EXECUTABLE} -m pybind11_stubgen --ignore-all-errors netgen)")
-  install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/../stubs/netgen/ DESTINATION ${NG_INSTALL_DIR_PYTHON}/netgen/ COMPONENT netgen)
+  install(CODE "\
+  set(ENV{PYTHONPATH} ${CMAKE_INSTALL_PREFIX}/${NG_INSTALL_DIR_PYTHON})\n \
+  execute_process(COMMAND ${Python3_EXECUTABLE} -m pybind11_stubgen --ignore-all-errors netgen)\n \
+  execute_process(COMMAND ${Python3_EXECUTABLE} -m pybind11_stubgen --ignore-all-errors pyngcore)\n \
+  ")
+  install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/../stubs/ DESTINATION ${NG_INSTALL_DIR_PYTHON} COMPONENT netgen)
   endif()
 endif()
 endif(BUILD_STUB_FILES)

setup.py

@@ -8,7 +8,7 @@ from skbuild import setup
 import skbuild.cmaker
 from subprocess import check_output
 
-setup_requires = []
+setup_requires = ['pybind11-stubgen==2.5']
 
 def install_filter(cmake_manifest):
     print(cmake_manifest)
@@ -98,7 +98,7 @@ cmake_args += [
         '-DUSE_OCC=ON',
         '-DBUILD_FOR_CONDA=ON',
         f'-DNETGEN_PYTHON_PACKAGE_NAME={name}',
-        '-DBUILD_STUB_FILES=OFF',
+        '-DBUILD_STUB_FILES=ON',
 ]
 
 pyprefix = pathlib.Path(sys.prefix).as_posix()

tests/pytest/test_meshclass.py

@@ -1,16 +1,19 @@
 import pyngcore
 import netgen
+import pytest
+import tempfile
 
 from meshes import unit_mesh_3d
 
+
 def test_element_arrays(unit_mesh_3d):
     mesh = unit_mesh_3d
-    el0 = mesh.Elements0D()
     el1 = mesh.Elements1D()
     el2 = mesh.Elements2D()
     el3 = mesh.Elements3D()
     p = mesh.Points()
 
+    assert len(el1) > 0
     assert len(el2) > 0
     assert len(el3) > 0
     assert len(p) > 0
@@ -20,3 +23,43 @@ def test_element_arrays(unit_mesh_3d):
 
     for el in el3:
         assert len(el.vertices) == 4
+
+
+def test_copy_mesh():
+    pytest.importorskip("netgen.occ")
+    import netgen.occ as occ
+
+    box1 = occ.Box((0, 0, 0), (1, 1, 1))
+    box2 = occ.Box((1, 0, 0), (2, 1, 1))
+    box1.faces.name = "bnd1"
+    box1.name = "mat1"
+    box2.faces.name = "bnd2"
+    box2.name = "mat1"
+
+    geo = occ.OCCGeometry(occ.Glue([box1, box2]))
+    m3d = geo.GenerateMesh(maxh=99)
+
+    plane1 = occ.WorkPlane(occ.Axes((0, 0, 0), occ.X, occ.Y)).Rectangle(1, 1).Face()
+    plane1.name = "mat1"
+    plane1.edges.name = "bnd1"
+
+    plane2 = occ.WorkPlane(occ.Axes((0, 0, 0), occ.X, occ.Y)).Rectangle(2, 2).Face()
+    plane2.name = "mat2"
+    plane2.edges.name = "bnd2"
+
+    geo2 = occ.OCCGeometry(occ.Glue([plane1, plane2]), dim=2)
+    m2d = geo2.GenerateMesh(maxh=99)
+
+    for mesh in [m2d, m3d]:
+        copy = mesh.Copy()
+
+        assert copy.dim == mesh.dim
+        assert len(copy.Elements0D()) == len(mesh.Elements0D())
+        assert len(copy.Elements1D()) == len(mesh.Elements1D())
+        assert len(copy.Elements2D()) == len(mesh.Elements2D())
+        assert len(copy.Elements3D()) == len(mesh.Elements3D())
+        assert copy.GetNDomains() == mesh.GetNDomains()
+        assert copy.GetNFaceDescriptors() == mesh.GetNFaceDescriptors()
+        for dim in range(1, mesh.dim + 1):
+            assert copy.GetRegionNames(dim) == mesh.GetRegionNames(dim)
+        assert copy.GetIdentifications() == mesh.GetIdentifications()