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

.gitlab-ci.yml

@@ -32,13 +32,23 @@ push_github:
     - "echo off"
     - call "C:\Program Files (x86)\Microsoft Visual Studio\2019\Community\VC\Auxiliary\Build\vcvars64"
     - set CI_DIR=C:\ci\%CI_PIPELINE_ID%
-    - set CLCACHE_BASEDIR=C:\ci\%CI_PIPELINE_ID%
+    - set CCACHE_BASEDIR=C:\ci\%CI_PIPELINE_ID%
     - set NETGEN_BUILD_DIR=%CI_DIR%\build
     - set INSTALL_DIR=%CI_DIR%\install
     - set SRC_DIR=%CI_DIR%\src
     - set NETGENDIR=%INSTALL_DIR%\bin
     - set PYTHONPATH=%INSTALL_DIR%\lib\site-packages
-    - set PATH=%NETGENDIR%;%PATH%
+    - echo %PATH%
+    - set PATH=%INSTALL_DIR%\bin;C:\python312;C:\python312\bin;C:\python312\Scripts;C:\tools\;%PATH%
+    - echo %PATH%
+    - set CCACHE_HARDLINK=1
+    - set CCACHE_NOHASHDIR=1
+    - C:\tools\ccache -s
+    - C:\tools\ccache -M 20G
+    - dir C:\python312
+    - python.exe --version
+    - python.exe -m pip install -U netgen-occt netgen-occt-devel
+    - cmake --version
 
 build_win:
   <<: *win
@@ -54,12 +64,14 @@ build_win:
     - >-
         cmake %SRC_DIR%
         -G Ninja
+        -DCMAKE_PREFIX=C:/python312
+        -DPython3_ROOT_DIR=C:/python312
         -DCMAKE_INSTALL_PREFIX=%INSTALL_DIR%
         -DCHECK_RANGE=ON
         -DUSE_CGNS=ON
         -DUSE_OCC=ON
         -DUSE_CCACHE=ON
-        -DENABLE_UNIT_TESTS=ON
+        -DENABLE_UNIT_TESTS=OFF
         -DCMAKE_BUILD_TYPE=Release
     - cmake --build . --target install --config Release
 
@@ -69,7 +81,7 @@ test_win:
   script:
     - pip install pytest-check
     - cd tests\pytest
-    - python test_tutorials.py new_results.json
+    - REM python test_tutorials.py new_results.json
     - cd %NETGEN_BUILD_DIR%\netgen
     - ctest -C Release -V --output-on-failure
     - cd ..
@@ -183,15 +195,6 @@ test_build_ngsolve:
       netgen_${CI_PIPELINE_ID}_installed:${UBUNTU_VERSION}
       bash -c 'cd /root/src/netgen/tests/ && ./build_ngsolve.sh'
 
-# cpp guideline checks
-# test_guidelines:
-#   <<: *ubuntu
-#   stage: test
-#   script:
-#     - docker run -e CCACHE_DIR=/ccache -v /mnt/ccache:/ccache netgen_${CI_PIPELINE_ID}:${UBUNTU_VERSION} bash /root/src/netgen/tests/build_guidelines.sh
-#   when: always
-#   allow_failure: true
-
 cleanup_ubuntu:
   stage: cleanup
   tags:
@@ -285,12 +288,11 @@ pip_windows:
     - pip
     - windows
   script:
-    - .\tests\build_pip.ps1 C:\Python37
-    - .\tests\build_pip.ps1 C:\Python38
-    - .\tests\build_pip.ps1 C:\Python39
-    - .\tests\build_pip.ps1 C:\Python310
-    - .\tests\build_pip.ps1 C:\Python311
     - .\tests\build_pip.ps1 C:\Python312
+    - .\tests\build_pip.ps1 C:\Python311
+    - .\tests\build_pip.ps1 C:\Python310
+    - .\tests\build_pip.ps1 C:\Python39
+    - .\tests\build_pip.ps1 C:\Python38
   when: manual
 
 pip_macos:
@@ -300,9 +302,9 @@ pip_macos:
     - macosx
     - m1
   script:
-    - ./tests/build_pip_mac.sh 3.8
-    - ./tests/build_pip_mac.sh 3.9
-    - ./tests/build_pip_mac.sh 3.10
-    - ./tests/build_pip_mac.sh 3.11
     - ./tests/build_pip_mac.sh 3.12
+    - ./tests/build_pip_mac.sh 3.11
+    - ./tests/build_pip_mac.sh 3.10
+    - ./tests/build_pip_mac.sh 3.9
+    - ./tests/build_pip_mac.sh 3.8
   when: manual

CMakeLists.txt

@@ -164,6 +164,7 @@ if(USE_CCACHE)
   find_program(CCACHE_FOUND NAMES ccache ccache.bat)
   if(CCACHE_FOUND)
     set_property(GLOBAL PROPERTY RULE_LAUNCH_COMPILE ${CCACHE_FOUND})
+    message(STATUS "Using ccache ${CCACHE_FOUND}")
   endif(CCACHE_FOUND)
 endif(USE_CCACHE)
 
@@ -398,7 +399,12 @@ if (USE_OCC)
     endif()
 
     target_link_libraries(occ_libs INTERFACE ${OCC_LIBRARIES})
-    include_directories(${OpenCASCADE_INCLUDE_DIR})
+    get_target_property(occ_include_dir TKernel INTERFACE_INCLUDE_DIRECTORIES)
+    if(NOT occ_include_dir)
+      set(occ_include_dir ${OpenCASCADE_INCLUDE_DIR})
+    endif()
+    target_include_directories(occ_libs INTERFACE ${occ_include_dir})
+    message(STATUS "OpenCasCade include dirs: ${occ_include_dir}")
     if(NOT OpenCASCADE_BUILD_SHARED_LIBS)
       if(OpenCASCADE_WITH_FREETYPE)
         find_library( FREETYPE NAMES freetype HINTS ${OpenCASCADE_LIBRARY_DIR})
@@ -413,10 +419,12 @@ if (USE_OCC)
       find_package(Threads REQUIRED)
       target_link_libraries(occ_libs INTERFACE Threads::Threads)
     endif()
-    message(STATUS "OCC DIRS ${OpenCASCADE_INCLUDE_DIR}")
     if(WIN32 AND USE_GUI)
-      target_link_libraries(nggui PRIVATE occ_libs Ws2_32.lib)
+      target_link_libraries(nggui PRIVATE Ws2_32.lib)
     endif(WIN32 AND USE_GUI)
+    if(USE_GUI)
+      target_link_libraries(nggui PRIVATE occ_libs)
+    endif(USE_GUI)
 endif (USE_OCC)
 
 #######################################################################
@@ -625,6 +633,16 @@ if(UNIX)
     endif(temp)
 endif(UNIX)
 
+if(USE_PYTHON AND NOT SKBUILD)
+  # install egg file to let python/pip know that Netgen ist installed
+  file( WRITE "netgen_mesher-py3.egg-info"
+"Metadata-Version: 2.1
+Name: netgen-mesher
+Version: ${NETGEN_VERSION_MAJOR}.${NETGEN_VERSION_MINOR}.${NETGEN_VERSION_PATCH}.post${NETGEN_VERSION_TWEAK}
+")
+  install(FILES netgen_mesher-py3.egg-info DESTINATION ${NG_INSTALL_DIR_PYTHON} COMPONENT netgen)
+endif()
+
 if(APPLE AND NOT SKBUILD)
     # create some auxiliary files
     set(mac_startup ${CMAKE_CURRENT_BINARY_DIR}/startup.sh)

README.md

@@ -4,3 +4,4 @@ NETGEN is an automatic 3d tetrahedral mesh generator. It accepts input from cons
 
 Find the Open Source Community on https://ngsolve.org
 Support & Services: https://cerbsim.com
+

cmake/SuperBuild.cmake

@@ -120,28 +120,15 @@ if(BUILD_OCC)
   list(APPEND NETGEN_DEPENDENCIES project_occ)
   set(OpenCascade_ROOT ${OCC_DIR})
 else(BUILD_OCC)
-    if(WIN32 AND NOT OCC_INCLUDE_DIR AND NOT OpenCASCADE_DIR)
-        # we can download prebuilt occ binaries for windows
-        ExternalProject_Add(win_download_occ
-          ${SUBPROJECT_ARGS}
-          URL ${OCC_DOWNLOAD_URL_WIN}
-          UPDATE_COMMAND "" # Disable update
-          BUILD_IN_SOURCE 1
-          CONFIGURE_COMMAND ""
-          BUILD_COMMAND ""
-          INSTALL_COMMAND ${CMAKE_COMMAND} -E copy_directory . ${CMAKE_INSTALL_PREFIX}
-          )
-        list(APPEND NETGEN_DEPENDENCIES win_download_occ)
-    else()
     find_package(OpenCascade NAMES OpenCasCade OpenCASCADE opencascade)
     if(NOT OpenCascade_FOUND)
       message(FATAL_ERROR "Opencascade not found, either\n\
+      - install pip packages netgen-occt-devel netgen-occ\n\
       - set OpenCascade_DIR to a directory containting opencascadeConfig.cmake\n\
       - build OpenCascade automatically by passing -DBUILD_OCC=ON\n\
       - disable OpenCascade by passing -DUSE_OCC=OFF\n\
       ")
     endif()
-    endif()
 endif(BUILD_OCC)
 endif(USE_OCC)
 
@@ -164,9 +151,11 @@ if(BUILD_ZLIB)
     # force linking the static library
     set(ZLIB_INCLUDE_DIRS ${ZLIB_ROOT}/include)
     set(ZLIB_LIBRARIES ${ZLIB_ROOT}/lib/zlibstatic.lib)
-  elseif(EMSCRIPTEN)
+    set(ZLIB_LIBRARY_RELEASE ${ZLIB_ROOT}/lib/zlibstatic.lib)
+  elseif(WIN32)
     set(ZLIB_INCLUDE_DIRS ${ZLIB_ROOT}/include)
     set(ZLIB_LIBRARIES ${ZLIB_ROOT}/lib/libz.a)
+    set(ZLIB_LIBRARY_RELEASE ${ZLIB_ROOT}/lib/libz.a)
   endif(WIN32)
 else()
     include(cmake/external_projects/zlib.cmake)
@@ -270,6 +259,7 @@ set_vars( NETGEN_CMAKE_ARGS
   OpenCascade_ROOT
   ZLIB_INCLUDE_DIRS
   ZLIB_LIBRARIES
+  ZLIB_LIBRARY_RELEASE
   ZLIB_ROOT
 
   NGLIB_LIBRARY_TYPE

external_dependencies/pybind11

@@ -1 +1 @@
-Subproject commit 80dc998efced8ceb2be59756668a7e90e8bef917
+Subproject commit 38bf7b174875c27c1ba98bdf5a9bf13d967f14d4

libsrc/core/CMakeLists.txt

@@ -93,7 +93,7 @@ install(FILES ngcore.hpp archive.hpp type_traits.hpp version.hpp ngcore_api.hpp
   xbool.hpp signal.hpp bitarray.hpp table.hpp hashtable.hpp ranges.hpp ngstream.hpp
   simd.hpp simd_avx.hpp simd_avx512.hpp simd_generic.hpp simd_sse.hpp simd_arm64.hpp
   register_archive.hpp autodiff.hpp autodiffdiff.hpp
-  ng_mpi.hpp ng_mpi_generated_declarations.hpp mpi4py_pycapi.h
+  ng_mpi.hpp ng_mpi_generated_declarations.hpp mpi4py_pycapi.h ng_mpi_native.hpp
   DESTINATION ${NG_INSTALL_DIR_INCLUDE}/core COMPONENT netgen_devel)
 
 if(ENABLE_CPP_CORE_GUIDELINES_CHECK)

libsrc/core/archive.hpp

@@ -77,7 +77,8 @@ namespace ngcore
   {
     template <class T, class Tuple, size_t... Is>
     T* construct_from_tuple(Tuple&& tuple, std::index_sequence<Is...> ) {
-      return new T{std::get<Is>(std::forward<Tuple>(tuple))...};
+      // return new T{std::get<Is>(std::forward<Tuple>(tuple))...};
+      return new T{std::get<Is>(std::move(tuple))...};
     }
 
     template <class T, class Tuple>

libsrc/core/array.hpp

@@ -733,7 +733,7 @@ namespace ngcore
 
     NETGEN_INLINE Array (Array && a2) 
     {
-      mt.Swap(sizeof(T) * allocsize, a2.mt, sizeof(T) * a2.allocsize);
+      mt.Swap(0., a2.mt, sizeof(T) * a2.allocsize);
 
       size = a2.size; 
       data = a2.data;
@@ -1528,6 +1528,8 @@ namespace ngcore
   }
                                
 
+  struct HTAHelp { };
+  
   // head-tail array
   template <size_t S, typename T>
   class HTArray
@@ -1535,10 +1537,22 @@ namespace ngcore
     HTArray<S-1,T> tail;
     T head;
   public:
-    HTArray () = default;
+    constexpr HTArray () = default;
-    HTArray (const HTArray &) = default;
+    constexpr HTArray (const HTArray &) = default;
     template <typename T2>
-    HTArray (const HTArray<S,T2> & a2) : tail(a2.Tail()), head(a2.Head()) { ; }
+    constexpr HTArray (const HTArray<S,T2> & a2) : tail(a2.Tail()), head(a2.Head()) { ; }
+
+    constexpr HTArray (T v) : tail(v), head(v) { } // all the same
+    
+    template <class... T2,
+              std::enable_if_t<S==1+sizeof...(T2),bool> = true>
+    constexpr HTArray (const T &v, T2... rest)
+      : tail{HTAHelp(), v,rest...}, head(std::get<S-2>(std::tuple(rest...))) { }
+
+    template <class... T2>
+    constexpr HTArray (HTAHelp h, const T &v, T2... rest)
+      : tail{h, v,rest...}, head(std::get<S-2>(std::tuple(rest...))) { }
+
     
     HTArray & operator= (const HTArray &) = default;
 
@@ -1559,10 +1573,15 @@ namespace ngcore
   {
     T head;
   public:
-    HTArray () = default;
+    constexpr HTArray () = default;
-    HTArray (const HTArray &) = default;
+    constexpr HTArray (const HTArray &) = default;
     template <typename T2>
-    HTArray (const HTArray<1,T2> & a2) : head(a2.Head()) { ; }
+    constexpr HTArray (const HTArray<1,T2> & a2) : head(a2.Head()) { ; }
+    constexpr HTArray (T v) : head(v) { } // all the same
+    template <class... T2>    
+    constexpr HTArray (HTAHelp h, const T &v, T2... rest)
+      : head(v) { } 
+
     
     HTArray & operator= (const HTArray &) = default;
 
@@ -1590,7 +1609,7 @@ namespace ngcore
     HTArray (const HTArray &) = default;
     template <typename T2>
     HTArray (const HTArray<0,T2> & a2) { ; }
-    
+    constexpr HTArray (T v) { } // all the same    
     HTArray & operator= (const HTArray &) = default;
 
     /*

libsrc/core/exception.cpp

@@ -23,6 +23,47 @@ namespace ngcore
 }
 
 
+  Exception :: Exception(std::string_view s1, std::string_view s2)
+    : Exception(std::string(s1)+std::string(s2))
+  { }
+
+  Exception :: Exception(std::string_view s1, std::string_view s2, std::string_view s3)
+    : Exception(std::string(s1)+std::string(s2)+std::string(s3))
+  { }
+
+
+  void Exception :: Throw (std::string_view s1)
+  {
+    throw Exception(std::string(s1));
+  }
+
+  void Exception :: Throw (std::string_view s1, std::string_view s2)
+  {
+    throw Exception(std::string(s1)+std::string(s2));
+  }
+    
+  void Exception :: Throw (std::string_view s1, std::string_view s2, std::string_view s3)
+  {
+    throw Exception(std::string(s1)+std::string(s2)+std::string(s3));
+  }
+  
+  
+  RangeException :: RangeException (// const std::string & where,
+                                    const char * where,
+                                    int ind, int imin, int imax) : Exception("")
+  {
+    std::stringstream str;
+    str << where << ": index " << ind << " out of range [" << imin << "," << imax << ")\n";
+    Append (str.str());
+    Append (GetBackTrace());
+  }
+  
+
+  void ThrowRangeException(const char * s, int ind, int imin, int imax)
+  {
+    throw RangeException(s, ind, imin, imax);
+  }
+  
   void ThrowException(const std::string & s)
   {
     throw Exception (s);
@@ -32,6 +73,13 @@ namespace ngcore
   {
     throw Exception (s);
   }
+
+
+  void ThrowNotTheSameException(const char * s, long int a, long int b)
+  {
+    throw ngcore::Exception(std::string(s) + ", a="+ToString(a) + ", b="+ToString(b) + GetBackTrace());     
+  }
+  
 } // namespace ngcore
 
 

libsrc/core/exception.hpp

@@ -32,8 +32,14 @@ namespace ngcore
     Exception(Exception&&) = default;
     Exception(const std::string& s); //  : m_what(s) {}
     Exception(const char* s); //  : m_what(s) {}
+    Exception(std::string_view s1, std::string_view s2);
+    Exception(std::string_view s1, std::string_view s2, std::string_view s3);
     ~Exception() override = default;
 
+    [[noreturn]] static void Throw (std::string_view s1);
+    [[noreturn]] static void Throw (std::string_view s1, std::string_view s2);
+    [[noreturn]] static void Throw (std::string_view s1, std::string_view s2, std::string_view s3);
+    
     Exception& operator =(const Exception&) = default;
     Exception& operator =(Exception&&) noexcept = default;
 
@@ -49,23 +55,26 @@ namespace ngcore
     const char* what() const noexcept override { return m_what.c_str(); }
   };
   
-  NGCORE_API void ThrowException(const std::string & s);
+  [[noreturn]] NGCORE_API void ThrowException(const std::string & s);
-  NGCORE_API void ThrowException(const char * s);
+  [[noreturn]] NGCORE_API void ThrowException(const char * s);
   
   // Out of Range exception
   class NGCORE_API RangeException : public Exception
   {
   public:
     /// where it occurs, index, minimal and maximal indices
-    RangeException (const std::string & where,
+    RangeException (// const std::string & where,
-                    int ind, int imin, int imax) : Exception("")
+                    const char * where,
+                    int ind, int imin, int imax);
+    /*
+    : Exception("")
       {
         std::stringstream str;
         str << where << ": index " << ind << " out of range [" << imin << "," << imax << ")\n";
         Append (str.str());
         Append (GetBackTrace());
       }
-
+    */
     template<typename T>
     RangeException(const std::string& where, const T& value)
     {
@@ -75,6 +84,10 @@ namespace ngcore
     }
   };
 
+  [[noreturn]] NGCORE_API void ThrowRangeException(const char * s, int ind, int imin, int imax);
+  [[noreturn]] NGCORE_API void ThrowNotTheSameException(const char * s, long int a, long int b);
+  
+  
   // Exception used if no simd implementation is available to fall back to standard evaluation
   class NGCORE_API ExceptionNOSIMD : public Exception
   { public: using Exception::Exception; };
@@ -86,21 +99,29 @@ namespace ngcore
 // Convenience macro to append file name and line of exception origin to the string
 #define NG_EXCEPTION(s) ngcore::Exception(__FILE__ ":" NETGEN_CORE_NGEXEPTION_STR(__LINE__) "\t"+std::string(s))
 
+
+template <typename T>
+struct IsSafe {
+  constexpr operator bool() const { return false; } };
+
 #if defined(NETGEN_ENABLE_CHECK_RANGE) && !defined(__CUDA_ARCH__)
 #define NETGEN_CHECK_RANGE(value, min, max_plus_one) \
-  { if ((value)<(min) ||  (value)>=(max_plus_one)) \
+  { if constexpr (!IsSafe<decltype(value)>()) { \
-      throw ngcore::RangeException(__FILE__ ":" NETGEN_CORE_NGEXEPTION_STR(__LINE__) "\t", int(value), int(min), int(max_plus_one)); }
+  if ((value)<(min) ||  (value)>=(max_plus_one))                        \
-#define NETGEN_CHECK_SHAPE(a,b) \
+    ThrowRangeException(__FILE__ ":" NETGEN_CORE_NGEXEPTION_STR(__LINE__) "\t", int(value), int(min), int(max_plus_one)); }  }
-  { if(a.Shape() != b.Shape()) \
+
-      throw ngcore::Exception(__FILE__ ":" NETGEN_CORE_NGEXEPTION_STR(__LINE__) "\t: shape don't match"); }
 #define NETGEN_CHECK_SAME(a,b) \
-  { if(a != b) \
+  { if(a != b) {                                                        \
-      throw ngcore::Exception(__FILE__ ":" NETGEN_CORE_NGEXEPTION_STR(__LINE__) "\t: not the same, a="+ToString(a) + ", b="+ToString(b) + GetBackTrace()); }
+      if constexpr(std::is_same<decltype(a),size_t>() && std::is_same<decltype(b),size_t>()) \
+        ThrowNotTheSameException(__FILE__ ":" NETGEN_CORE_NGEXEPTION_STR(__LINE__) "\t: not the same, a=", long(a), long(b)); \
+      else \
+        throw ngcore::Exception(__FILE__ ":" NETGEN_CORE_NGEXEPTION_STR(__LINE__) "\t: not the same, a="+ToString(a) + ", b="+ToString(b) + GetBackTrace()); \
+    } }
 #define NETGEN_NOEXCEPT 
 #else // defined(NETGEN_ENABLE_CHECK_RANGE) && !defined(__CUDA_ARCH__)
 #define NETGEN_CHECK_RANGE(value, min, max)
 #define NETGEN_CHECK_SAME(a,b)
-#define NETGEN_CHECK_SHAPE(a,b)
+// #define NETGEN_CHECK_SHAPE(a,b)
 #define NETGEN_NOEXCEPT noexcept
 #endif // defined(NETGEN_ENABLE_CHECK_RANGE) && !defined(__CUDA_ARCH__)
 

libsrc/core/flags.cpp

@@ -16,6 +16,7 @@
 namespace ngcore
 {
   using std::string;
+  using std::string_view;
   using std::endl;
   Flags :: Flags () { ; }
 
@@ -209,18 +210,18 @@ namespace ngcore
   }
 
 
-  double Flags :: GetNumFlag (const string & name, double def) const
+  double Flags :: GetNumFlag (string_view name, double def) const
   {
     if (numflags.Used (name))
-      return numflags[name];
+      return numflags[string(name)];
     else
       return def;
   }
   
-  const double * Flags :: GetNumFlagPtr (const string & name) const
+  const double * Flags :: GetNumFlagPtr (string_view name) const
   {
     if (numflags.Used (name))
-      return & ((SymbolTable<double>&)numflags)[name];
+      return & ((SymbolTable<double>&)numflags)[string(name)];
     else
       return NULL;
   }
@@ -239,16 +240,16 @@ namespace ngcore
     return defflags.Used (name);
   }
   */
-  bool Flags :: GetDefineFlag (const string & name) const throw()
+  bool Flags :: GetDefineFlag (string_view name) const throw()
   {
-    if (!defflags.Used (name)) return false;
+    if (!defflags.Used (string(name))) return false;
-    return defflags[name];
+    return defflags[string(name)];
   }
 
-  xbool Flags :: GetDefineFlagX (const string & name) const throw()
+  xbool Flags :: GetDefineFlagX (string_view name) const throw()
   {
-    if (!defflags.Used (name)) return maybe;
+    if (!defflags.Used (string(name))) return maybe;
-    return bool(defflags[name]);
+    return bool(defflags[string(name)]);
   }
 
 
@@ -296,32 +297,32 @@ namespace ngcore
     return empty;
   }
 
-  bool Flags :: StringFlagDefined (const string & name) const
+  bool Flags :: StringFlagDefined (string_view name) const noexcept
   {
     return strflags.Used (name);
   }
 
-  bool Flags :: NumFlagDefined (const string &name) const
+  bool Flags :: NumFlagDefined (string_view name) const noexcept
   {
     return numflags.Used (name);
   }
 
-  bool Flags :: FlagsFlagDefined (const string &name) const
+  bool Flags :: FlagsFlagDefined (string_view name) const noexcept
   {
     return flaglistflags.Used (name);
   }
   
-  bool Flags :: StringListFlagDefined (const string & name) const
+  bool Flags :: StringListFlagDefined (string_view name) const noexcept
   {
     return strlistflags.Used (name);
   }
 
-  bool Flags :: NumListFlagDefined (const string & name) const
+  bool Flags :: NumListFlagDefined (string_view name) const noexcept
   {
     return numlistflags.Used (name);
   }
 
-  bool Flags :: AnyFlagDefined (const string& name) const
+  bool Flags :: AnyFlagDefined (string_view name) const noexcept
   {
     return anyflags.Used(name);
   }

libsrc/core/flags.hpp

@@ -125,15 +125,15 @@ namespace ngcore
     /// Returns std::string flag, default value if not exists
     std::string GetStringFlag (const std::string & name, std::string def = "") const;
     /// Returns numerical flag, default value if not exists
-    double GetNumFlag (const std::string & name, double def) const;
+    double GetNumFlag (std::string_view name, double def) const;
     /// Returns address of numerical flag, null if not exists
-    const double * GetNumFlagPtr (const std::string & name) const;
+    const double * GetNumFlagPtr (std::string_view name) const;
     /// Returns address of numerical flag, null if not exists
     double * GetNumFlagPtr (const std::string & name);
     /// Returns boolean flag
     // int GetDefineFlag (const char * name) const;
-    bool GetDefineFlag (const std::string & name) const throw();
+    bool GetDefineFlag (std::string_view name) const  noexcept;
-    xbool GetDefineFlagX (const std::string & name) const throw();
+    xbool GetDefineFlagX (std::string_view name) const  noexcept;
     /// Returns string list flag, empty array if not exist
     const Array<std::string> & GetStringListFlag (const std::string & name) const;
     /// Returns num list flag, empty array if not exist
@@ -144,16 +144,16 @@ namespace ngcore
 
 
     /// Test, if string flag is defined
-    bool StringFlagDefined (const std::string & name) const;
+    bool StringFlagDefined (std::string_view name) const noexcept;
     /// Test, if num flag is defined
-    bool NumFlagDefined (const std::string & name) const;
+    bool NumFlagDefined (std::string_view name) const noexcept;
     /// Test, if num flag is defined
-    bool FlagsFlagDefined (const std::string & name) const;
+    bool FlagsFlagDefined (std::string_view name) const noexcept;
     /// Test, if string list flag is defined
-    bool StringListFlagDefined (const std::string & name) const;
+    bool StringListFlagDefined (std::string_view name) const noexcept;
     /// Test, if num list flag is defined
-    bool NumListFlagDefined (const std::string & name) const;
+    bool NumListFlagDefined (std::string_view name) const noexcept;
-    bool AnyFlagDefined (const std::string& name) const;
+    bool AnyFlagDefined (std::string_view name) const noexcept;
 
     /// number of string flags
     int GetNStringFlags () const { return strflags.Size(); }

libsrc/core/generate_mpi_sources.py

@@ -5,6 +5,8 @@ functions = [
         ("int", "MPI_Alltoall", "void*", "int", "MPI_Datatype", "void*", "int", "MPI_Datatype", "MPI_Comm"),
         ("int", "MPI_Barrier", "MPI_Comm"),
         ("int", "MPI_Bcast", "void*", "int", "MPI_Datatype", "int", "MPI_Comm"),
+        ("int", "MPI_Comm_c2f", "MPI_Comm"),
+        ("int", "MPI_Comm_create", "MPI_Comm", "MPI_Group", "MPI_Comm*"),
         ("int", "MPI_Comm_create_group", "MPI_Comm", "MPI_Group", "int", "MPI_Comm*"),
         ("int", "MPI_Comm_free", "MPI_Comm*"),
         ("int", "MPI_Comm_group", "MPI_Comm", "MPI_Group*"),
@@ -12,6 +14,7 @@ functions = [
         ("int", "MPI_Comm_size", "MPI_Comm", "int*"),
         ("int", "MPI_Finalize"),
         ("int", "MPI_Gather", "void*", "int", "MPI_Datatype", "void*", "int", "MPI_Datatype", "int", "MPI_Comm"),
+        ("int", "MPI_Gatherv", "void*", "int", "MPI_Datatype", "void*", "int*", "int*", "MPI_Datatype", "int", "MPI_Comm"),
         ("int", "MPI_Get_count", "MPI_Status*", "MPI_Datatype", "int*"),
         ("int", "MPI_Get_processor_name", "char*", "int*"),
         ("int", "MPI_Group_incl", "MPI_Group", "int", "int*", "MPI_Group*"),
@@ -24,11 +27,14 @@ functions = [
         ("int", "MPI_Probe", "int", "int", "MPI_Comm", "MPI_Status*"),
         ("int", "MPI_Query_thread", "int*"),
         ("int", "MPI_Recv", "void*", "int", "MPI_Datatype", "int", "int", "MPI_Comm", "MPI_Status*"),
+        ("int", "MPI_Recv_init", "void*", "int", "MPI_Datatype", "int", "int", "MPI_Comm", "MPI_Request*"),
         ("int", "MPI_Reduce", "void*", "void*", "int", "MPI_Datatype", "MPI_Op", "int", "MPI_Comm"),
         ("int", "MPI_Reduce_local", "void*", "void*", "int", "MPI_Datatype", "MPI_Op"),
         ("int", "MPI_Request_free", "MPI_Request*"),
         ("int", "MPI_Scatter", "void*", "int", "MPI_Datatype", "void*", "int", "MPI_Datatype", "int", "MPI_Comm"),
         ("int", "MPI_Send", "void*", "int", "MPI_Datatype", "int", "int", "MPI_Comm"),
+        ("int", "MPI_Send_init", "void*", "int", "MPI_Datatype", "int", "int", "MPI_Comm", "MPI_Request*"),
+        ("int", "MPI_Startall", "int", "MPI_Request*:0"),
         ("int", "MPI_Type_commit", "MPI_Datatype*"),
         ("int", "MPI_Type_contiguous", "int", "MPI_Datatype", "MPI_Datatype*"),
         ("int", "MPI_Type_create_resized", "MPI_Datatype", "MPI_Aint", "MPI_Aint", "MPI_Datatype*"),
@@ -40,7 +46,6 @@ functions = [
         ("int", "MPI_Wait", "MPI_Request*", "MPI_Status*"),
         ("int", "MPI_Waitall", "int", "MPI_Request*:0", "MPI_Status*"),
         ("int", "MPI_Waitany", "int", "MPI_Request*:0", "int*", "MPI_Status*"),
-        ("int", "MPI_Comm_c2f", "MPI_Comm"),
         ]
 
 constants = [
@@ -51,6 +56,7 @@ constants = [
         ("MPI_Datatype", "MPI_C_BOOL"),
         ("MPI_Datatype", "MPI_DATATYPE_NULL"),
         ("MPI_Datatype", "MPI_DOUBLE"),
+        ("MPI_Datatype", "MPI_FLOAT"),
         ("MPI_Datatype", "MPI_INT"),
         ("MPI_Datatype", "MPI_SHORT"),
         ("MPI_Datatype", "MPI_UINT64_T"),
@@ -58,6 +64,7 @@ constants = [
         ("MPI_Op", "MPI_MAX"),
         ("MPI_Op", "MPI_MIN"),
         ("MPI_Op", "MPI_SUM"),
+        ("MPI_Request", "MPI_REQUEST_NULL"),
         ("MPI_Status*", "MPI_STATUSES_IGNORE"),
         ("MPI_Status*", "MPI_STATUS_IGNORE"),
         ("int", "MPI_ANY_SOURCE"),

libsrc/core/hashtable.hpp

@@ -46,12 +46,23 @@ namespace ngcore
   class IVec
   {
     /// data
-    T i[(N>0)?N:1];
+    // T i[(N>0)?N:1];
+
+    HTArray<N,T> i;
     
   public:
     ///
-    NETGEN_INLINE IVec () { }
+    constexpr NETGEN_INLINE IVec () = default;
+    constexpr NETGEN_INLINE IVec (const IVec & i1) : i(i1.i) { }
 
+    constexpr NETGEN_INLINE IVec (T ai1) : i(ai1) { }
+    
+    template <class... T2,
+              std::enable_if_t<N==1+sizeof...(T2),bool> = true>
+    constexpr IVec (const T &v, T2... rest)
+      : i{v,rest...} { } 
+
+    /*
     /// init all
     NETGEN_INLINE IVec (T ai1)
     { 
@@ -77,11 +88,13 @@ namespace ngcore
     /// init i[0], i[1], i[2]
     NETGEN_INLINE IVec (T ai1, T ai2, T ai3, T ai4, T ai5, T ai6, T ai7, T ai8, T ai9)
       : i{ai1, ai2, ai3, ai4, ai5, ai6, ai7, ai8, ai9 } { ; }            
+    */
     
     template <typename ARCHIVE>
     void DoArchive(ARCHIVE& ar)
     {
-      ar.Do(i, N);
+      // ar.Do(i.begin(), N);
+      ar.Do(i.Ptr(), N);
     }
 
     template <int N2, typename T2>

libsrc/core/memtracer.hpp

@@ -35,7 +35,7 @@ namespace ngcore
 
   class MemoryTracer
   {
-    #ifdef NETGEN_TRACE_MEMORY
+    #if defined(NETGEN_TRACE_MEMORY) && !defined(__CUDA_ARCH__)
     NGCORE_API static std::vector<std::string> names;
     NGCORE_API static std::vector<int> parents;
 
@@ -148,7 +148,7 @@ namespace ngcore
 
     static const std::vector<std::string> & GetNames() { return names; }
     static const std::vector<int> & GetParents() { return parents; }
-#else // NETGEN_TRACE_MEMORY
+#else // defined(NETGEN_TRACE_MEMORY) && !defined(__CUDA_ARCH__)
   public:
     MemoryTracer() {}
     MemoryTracer( std::string /* name */ ) {}

libsrc/core/ng_mpi.cpp

@@ -7,16 +7,16 @@
 
 #include "array.hpp"
 #include "ngcore_api.hpp"
-#include "pybind11/pytypes.h"
 
 #ifdef NG_PYTHON
+#include "pybind11/pytypes.h"
 #include "python_ngcore.hpp"
-#endif
 
 #define MPI4PY_LIMITED_API 1
 #define MPI4PY_LIMITED_API_SKIP_MESSAGE 1
 #define MPI4PY_LIMITED_API_SKIP_SESSION 1
 #include "mpi4py_pycapi.h"  // mpi4py < 4.0.0
+#endif
 
 #ifdef MSMPI_VER
 int MPI_Comm_create_group(MPI_Comm arg0, MPI_Group arg1, int arg2,
@@ -49,7 +49,7 @@ void gather_strided_array(size_t count, char* data) {
   if constexpr (size < stride) {
     char* dst = data;
     char* src = data;
-    for (auto i : Range(count)) {
+    for ( [[maybe_unused]] auto i : Range(count)) {
       memcpy(dst, src, size);
       dst += size;
       src += stride;
@@ -164,7 +164,7 @@ void ng_init_mpi() {
       imported_mpi4py = true;
     }
     PyObject* py_src = src.ptr();
-    auto type = Py_TYPE(py_src);
+    [[maybe_unused]] auto type = Py_TYPE(py_src);
     if (PyObject_TypeCheck(py_src, &PyMPIComm_Type)) {
       dst = mpi2ng(*PyMPIComm_Get(py_src));
       return !PyErr_Occurred();

libsrc/core/ng_mpi.hpp

@@ -54,6 +54,7 @@ struct NG_MPI_Request {
   NG_MPI_Request() = default;
   NG_MPI_Request(uintptr_t value_) : value(value_) {}
   void operator=(uintptr_t value_) { value = value_; }
+  void operator=(void *value_) { value = reinterpret_cast<uintptr_t>(value_); }
 };
 
 struct NG_MPI_Op {

libsrc/core/ng_mpi_generated_declarations.hpp

@@ -5,6 +5,8 @@ NGCORE_API extern int (*NG_MPI_Allreduce)(void*, void*, int, NG_MPI_Datatype, NG
 NGCORE_API extern int (*NG_MPI_Alltoall)(void*, int, NG_MPI_Datatype, void*, int, NG_MPI_Datatype, NG_MPI_Comm);
 NGCORE_API extern int (*NG_MPI_Barrier)(NG_MPI_Comm);
 NGCORE_API extern int (*NG_MPI_Bcast)(void*, int, NG_MPI_Datatype, int, NG_MPI_Comm);
+NGCORE_API extern int (*NG_MPI_Comm_c2f)(NG_MPI_Comm);
+NGCORE_API extern int (*NG_MPI_Comm_create)(NG_MPI_Comm, NG_MPI_Group, NG_MPI_Comm*);
 NGCORE_API extern int (*NG_MPI_Comm_create_group)(NG_MPI_Comm, NG_MPI_Group, int, NG_MPI_Comm*);
 NGCORE_API extern int (*NG_MPI_Comm_free)(NG_MPI_Comm*);
 NGCORE_API extern int (*NG_MPI_Comm_group)(NG_MPI_Comm, NG_MPI_Group*);
@@ -12,6 +14,7 @@ NGCORE_API extern int (*NG_MPI_Comm_rank)(NG_MPI_Comm, int*);
 NGCORE_API extern int (*NG_MPI_Comm_size)(NG_MPI_Comm, int*);
 NGCORE_API extern int (*NG_MPI_Finalize)();
 NGCORE_API extern int (*NG_MPI_Gather)(void*, int, NG_MPI_Datatype, void*, int, NG_MPI_Datatype, int, NG_MPI_Comm);
+NGCORE_API extern int (*NG_MPI_Gatherv)(void*, int, NG_MPI_Datatype, void*, int*, int*, NG_MPI_Datatype, int, NG_MPI_Comm);
 NGCORE_API extern int (*NG_MPI_Get_count)(NG_MPI_Status*, NG_MPI_Datatype, int*);
 NGCORE_API extern int (*NG_MPI_Get_processor_name)(char*, int*);
 NGCORE_API extern int (*NG_MPI_Group_incl)(NG_MPI_Group, int, int*, NG_MPI_Group*);
@@ -24,11 +27,14 @@ NGCORE_API extern int (*NG_MPI_Isend)(void*, int, NG_MPI_Datatype, int, int, NG_
 NGCORE_API extern int (*NG_MPI_Probe)(int, int, NG_MPI_Comm, NG_MPI_Status*);
 NGCORE_API extern int (*NG_MPI_Query_thread)(int*);
 NGCORE_API extern int (*NG_MPI_Recv)(void*, int, NG_MPI_Datatype, int, int, NG_MPI_Comm, NG_MPI_Status*);
+NGCORE_API extern int (*NG_MPI_Recv_init)(void*, int, NG_MPI_Datatype, int, int, NG_MPI_Comm, NG_MPI_Request*);
 NGCORE_API extern int (*NG_MPI_Reduce)(void*, void*, int, NG_MPI_Datatype, NG_MPI_Op, int, NG_MPI_Comm);
 NGCORE_API extern int (*NG_MPI_Reduce_local)(void*, void*, int, NG_MPI_Datatype, NG_MPI_Op);
 NGCORE_API extern int (*NG_MPI_Request_free)(NG_MPI_Request*);
 NGCORE_API extern int (*NG_MPI_Scatter)(void*, int, NG_MPI_Datatype, void*, int, NG_MPI_Datatype, int, NG_MPI_Comm);
 NGCORE_API extern int (*NG_MPI_Send)(void*, int, NG_MPI_Datatype, int, int, NG_MPI_Comm);
+NGCORE_API extern int (*NG_MPI_Send_init)(void*, int, NG_MPI_Datatype, int, int, NG_MPI_Comm, NG_MPI_Request*);
+NGCORE_API extern int (*NG_MPI_Startall)(int, NG_MPI_Request*);
 NGCORE_API extern int (*NG_MPI_Type_commit)(NG_MPI_Datatype*);
 NGCORE_API extern int (*NG_MPI_Type_contiguous)(int, NG_MPI_Datatype, NG_MPI_Datatype*);
 NGCORE_API extern int (*NG_MPI_Type_create_resized)(NG_MPI_Datatype, NG_MPI_Aint, NG_MPI_Aint, NG_MPI_Datatype*);
@@ -40,7 +46,6 @@ NGCORE_API extern int (*NG_MPI_Type_size)(NG_MPI_Datatype, int*);
 NGCORE_API extern int (*NG_MPI_Wait)(NG_MPI_Request*, NG_MPI_Status*);
 NGCORE_API extern int (*NG_MPI_Waitall)(int, NG_MPI_Request*, NG_MPI_Status*);
 NGCORE_API extern int (*NG_MPI_Waitany)(int, NG_MPI_Request*, int*, NG_MPI_Status*);
-NGCORE_API extern int (*NG_MPI_Comm_c2f)(NG_MPI_Comm);
 NGCORE_API extern NG_MPI_Comm NG_MPI_COMM_NULL;
 NGCORE_API extern NG_MPI_Comm NG_MPI_COMM_WORLD;
 NGCORE_API extern NG_MPI_Datatype NG_MPI_CHAR;
@@ -48,6 +53,7 @@ NGCORE_API extern NG_MPI_Datatype NG_MPI_CXX_DOUBLE_COMPLEX;
 NGCORE_API extern NG_MPI_Datatype NG_MPI_C_BOOL;
 NGCORE_API extern NG_MPI_Datatype NG_MPI_DATATYPE_NULL;
 NGCORE_API extern NG_MPI_Datatype NG_MPI_DOUBLE;
+NGCORE_API extern NG_MPI_Datatype NG_MPI_FLOAT;
 NGCORE_API extern NG_MPI_Datatype NG_MPI_INT;
 NGCORE_API extern NG_MPI_Datatype NG_MPI_SHORT;
 NGCORE_API extern NG_MPI_Datatype NG_MPI_UINT64_T;
@@ -55,6 +61,7 @@ NGCORE_API extern NG_MPI_Op NG_MPI_LOR;
 NGCORE_API extern NG_MPI_Op NG_MPI_MAX;
 NGCORE_API extern NG_MPI_Op NG_MPI_MIN;
 NGCORE_API extern NG_MPI_Op NG_MPI_SUM;
+NGCORE_API extern NG_MPI_Request NG_MPI_REQUEST_NULL;
 NGCORE_API extern NG_MPI_Status* NG_MPI_STATUSES_IGNORE;
 NGCORE_API extern NG_MPI_Status* NG_MPI_STATUS_IGNORE;
 NGCORE_API extern int NG_MPI_ANY_SOURCE;
@@ -74,6 +81,8 @@ static const auto NG_MPI_Allreduce = MPI_Allreduce;
 static const auto NG_MPI_Alltoall = MPI_Alltoall;
 static const auto NG_MPI_Barrier = MPI_Barrier;
 static const auto NG_MPI_Bcast = MPI_Bcast;
+static const auto NG_MPI_Comm_c2f = MPI_Comm_c2f;
+static const auto NG_MPI_Comm_create = MPI_Comm_create;
 static const auto NG_MPI_Comm_create_group = MPI_Comm_create_group;
 static const auto NG_MPI_Comm_free = MPI_Comm_free;
 static const auto NG_MPI_Comm_group = MPI_Comm_group;
@@ -81,6 +90,7 @@ static const auto NG_MPI_Comm_rank = MPI_Comm_rank;
 static const auto NG_MPI_Comm_size = MPI_Comm_size;
 static const auto NG_MPI_Finalize = MPI_Finalize;
 static const auto NG_MPI_Gather = MPI_Gather;
+static const auto NG_MPI_Gatherv = MPI_Gatherv;
 static const auto NG_MPI_Get_count = MPI_Get_count;
 static const auto NG_MPI_Get_processor_name = MPI_Get_processor_name;
 static const auto NG_MPI_Group_incl = MPI_Group_incl;
@@ -93,11 +103,14 @@ static const auto NG_MPI_Isend = MPI_Isend;
 static const auto NG_MPI_Probe = MPI_Probe;
 static const auto NG_MPI_Query_thread = MPI_Query_thread;
 static const auto NG_MPI_Recv = MPI_Recv;
+static const auto NG_MPI_Recv_init = MPI_Recv_init;
 static const auto NG_MPI_Reduce = MPI_Reduce;
 static const auto NG_MPI_Reduce_local = MPI_Reduce_local;
 static const auto NG_MPI_Request_free = MPI_Request_free;
 static const auto NG_MPI_Scatter = MPI_Scatter;
 static const auto NG_MPI_Send = MPI_Send;
+static const auto NG_MPI_Send_init = MPI_Send_init;
+static const auto NG_MPI_Startall = MPI_Startall;
 static const auto NG_MPI_Type_commit = MPI_Type_commit;
 static const auto NG_MPI_Type_contiguous = MPI_Type_contiguous;
 static const auto NG_MPI_Type_create_resized = MPI_Type_create_resized;
@@ -109,7 +122,6 @@ static const auto NG_MPI_Type_size = MPI_Type_size;
 static const auto NG_MPI_Wait = MPI_Wait;
 static const auto NG_MPI_Waitall = MPI_Waitall;
 static const auto NG_MPI_Waitany = MPI_Waitany;
-static const auto NG_MPI_Comm_c2f = MPI_Comm_c2f;
 static const decltype(MPI_COMM_NULL) NG_MPI_COMM_NULL = MPI_COMM_NULL;
 static const decltype(MPI_COMM_WORLD) NG_MPI_COMM_WORLD = MPI_COMM_WORLD;
 static const decltype(MPI_CHAR) NG_MPI_CHAR = MPI_CHAR;
@@ -117,6 +129,7 @@ static const decltype(MPI_CXX_DOUBLE_COMPLEX) NG_MPI_CXX_DOUBLE_COMPLEX = MPI_CX
 static const decltype(MPI_C_BOOL) NG_MPI_C_BOOL = MPI_C_BOOL;
 static const decltype(MPI_DATATYPE_NULL) NG_MPI_DATATYPE_NULL = MPI_DATATYPE_NULL;
 static const decltype(MPI_DOUBLE) NG_MPI_DOUBLE = MPI_DOUBLE;
+static const decltype(MPI_FLOAT) NG_MPI_FLOAT = MPI_FLOAT;
 static const decltype(MPI_INT) NG_MPI_INT = MPI_INT;
 static const decltype(MPI_SHORT) NG_MPI_SHORT = MPI_SHORT;
 static const decltype(MPI_UINT64_T) NG_MPI_UINT64_T = MPI_UINT64_T;
@@ -124,6 +137,7 @@ static const decltype(MPI_LOR) NG_MPI_LOR = MPI_LOR;
 static const decltype(MPI_MAX) NG_MPI_MAX = MPI_MAX;
 static const decltype(MPI_MIN) NG_MPI_MIN = MPI_MIN;
 static const decltype(MPI_SUM) NG_MPI_SUM = MPI_SUM;
+static const decltype(MPI_REQUEST_NULL) NG_MPI_REQUEST_NULL = MPI_REQUEST_NULL;
 static const decltype(MPI_STATUSES_IGNORE) NG_MPI_STATUSES_IGNORE = MPI_STATUSES_IGNORE;
 static const decltype(MPI_STATUS_IGNORE) NG_MPI_STATUS_IGNORE = MPI_STATUS_IGNORE;
 static const decltype(MPI_ANY_SOURCE) NG_MPI_ANY_SOURCE = MPI_ANY_SOURCE;

libsrc/core/ng_mpi_generated_dummy_init.hpp

@@ -4,6 +4,8 @@ decltype(NG_MPI_Allreduce) NG_MPI_Allreduce = [](void*, void*, int, NG_MPI_Datat
 decltype(NG_MPI_Alltoall) NG_MPI_Alltoall = [](void*, int, NG_MPI_Datatype, void*, int, NG_MPI_Datatype, NG_MPI_Comm)->int { throw no_mpi(); };
 decltype(NG_MPI_Barrier) NG_MPI_Barrier = [](NG_MPI_Comm)->int { throw no_mpi(); };
 decltype(NG_MPI_Bcast) NG_MPI_Bcast = [](void*, int, NG_MPI_Datatype, int, NG_MPI_Comm)->int { throw no_mpi(); };
+decltype(NG_MPI_Comm_c2f) NG_MPI_Comm_c2f = [](NG_MPI_Comm)->int { throw no_mpi(); };
+decltype(NG_MPI_Comm_create) NG_MPI_Comm_create = [](NG_MPI_Comm, NG_MPI_Group, NG_MPI_Comm*)->int { throw no_mpi(); };
 decltype(NG_MPI_Comm_create_group) NG_MPI_Comm_create_group = [](NG_MPI_Comm, NG_MPI_Group, int, NG_MPI_Comm*)->int { throw no_mpi(); };
 decltype(NG_MPI_Comm_free) NG_MPI_Comm_free = [](NG_MPI_Comm*)->int { throw no_mpi(); };
 decltype(NG_MPI_Comm_group) NG_MPI_Comm_group = [](NG_MPI_Comm, NG_MPI_Group*)->int { throw no_mpi(); };
@@ -11,6 +13,7 @@ decltype(NG_MPI_Comm_rank) NG_MPI_Comm_rank = [](NG_MPI_Comm, int*)->int { throw
 decltype(NG_MPI_Comm_size) NG_MPI_Comm_size = [](NG_MPI_Comm, int*)->int { throw no_mpi(); };
 decltype(NG_MPI_Finalize) NG_MPI_Finalize = []()->int { throw no_mpi(); };
 decltype(NG_MPI_Gather) NG_MPI_Gather = [](void*, int, NG_MPI_Datatype, void*, int, NG_MPI_Datatype, int, NG_MPI_Comm)->int { throw no_mpi(); };
+decltype(NG_MPI_Gatherv) NG_MPI_Gatherv = [](void*, int, NG_MPI_Datatype, void*, int*, int*, NG_MPI_Datatype, int, NG_MPI_Comm)->int { throw no_mpi(); };
 decltype(NG_MPI_Get_count) NG_MPI_Get_count = [](NG_MPI_Status*, NG_MPI_Datatype, int*)->int { throw no_mpi(); };
 decltype(NG_MPI_Get_processor_name) NG_MPI_Get_processor_name = [](char*, int*)->int { throw no_mpi(); };
 decltype(NG_MPI_Group_incl) NG_MPI_Group_incl = [](NG_MPI_Group, int, int*, NG_MPI_Group*)->int { throw no_mpi(); };
@@ -23,11 +26,14 @@ decltype(NG_MPI_Isend) NG_MPI_Isend = [](void*, int, NG_MPI_Datatype, int, int,
 decltype(NG_MPI_Probe) NG_MPI_Probe = [](int, int, NG_MPI_Comm, NG_MPI_Status*)->int { throw no_mpi(); };
 decltype(NG_MPI_Query_thread) NG_MPI_Query_thread = [](int*)->int { throw no_mpi(); };
 decltype(NG_MPI_Recv) NG_MPI_Recv = [](void*, int, NG_MPI_Datatype, int, int, NG_MPI_Comm, NG_MPI_Status*)->int { throw no_mpi(); };
+decltype(NG_MPI_Recv_init) NG_MPI_Recv_init = [](void*, int, NG_MPI_Datatype, int, int, NG_MPI_Comm, NG_MPI_Request*)->int { throw no_mpi(); };
 decltype(NG_MPI_Reduce) NG_MPI_Reduce = [](void*, void*, int, NG_MPI_Datatype, NG_MPI_Op, int, NG_MPI_Comm)->int { throw no_mpi(); };
 decltype(NG_MPI_Reduce_local) NG_MPI_Reduce_local = [](void*, void*, int, NG_MPI_Datatype, NG_MPI_Op)->int { throw no_mpi(); };
 decltype(NG_MPI_Request_free) NG_MPI_Request_free = [](NG_MPI_Request*)->int { throw no_mpi(); };
 decltype(NG_MPI_Scatter) NG_MPI_Scatter = [](void*, int, NG_MPI_Datatype, void*, int, NG_MPI_Datatype, int, NG_MPI_Comm)->int { throw no_mpi(); };
 decltype(NG_MPI_Send) NG_MPI_Send = [](void*, int, NG_MPI_Datatype, int, int, NG_MPI_Comm)->int { throw no_mpi(); };
+decltype(NG_MPI_Send_init) NG_MPI_Send_init = [](void*, int, NG_MPI_Datatype, int, int, NG_MPI_Comm, NG_MPI_Request*)->int { throw no_mpi(); };
+decltype(NG_MPI_Startall) NG_MPI_Startall = [](int, NG_MPI_Request*)->int { throw no_mpi(); };
 decltype(NG_MPI_Type_commit) NG_MPI_Type_commit = [](NG_MPI_Datatype*)->int { throw no_mpi(); };
 decltype(NG_MPI_Type_contiguous) NG_MPI_Type_contiguous = [](int, NG_MPI_Datatype, NG_MPI_Datatype*)->int { throw no_mpi(); };
 decltype(NG_MPI_Type_create_resized) NG_MPI_Type_create_resized = [](NG_MPI_Datatype, NG_MPI_Aint, NG_MPI_Aint, NG_MPI_Datatype*)->int { throw no_mpi(); };
@@ -39,7 +45,6 @@ decltype(NG_MPI_Type_size) NG_MPI_Type_size = [](NG_MPI_Datatype, int*)->int { t
 decltype(NG_MPI_Wait) NG_MPI_Wait = [](NG_MPI_Request*, NG_MPI_Status*)->int { throw no_mpi(); };
 decltype(NG_MPI_Waitall) NG_MPI_Waitall = [](int, NG_MPI_Request*, NG_MPI_Status*)->int { throw no_mpi(); };
 decltype(NG_MPI_Waitany) NG_MPI_Waitany = [](int, NG_MPI_Request*, int*, NG_MPI_Status*)->int { throw no_mpi(); };
-decltype(NG_MPI_Comm_c2f) NG_MPI_Comm_c2f = [](NG_MPI_Comm)->int { throw no_mpi(); };
 NG_MPI_Comm NG_MPI_COMM_NULL = 0;
 NG_MPI_Comm NG_MPI_COMM_WORLD = 0;
 NG_MPI_Datatype NG_MPI_CHAR = 0;
@@ -47,6 +52,7 @@ NG_MPI_Datatype NG_MPI_CXX_DOUBLE_COMPLEX = 0;
 NG_MPI_Datatype NG_MPI_C_BOOL = 0;
 NG_MPI_Datatype NG_MPI_DATATYPE_NULL = 0;
 NG_MPI_Datatype NG_MPI_DOUBLE = 0;
+NG_MPI_Datatype NG_MPI_FLOAT = 0;
 NG_MPI_Datatype NG_MPI_INT = 0;
 NG_MPI_Datatype NG_MPI_SHORT = 0;
 NG_MPI_Datatype NG_MPI_UINT64_T = 0;
@@ -54,6 +60,7 @@ NG_MPI_Op NG_MPI_LOR = 0;
 NG_MPI_Op NG_MPI_MAX = 0;
 NG_MPI_Op NG_MPI_MIN = 0;
 NG_MPI_Op NG_MPI_SUM = 0;
+NG_MPI_Request NG_MPI_REQUEST_NULL = 0;
 NG_MPI_Status* NG_MPI_STATUSES_IGNORE = 0;
 NG_MPI_Status* NG_MPI_STATUS_IGNORE = 0;
 int NG_MPI_ANY_SOURCE = 0;

libsrc/core/ng_mpi_generated_init.hpp

@@ -4,6 +4,8 @@ NG_MPI_Allreduce = [](void* arg0, void* arg1, int arg2, NG_MPI_Datatype arg3, NG
 NG_MPI_Alltoall = [](void* arg0, int arg1, NG_MPI_Datatype arg2, void* arg3, int arg4, NG_MPI_Datatype arg5, NG_MPI_Comm arg6)->int { return MPI_Alltoall( arg0,  arg1,  ng2mpi(arg2),  arg3,  arg4,  ng2mpi(arg5),  ng2mpi(arg6)); };
 NG_MPI_Barrier = [](NG_MPI_Comm arg0)->int { return MPI_Barrier( ng2mpi(arg0)); };
 NG_MPI_Bcast = [](void* arg0, int arg1, NG_MPI_Datatype arg2, int arg3, NG_MPI_Comm arg4)->int { return MPI_Bcast( arg0,  arg1,  ng2mpi(arg2),  arg3,  ng2mpi(arg4)); };
+NG_MPI_Comm_c2f = [](NG_MPI_Comm arg0)->int { return MPI_Comm_c2f( ng2mpi(arg0)); };
+NG_MPI_Comm_create = [](NG_MPI_Comm arg0, NG_MPI_Group arg1, NG_MPI_Comm* arg2)->int { return MPI_Comm_create( ng2mpi(arg0),  ng2mpi(arg1),  ng2mpi(arg2)); };
 NG_MPI_Comm_create_group = [](NG_MPI_Comm arg0, NG_MPI_Group arg1, int arg2, NG_MPI_Comm* arg3)->int { return MPI_Comm_create_group( ng2mpi(arg0),  ng2mpi(arg1),  arg2,  ng2mpi(arg3)); };
 NG_MPI_Comm_free = [](NG_MPI_Comm* arg0)->int { return MPI_Comm_free( ng2mpi(arg0)); };
 NG_MPI_Comm_group = [](NG_MPI_Comm arg0, NG_MPI_Group* arg1)->int { return MPI_Comm_group( ng2mpi(arg0),  ng2mpi(arg1)); };
@@ -11,6 +13,7 @@ NG_MPI_Comm_rank = [](NG_MPI_Comm arg0, int* arg1)->int { return MPI_Comm_rank(
 NG_MPI_Comm_size = [](NG_MPI_Comm arg0, int* arg1)->int { return MPI_Comm_size( ng2mpi(arg0),  arg1); };
 NG_MPI_Finalize = []()->int { return MPI_Finalize(); };
 NG_MPI_Gather = [](void* arg0, int arg1, NG_MPI_Datatype arg2, void* arg3, int arg4, NG_MPI_Datatype arg5, int arg6, NG_MPI_Comm arg7)->int { return MPI_Gather( arg0,  arg1,  ng2mpi(arg2),  arg3,  arg4,  ng2mpi(arg5),  arg6,  ng2mpi(arg7)); };
+NG_MPI_Gatherv = [](void* arg0, int arg1, NG_MPI_Datatype arg2, void* arg3, int* arg4, int* arg5, NG_MPI_Datatype arg6, int arg7, NG_MPI_Comm arg8)->int { return MPI_Gatherv( arg0,  arg1,  ng2mpi(arg2),  arg3,  arg4,  arg5,  ng2mpi(arg6),  arg7,  ng2mpi(arg8)); };
 NG_MPI_Get_count = [](NG_MPI_Status* arg0, NG_MPI_Datatype arg1, int* arg2)->int { return MPI_Get_count( ng2mpi(arg0),  ng2mpi(arg1),  arg2); };
 NG_MPI_Get_processor_name = [](char* arg0, int* arg1)->int { return MPI_Get_processor_name( arg0,  arg1); };
 NG_MPI_Group_incl = [](NG_MPI_Group arg0, int arg1, int* arg2, NG_MPI_Group* arg3)->int { return MPI_Group_incl( ng2mpi(arg0),  arg1,  arg2,  ng2mpi(arg3)); };
@@ -23,11 +26,14 @@ NG_MPI_Isend = [](void* arg0, int arg1, NG_MPI_Datatype arg2, int arg3, int arg4
 NG_MPI_Probe = [](int arg0, int arg1, NG_MPI_Comm arg2, NG_MPI_Status* arg3)->int { return MPI_Probe( arg0,  arg1,  ng2mpi(arg2),  ng2mpi(arg3)); };
 NG_MPI_Query_thread = [](int* arg0)->int { return MPI_Query_thread( arg0); };
 NG_MPI_Recv = [](void* arg0, int arg1, NG_MPI_Datatype arg2, int arg3, int arg4, NG_MPI_Comm arg5, NG_MPI_Status* arg6)->int { return MPI_Recv( arg0,  arg1,  ng2mpi(arg2),  arg3,  arg4,  ng2mpi(arg5),  ng2mpi(arg6)); };
+NG_MPI_Recv_init = [](void* arg0, int arg1, NG_MPI_Datatype arg2, int arg3, int arg4, NG_MPI_Comm arg5, NG_MPI_Request* arg6)->int { return MPI_Recv_init( arg0,  arg1,  ng2mpi(arg2),  arg3,  arg4,  ng2mpi(arg5),  ng2mpi(arg6)); };
 NG_MPI_Reduce = [](void* arg0, void* arg1, int arg2, NG_MPI_Datatype arg3, NG_MPI_Op arg4, int arg5, NG_MPI_Comm arg6)->int { return MPI_Reduce( arg0,  arg1,  arg2,  ng2mpi(arg3),  ng2mpi(arg4),  arg5,  ng2mpi(arg6)); };
 NG_MPI_Reduce_local = [](void* arg0, void* arg1, int arg2, NG_MPI_Datatype arg3, NG_MPI_Op arg4)->int { return MPI_Reduce_local( arg0,  arg1,  arg2,  ng2mpi(arg3),  ng2mpi(arg4)); };
 NG_MPI_Request_free = [](NG_MPI_Request* arg0)->int { return MPI_Request_free( ng2mpi(arg0)); };
 NG_MPI_Scatter = [](void* arg0, int arg1, NG_MPI_Datatype arg2, void* arg3, int arg4, NG_MPI_Datatype arg5, int arg6, NG_MPI_Comm arg7)->int { return MPI_Scatter( arg0,  arg1,  ng2mpi(arg2),  arg3,  arg4,  ng2mpi(arg5),  arg6,  ng2mpi(arg7)); };
 NG_MPI_Send = [](void* arg0, int arg1, NG_MPI_Datatype arg2, int arg3, int arg4, NG_MPI_Comm arg5)->int { return MPI_Send( arg0,  arg1,  ng2mpi(arg2),  arg3,  arg4,  ng2mpi(arg5)); };
+NG_MPI_Send_init = [](void* arg0, int arg1, NG_MPI_Datatype arg2, int arg3, int arg4, NG_MPI_Comm arg5, NG_MPI_Request* arg6)->int { return MPI_Send_init( arg0,  arg1,  ng2mpi(arg2),  arg3,  arg4,  ng2mpi(arg5),  ng2mpi(arg6)); };
+NG_MPI_Startall = [](int arg0, NG_MPI_Request* arg1)->int { return MPI_Startall( arg0,  ng2mpi(arg1, arg0)); };
 NG_MPI_Type_commit = [](NG_MPI_Datatype* arg0)->int { return MPI_Type_commit( ng2mpi(arg0)); };
 NG_MPI_Type_contiguous = [](int arg0, NG_MPI_Datatype arg1, NG_MPI_Datatype* arg2)->int { return MPI_Type_contiguous( arg0,  ng2mpi(arg1),  ng2mpi(arg2)); };
 NG_MPI_Type_create_resized = [](NG_MPI_Datatype arg0, NG_MPI_Aint arg1, NG_MPI_Aint arg2, NG_MPI_Datatype* arg3)->int { return MPI_Type_create_resized( ng2mpi(arg0),  ng2mpi(arg1),  ng2mpi(arg2),  ng2mpi(arg3)); };
@@ -39,7 +45,6 @@ NG_MPI_Type_size = [](NG_MPI_Datatype arg0, int* arg1)->int { return MPI_Type_si
 NG_MPI_Wait = [](NG_MPI_Request* arg0, NG_MPI_Status* arg1)->int { return MPI_Wait( ng2mpi(arg0),  ng2mpi(arg1)); };
 NG_MPI_Waitall = [](int arg0, NG_MPI_Request* arg1, NG_MPI_Status* arg2)->int { return MPI_Waitall( arg0,  ng2mpi(arg1, arg0),  ng2mpi(arg2)); };
 NG_MPI_Waitany = [](int arg0, NG_MPI_Request* arg1, int* arg2, NG_MPI_Status* arg3)->int { return MPI_Waitany( arg0,  ng2mpi(arg1, arg0),  arg2,  ng2mpi(arg3)); };
-NG_MPI_Comm_c2f = [](NG_MPI_Comm arg0)->int { return MPI_Comm_c2f( ng2mpi(arg0)); };
 NG_MPI_COMM_NULL = mpi2ng(MPI_COMM_NULL);
 NG_MPI_COMM_WORLD = mpi2ng(MPI_COMM_WORLD);
 NG_MPI_CHAR = mpi2ng(MPI_CHAR);
@@ -47,6 +52,7 @@ NG_MPI_CXX_DOUBLE_COMPLEX = mpi2ng(MPI_CXX_DOUBLE_COMPLEX);
 NG_MPI_C_BOOL = mpi2ng(MPI_C_BOOL);
 NG_MPI_DATATYPE_NULL = mpi2ng(MPI_DATATYPE_NULL);
 NG_MPI_DOUBLE = mpi2ng(MPI_DOUBLE);
+NG_MPI_FLOAT = mpi2ng(MPI_FLOAT);
 NG_MPI_INT = mpi2ng(MPI_INT);
 NG_MPI_SHORT = mpi2ng(MPI_SHORT);
 NG_MPI_UINT64_T = mpi2ng(MPI_UINT64_T);
@@ -54,6 +60,7 @@ NG_MPI_LOR = mpi2ng(MPI_LOR);
 NG_MPI_MAX = mpi2ng(MPI_MAX);
 NG_MPI_MIN = mpi2ng(MPI_MIN);
 NG_MPI_SUM = mpi2ng(MPI_SUM);
+NG_MPI_REQUEST_NULL = mpi2ng(MPI_REQUEST_NULL);
 NG_MPI_STATUSES_IGNORE = mpi2ng(MPI_STATUSES_IGNORE);
 NG_MPI_STATUS_IGNORE = mpi2ng(MPI_STATUS_IGNORE);
 NG_MPI_ANY_SOURCE = mpi2ng(MPI_ANY_SOURCE);

libsrc/core/ng_mpi_native.hpp

@@ -0,0 +1,21 @@
+#ifndef NG_MPI_NATIVE_HPP
+#define NG_MPI_NATIVE_HPP
+
+#include <mpi.h>
+
+#include "mpi_wrapper.hpp"
+#include "ng_mpi.hpp"
+
+namespace ngcore {
+
+MPI_Comm NG_MPI_Native(NG_MPI_Comm comm) {
+  return reinterpret_cast<MPI_Comm>(comm.value);
+}
+
+MPI_Comm NG_MPI_Native(NgMPI_Comm comm) {
+  return reinterpret_cast<MPI_Comm>(static_cast<NG_MPI_Comm>(comm).value);
+}
+
+}  // namespace ngcore
+
+#endif  // NG_MPI_NATIVE_HPP

libsrc/core/ng_mpi_wrapper.cpp

@@ -6,7 +6,9 @@
 
 #include "ng_mpi.hpp"
 #include "ngstream.hpp"
+#ifdef NG_PYTHON
 #include "python_ngcore.hpp"
+#endif // NG_PYTHON
 #include "utils.hpp"
 
 using std::cerr;
@@ -14,10 +16,12 @@ using std::cout;
 using std::endl;
 
 #ifndef NG_MPI_WRAPPER
+#ifdef NG_PYTHON
 #define MPI4PY_LIMITED_API 1
 #define MPI4PY_LIMITED_API_SKIP_MESSAGE 1
 #define MPI4PY_LIMITED_API_SKIP_SESSION 1
 #include "mpi4py_pycapi.h"  // mpi4py < 4.0.0
+#endif // NG_PYTHON
 #endif // NG_MPI_WRAPPER
 
 namespace ngcore {
@@ -94,6 +98,7 @@ void InitMPI(std::optional<std::filesystem::path> mpi_lib_path) {
       throw e;
     }
   } else {
+#ifdef NG_PYTHON
     // Use mpi4py to init MPI library and get the vendor name
     auto mpi4py = py::module::import("mpi4py.MPI");
     vendor = mpi4py.attr("get_vendor")()[py::int_(0)].cast<std::string>();
@@ -106,6 +111,7 @@ void InitMPI(std::optional<std::filesystem::path> mpi_lib_path) {
     mpi_lib =
         std::make_unique<SharedLibrary>(mpi4py_lib_file, std::nullopt, true);
 #endif  // WIN32
+#endif // NG_PYTHON
   }
 
   std::string ng_lib_name = "";

libsrc/core/python_ngcore.cpp

@@ -14,13 +14,11 @@ namespace ngcore
   {
     if (py::isinstance<py::dict>(value))
       {             
-        py::dict vdd(value);
+        Flags nested_flags;
-        // call recursively to set dictionary
+        for(auto item : value.cast<py::dict>())
-        for (auto item : vdd) {
+          SetFlag(nested_flags, item.first.cast<string>(),
-          string name = item.first.cast<string>();
+                  item.second.cast<py::object>());
-          py::object val = py::reinterpret_borrow<py::object>(item.second);
+        flags.SetFlag(s, nested_flags);
-          SetFlag(flags, name, val);
-        }
         return;
       }
 
@@ -103,7 +101,9 @@ namespace ngcore
           }
       }
 
-    auto flags = py::cast<Flags>(flags_dict);
+    Flags flags;
+    for(auto item : flags_dict)
+      SetFlag(flags, item.first.cast<string>(), item.second.cast<py::object>());
 
     for (auto item : kwargs)
       {

libsrc/core/python_ngcore_export.cpp

@@ -141,20 +141,28 @@ PYBIND11_MODULE(pyngcore, m) // NOLINT
 
     .def(py::self | py::self)
     .def(py::self & py::self)
-    .def(py::self |= py::self)
+    // .def(py::self |= py::self)   // false clang warnings,
-    .def(py::self &= py::self)
+    // .def(py::self &= py::self)   // see https://github.com/pybind/pybind11/issues/1893
+    .def("__ior__", [](BitArray& lhs, const BitArray& rhs) { return lhs |= rhs; }, py::is_operator())
+    .def("__iand__", [](BitArray& lhs, const BitArray& rhs) { return lhs &= rhs; }, py::is_operator())    
     .def(~py::self)
     ;
 
   py::class_<Flags>(m, "Flags")
     .def(py::init<>())
     .def("__str__", &ToString<Flags>)
-    .def(py::init([](py::object & obj) {
+    .def(py::init([](py::dict kwargs) {
           Flags flags;
-          py::dict d(obj);          
+          for (auto d : kwargs)
-          SetFlag (flags, "", d);
+            SetFlag(flags, d.first.cast<string>(), d.second.cast<py::object>());
           return flags;
-        }), py::arg("obj"), "Create Flags by given object")
+    }), "Create flags from dict")
+    .def(py::init([](py::kwargs kwargs) {
+          Flags flags;
+          for (auto d : kwargs)
+            SetFlag(flags, d.first.cast<string>(), d.second.cast<py::object>());
+          return flags;
+        }), "Create flags from kwargs")
     .def(py::pickle([] (const Flags& self)
         {
           std::stringstream str;

libsrc/core/register_archive.hpp

@@ -34,6 +34,8 @@ namespace ngcore {
     return *this;
   }
 
+  /*
+    // now using has_shared_from_this2 in archive.hpp
   template <typename T>
   struct has_shared_from_this
   {
@@ -42,6 +44,7 @@ namespace ngcore {
     typedef decltype( check<T>(sizeof(char)) ) type;
     static constexpr type value = type();
   };
+  */
 #endif // NETGEN_PYTHON
 
 
@@ -59,7 +62,7 @@ namespace ngcore {
       {
         detail::TCargs<T> args;
         ar &args;
-        auto nT = detail::constructIfPossible<T>(args);
+        auto nT = detail::constructIfPossible<T>(std::move(args));
         return typeid(T) == ti ? nT
           : Archive::Caster<T, Bases>::tryUpcast(ti, nT);
       };

libsrc/core/simd_avx512.hpp

@@ -111,7 +111,7 @@ namespace ngcore
     template <typename Function>
     void SIMD_function (const Function & func, std::true_type)
     {
-      data = (__m512){ func(7), func(6), func(5), func(4),
+      data = (__m512d){ func(7), func(6), func(5), func(4),
                        func(3), func(2), func(1), func(0) };
     }
 

libsrc/core/symboltable.hpp

@@ -45,7 +45,7 @@ namespace ngcore
     }
 
     /// INDEX of symbol name, throws exception if unused
-    size_t Index (const std::string & name) const
+    size_t Index (std::string_view name) const
     {
       for (size_t i = 0; i < names.size(); i++)
         if (names[i] == name) return i;
@@ -53,7 +53,7 @@ namespace ngcore
     }
 
     /// Index of symbol name, returns -1 if unused
-    int CheckIndex (const std::string & name) const
+    int CheckIndex (std::string_view name) const
     {
       for (int i = 0; i < names.size(); i++)
         if (names[i] == name) return i;
@@ -67,12 +67,12 @@ namespace ngcore
     }
 
     /// Returns reference to element. exception for unused identifier
-    reference operator[] (const std::string & name)
+    reference operator[] (std::string_view name)
     {
       return data[Index (name)];
     }
 
-    const_reference operator[] (const std::string & name) const
+    const_reference operator[] (std::string_view name) const
     {
       return data[Index (name)];
     }
@@ -99,7 +99,7 @@ namespace ngcore
     }
 
     /// Associates el to the string name, overrides if name is used
-    void Set (const std::string & name, const T & el)
+    void Set (std::string_view name, const T & el)
     {
       int i = CheckIndex (name);
       if (i >= 0)
@@ -107,14 +107,23 @@ namespace ngcore
       else
         {
           data.push_back(el);
-          names.push_back(name);
+          names.push_back(std::string(name));
         }
     }
 
+
+    
+    /*
     bool Used (const std::string & name) const
     {
       return CheckIndex(name) >= 0;
     }
+    */
+    
+    bool Used (std::string_view name) const
+    {
+      return CheckIndex(name) >= 0;
+    }
     
     /// Deletes symboltable
     inline void DeleteAll ()

libsrc/core/taskmanager.cpp

@@ -168,6 +168,12 @@ namespace ngcore
         trace = nullptr;
       }
     num_threads = 1;
+#ifdef USE_NUMA
+      for (int j = 0; j < num_nodes; j++)
+          numa_free (nodedata[j], sizeof(NodeData));
+#else
+      delete nodedata[0];
+#endif
   }
 
 #ifdef WIN32

libsrc/core/taskmanager.hpp

@@ -10,6 +10,7 @@
 #include <atomic>
 #include <functional>
 #include <list>
+#include <cmath>
 #include <ostream>
 #include <thread>
 
@@ -1010,7 +1011,7 @@ public:
     int num_nodes = numa_num_configured_nodes();
     size_t pagesize = numa_pagesize();
     
-    int npages = ceil ( double(s)*sizeof(T) / pagesize );
+    int npages = std::ceil ( double(s)*sizeof(T) / pagesize );
 
     // cout << "size = " << numa_size << endl;
     // cout << "npages = " << npages << endl;

libsrc/core/utils.hpp

@@ -182,7 +182,7 @@ namespace ngcore
   
   /// square element 
   template <class T>
-  NETGEN_INLINE T sqr (const T a)
+  NETGEN_INLINE constexpr T sqr (const T a)
   {
     return a * a; 
   }

libsrc/csg/python_csg.cpp

@@ -756,10 +756,8 @@ However, when r = 0, the top part becomes a point(tip) and meshing fails!
            {
              MeshingParameters mp;
              if(pars) mp = *pars;
-             {
-               py::gil_scoped_acquire aq;
              CreateMPfromKwargs(mp, kwargs);
-             }
+             py::gil_scoped_release gil_rel;
              auto mesh = make_shared<Mesh>();
              SetGlobalMesh (mesh);
              mesh->SetGeometry(geo);
@@ -770,8 +768,7 @@ However, when r = 0, the top part becomes a point(tip) and meshing fails!
                throw Exception("Meshing failed!");
              return mesh;
            }, py::arg("mp") = nullptr,
-       meshingparameter_description.c_str(),
+       meshingparameter_description.c_str())
-    py::call_guard<py::gil_scoped_release>())
     ;
 
   m.def("Save", FunctionPointer 

libsrc/geom2d/python_geom2d.cpp

@@ -401,10 +401,8 @@ NGCORE_API_EXPORT void ExportGeom2d(py::module &m)
 		{
                   MeshingParameters mp;
                   if(pars) mp = *pars;
-                  {
-                    py::gil_scoped_acquire aq;
                   CreateMPfromKwargs(mp, kwargs);
-                  }
+                  py::gil_scoped_release gil_release;
 		  auto mesh = make_shared<Mesh>();
                   mesh->SetGeometry(self);
                   SetGlobalMesh (mesh);
@@ -414,7 +412,6 @@ NGCORE_API_EXPORT void ExportGeom2d(py::module &m)
                     throw Exception("Meshing failed!");
 		  return mesh;
                 }, py::arg("mp") = nullopt,
-      py::call_guard<py::gil_scoped_release>(),
       meshingparameter_description.c_str())
     .def("_SetDomainTensorMeshing", &SplineGeometry2d::SetDomainTensorMeshing)
     ;
@@ -427,7 +424,8 @@ NGCORE_API_EXPORT void ExportGeom2d(py::module &m)
     .def(py::self-py::self)
     .def(py::self*py::self)
     .def(py::self+=py::self)
-    .def(py::self-=py::self)
+    // .def(py::self-=py::self) // false clange warning, see  https://github.com/pybind/pybind11/issues/1893
+    .def("__isub__", [](Solid2d& lhs, const Solid2d& rhs) { return lhs -= rhs; }, py::is_operator())
     .def(py::self*=py::self)
 
     .def("Mat", &Solid2d::Mat)
@@ -466,10 +464,8 @@ NGCORE_API_EXPORT void ExportGeom2d(py::module &m)
 		{
                   MeshingParameters mp;
                   if(pars) mp = *pars;
-                  {
-                    py::gil_scoped_acquire aq;
                     CreateMPfromKwargs(mp, kwargs);
-                  }
+                  py::gil_scoped_release gil_release;
 		  auto mesh = make_shared<Mesh>();
                   auto geo = self.GenerateSplineGeometry();
                   mesh->SetGeometry(geo);
@@ -480,7 +476,6 @@ NGCORE_API_EXPORT void ExportGeom2d(py::module &m)
                     throw Exception("Meshing failed!");
 		  return mesh;
                 }, py::arg("mp") = nullopt,
-      py::call_guard<py::gil_scoped_release>(),
       meshingparameter_description.c_str())
     ;
 

libsrc/include/nginterface_v2.hpp

@@ -160,6 +160,7 @@ namespace netgen
     FlatArray<T_FACE2> faces;    
     Ng_Facets facets;
     bool is_curved;
+    int8_t newest_vertex;
   };
 
   

libsrc/include/nginterface_v2_impl.hpp

@@ -191,6 +191,7 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<2> (size_t nr) const
       ret.facets.ptr = ret.edges.Data();
     }
   ret.is_curved = el.IsCurved();
+  ret.newest_vertex = el.NewestVertex();
   return ret;
 }
 
@@ -226,6 +227,7 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<3> (size_t nr) const
   ret.facets.ptr = ret.faces.Data();
 
   ret.is_curved = el.IsCurved();
+  ret.newest_vertex = el.NewestVertex();
   return ret;
 }
 

libsrc/meshing/adfront3.cpp

@@ -787,6 +787,54 @@ void AdFront3 :: SetStartFront (int /* baseelnp */)
     */
 }
 
+bool AdFront3 :: PointInsideGroup(const NgArray<PointIndex, PointIndex::BASE> &grouppindex,
+                                  const NgArray<MiniElement2d> &groupfaces) const
+{
+  for(auto pi : Range(points))
+    {
+      const auto& p = points[pi].P();
+      bool found = false;
+      for(const auto& f : groupfaces)
+        {
+          for(auto i : Range(3))
+            if(grouppindex.Get(f.PNum(i+1)) == pi)
+              {
+                found = true;
+                break;
+              }
+        }
+      if(found)
+        continue;
+
+      // "random" direction
+      Vec<3> dir = { 0.123871, 0.15432,-0.43989 };
+      DenseMatrix a(3), ainv(3);
+      Vector b(3), u(3);
+
+      int count = 0;
+      for(const auto& f : groupfaces)
+        {
+          const auto& p1 = points[grouppindex.Get(f.PNum(1))].P();
+          auto v1 = points[grouppindex.Get(f.PNum(2))].P() - p1;
+          auto v2 = points[grouppindex.Get(f.PNum(3))].P() - p1;
+          for(auto i : Range(3))
+            {
+              a(i,0) = v1[i];
+              a(i,1) = v2[i];
+              a(i,2) = -dir[i];
+              b(i) = p[i] - p1[i];
+            }
+          CalcInverse (a, ainv);
+          ainv.Mult (b, u);
+          if (u(0) >= 0 && u(1) >= 0 && u(0)+u(1) <= 1 &&
+	    u(2) > 0)
+	    count++;
+        }
+        if (count % 2 == 1)
+          return true;
+    }
+  return false;
+}
 
 bool AdFront3 :: Inside (const Point<3> & p) const
 {

libsrc/meshing/adfront3.hpp

@@ -262,6 +262,9 @@ public:
   void GetIntersectingFaces (const Point<3> & pmin, const Point<3> & pmax, 
 			     NgArray<int> & ifaces) const;
 
+  bool PointInsideGroup(const NgArray<PointIndex, PointIndex::BASE> &grouppindex,
+                        const NgArray<MiniElement2d>& groupfaces) const;
+
   ///
   void GetFaceBoundingBox (int i, Box3d & box) const;
 

libsrc/meshing/basegeom.cpp

@@ -8,25 +8,26 @@ namespace netgen
 {
   struct PointTree
   {
-      BoxTree<3> tree;
+      std::map<int, BoxTree<3>> tree;
+      Box<3> bounding_box;
 
-      PointTree( Box<3> bb ) : tree(bb) {}
+      PointTree( Box<3> bb ) : bounding_box(bb) {}
 
-      void Insert(Point<3> p, PointIndex n)
+      void Insert(Point<3> p, PointIndex n, int index)
       {
-          tree.Insert(p, p, n);
+        if(tree.count(index) == 0)
+          tree.emplace(index, bounding_box);
+        tree.at(index).Insert(p, p, n);
       }
 
-      PointIndex Find(Point<3> p) const
+      PointIndex Find(Point<3> p, int index) const
       {
           ArrayMem<int, 1> points;
-          tree.GetIntersecting(p, p, points);
+          tree.at(index).GetIntersecting(p, p, points);
           if(points.Size()==0)
               throw Exception("cannot find mapped point " + ToString(p));
           return points[0];
       }
-
-      double GetTolerance() { return tree.GetTolerance(); }
   };
 
   DLL_HEADER GeometryRegisterArray geometryregister;
@@ -266,7 +267,7 @@ namespace netgen
         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))
+              if(ident.trafo && e->IsMappedShape(*e_other, *ident.trafo, tol))
                 e->identifications.Append( {e, e_other, ident.trafo, ident.type, ident.name} );
 
       for(auto & e : edges)
@@ -278,9 +279,11 @@ namespace netgen
               GeometryVertex * pfrom[] = { &from.GetStartVertex(), &from.GetEndVertex() };
               GeometryVertex * pto[] = { &to.GetStartVertex(), &to.GetEndVertex() };
 
+              if(!ident.trafo) continue;
+
               // swap points of other edge if necessary
-              Point<3> p_from0 = ident.trafo(from.GetStartVertex().GetPoint());
+              Point<3> p_from0 = (*ident.trafo)(from.GetStartVertex().GetPoint());
-              Point<3> p_from1 = ident.trafo(from.GetEndVertex().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))
@@ -298,10 +301,7 @@ namespace netgen
       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;
 
@@ -315,16 +315,23 @@ namespace netgen
                   auto other = need_inverse ? ident.to : ident.from;
                   if(other->primary->nr < s->primary->nr)
                   {
-                      auto trafo = ident.trafo;
+                      s->primary = other->primary;
+                      if(ident.trafo)
+                      {
+                        auto trafo = *ident.trafo;
                         if(need_inverse)
                             trafo = trafo.CalcInverse();
-                      s->primary = other->primary;
+                        if(!s->primary_to_me)
-                      s->primary_to_me.Combine(trafo, other->primary_to_me);
+                          s->primary_to_me = Transformation<3>( Vec<3>{0., 0., 0.} );
+                        if(!other->primary_to_me)
+                          other->primary_to_me = Transformation<3>( Vec<3>{0., 0., 0.} );
+                        s->primary_to_me->Combine(trafo, *other->primary_to_me);
                         changed = true;
                     }
                   }
               }
             }
+          }
       };
 
       find_primary(vertices);
@@ -578,7 +585,6 @@ namespace netgen
     for(auto & vert : vertices)
       {
         auto pi = mesh.AddPoint(vert->GetPoint(), vert->properties.layer);
-        tree.Insert(mesh[pi], pi);
         vert2meshpt[vert->nr] = pi;
         mesh[pi].Singularity(vert->properties.hpref);
         mesh[pi].SetType(FIXEDPOINT);
@@ -610,7 +616,7 @@ namespace netgen
         endp = vert2meshpt[edge->GetEndVertex().nr];
 
         // ignore collapsed edges
-        if(startp == endp && edge->GetLength() < 1e-10 * bounding_box.Diam())
+        if(edge->IsDegenerated())
           continue;
 
         // ----------- Add Points to mesh and create segments -----
@@ -658,7 +664,9 @@ namespace netgen
             edge_params.SetSize(np-2);
             for(auto i : Range(np-2))
             {
-                edge_points[i] = trafo(mesh[pnums_primary[i+1]]);
+                edge_points[i] = mesh[pnums_primary[i+1]];
+                if(trafo)
+                  edge_points[i] = (*trafo)(edge_points[i]);
                 EdgePointGeomInfo gi;
                 edge->ProjectPoint(edge_points[i], &gi);
                 edge_params[i] = gi.dist;
@@ -689,7 +697,8 @@ namespace netgen
             {
               for(size_t i : std::vector{0UL, pnums_primary.Size()-1})
               {
-                auto p_mapped = trafo(mesh[pnums_primary[i]]);
+                auto p_mapped = mesh[pnums_primary[i]];
+                if(trafo) p_mapped = (*trafo)(p_mapped);
                 EdgePointGeomInfo gi;
                 edge->ProjectPoint(p_mapped, &gi);
                 params[i] = gi.dist;
@@ -707,7 +716,8 @@ namespace netgen
         for(auto i : Range(edge_points))
         {
             auto pi = mesh.AddPoint(edge_points[i], edge->properties.layer);
-            tree.Insert(mesh[pi], pi);
+            if(edge->identifications.Size())
+              tree.Insert(mesh[pi], pi, edge->nr);
             pnums[i+1] = pi;
         }
 
@@ -739,10 +749,16 @@ namespace netgen
           if(ident.from == edge.get())
           {
             auto & pnums = all_pnums[edge->nr];
+            if(pnums.Size() < 2) continue; // degenerated edge
             // 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]));
+                Point<3> p_other = mesh[pi];
+                if(ident.trafo)
+                  p_other = (*ident.trafo)(mesh[pi]);
+                else
+                  static_cast<GeometryEdge*>(ident.to)->ProjectPoint(p_other, nullptr);
+                auto pi_other = tree.Find(p_other, ident.to->nr);
                 identifications.Add(pi, pi_other, ident.name, ident.type);
             }
           }
@@ -851,7 +867,7 @@ namespace netgen
                 for(auto pi : s.PNums())
                     if(!is_point_in_tree[pi])
                     {
-                        tree.Insert(mesh[pi], pi);
+                        tree.Insert(mesh[pi], pi, -1);
                         is_point_in_tree[pi] = true;
                     }
 
@@ -860,7 +876,7 @@ namespace netgen
             constexpr int NOT_MAPPED = -1;
             mapped_edges = UNINITIALIZED;
 
-            Transformation<3> trafo;
+            optional<Transformation<3>> trafo;
 
             if(face.IsConnectingCloseSurfaces())
             {
@@ -902,8 +918,6 @@ namespace netgen
                         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]]));
                         auto gis = sel.GeomInfo();
                         for(auto i : Range(2))
                         {
@@ -911,6 +925,21 @@ namespace netgen
                             gis[i].v = s.epgeominfo[i].v;
                         }
 
+                        Point<3> p2 = mesh[s[1]];
+                        Point<3> p3 = mesh[s[0]];
+                        if(trafo)
+                        {
+                          p2 = (*trafo)(p2);
+                          p3 = (*trafo)(p3);
+                        }
+                        else 
+                        {
+                          edges[mapped_edges[edgenr]]->ProjectPoint(p2, nullptr);
+                          edges[mapped_edges[edgenr]]->ProjectPoint(p3, nullptr);
+                        }
+                        sel[2] = tree.Find(p2, -1);
+                        sel[3] = tree.Find(p3, -1);
+
                         // find mapped segment to set PointGeomInfo correctly
                         Segment s_other;
                         for(auto si_other : p2seg[sel[2]])
@@ -989,7 +1018,7 @@ namespace netgen
                         if(mesh[pi].Type() == SURFACEPOINT && pi_to_face[pi]==-1)
                         {
                             pi_to_face[pi] = face->nr;
-                            tree.Insert(mesh[pi], pi);
+                            tree.Insert(mesh[pi], pi, -1);
                             pi_of_face[face->nr].Append(pi);
                         }
                     }
@@ -1040,7 +1069,24 @@ namespace netgen
             auto pi = seg[i];
             if(!is_point_in_tree[pi])
             {
-              tree.Insert(trafo(mesh[pi]), pi);
+              auto p = mesh[pi];
+              if(trafo)
+                p = (*trafo)(p);
+              else
+                for(auto& edge: dst.edges)
+                  if (edge->primary->nr == seg.edgenr-1)
+                    {
+                      if (mesh[pi].Type() == FIXEDPOINT) {
+                        if((edge->GetStartVertex().GetPoint() - p).Length2() >\
+                           (edge->GetEndVertex().GetPoint() - p).Length2())
+                          p = edge->GetEndVertex().GetPoint();
+                        else
+                          p = edge->GetStartVertex().GetPoint();
+                      }
+                      else
+                        edge->ProjectPoint(p, nullptr);
+                    }
+              tree.Insert(p, pi, -1);
               is_point_in_tree[pi] = true;
             }
           }
@@ -1052,7 +1098,7 @@ namespace netgen
           {
             auto pi = seg[i];
             if(!pmap[pi].IsValid())
-                pmap[tree.Find(mesh[pi])] = pi;
+              pmap[tree.Find(mesh[pi], -1)] = pi;
 
             // store uv values (might be different values for same point in case of internal edges)
             double u = seg.epgeominfo[i].u;
@@ -1068,6 +1114,7 @@ namespace netgen
     }
 
     xbool do_invert = maybe;
+    if(!trafo) do_invert = true;
 
     // now insert mapped surface elements
     for(auto sei : mesh.SurfaceElements().Range())
@@ -1076,14 +1123,6 @@ namespace netgen
         if(sel.GetIndex() != src.nr+1)
           continue;
 
-        if(do_invert.IsMaybe())
-        {
-            auto n_src = src.GetNormal(mesh[sel[0]]);
-            auto n_dist = dst.GetNormal(trafo(mesh[sel[0]]));
-            Mat<3> normal_matrix;
-            CalcInverse(Trans(trafo.GetMatrix()), normal_matrix);
-            do_invert = (normal_matrix * n_src) * n_dist < 0.0;
-        }
         auto sel_new = sel;
         sel_new.SetIndex(dst.nr+1);
         for(auto i : Range(sel.PNums()))
@@ -1091,12 +1130,25 @@ namespace netgen
             auto pi = sel[i];
             if(!pmap[pi].IsValid())
               {
-                pmap[pi] = mesh.AddPoint(trafo(mesh[pi]), 1, SURFACEPOINT);
+                auto p = mesh[pi];
+                if(trafo)
+                  p = (*trafo)(p);
+                else
+                  dst.Project(p);
+                pmap[pi] = mesh.AddPoint(p, 1, SURFACEPOINT);
               }
               sel_new[i] = pmap[pi];
               mapto[{pi, dst.nr}] = pmap[pi];
               mapto[{pmap[pi], src.nr}] = pi;
           }
+        if(do_invert.IsMaybe())
+        {
+            auto n_src = src.GetNormal(mesh[sel[0]]);
+            auto n_dist = dst.GetNormal(mesh[sel_new[0]]);
+            Mat<3> normal_matrix;
+            CalcInverse(Trans(trafo->GetMatrix()), normal_matrix);
+            do_invert = (normal_matrix * n_src) * n_dist < 0.0;
+        }
         if(do_invert.IsTrue())
             sel_new.Invert();
 

libsrc/meshing/basegeom.hpp

@@ -54,7 +54,7 @@ namespace netgen
   {
     GeometryShape * from;
     GeometryShape * to;
-    Transformation<3> trafo;
+    optional<Transformation<3>> trafo;
     Identifications::ID_TYPE type;
     string name = "";
   };
@@ -67,7 +67,7 @@ namespace netgen
     ShapeProperties properties;
     Array<ShapeIdentification> identifications;
     GeometryShape * primary;
-    Transformation<3> primary_to_me;
+    optional<Transformation<3>> primary_to_me = nullopt;
 
     virtual ~GeometryShape() {}
     virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const;

libsrc/meshing/bisect.cpp

@@ -3,6 +3,9 @@
 #include "bisect.hpp"
 #include "validate.hpp"
 
+#include "meshing.hpp"  // quickfix for parallel 
+
+
 #define noDEBUG
 
 
@@ -31,6 +34,7 @@ namespace netgen
     // unsigned char faceedges[4];
     bool incorder;
     unsigned int order:6;
+    int8_t newest_vertex;
 
     MarkedTet() = default;
     /*
@@ -192,6 +196,7 @@ namespace netgen
 
     bool incorder;
     unsigned int order:6;
+    int8_t newest_vertex;
   };
   
   ostream & operator<< (ostream & ost, const MarkedTri & mt)
@@ -1255,6 +1260,8 @@ namespace netgen
     newtet1.marked = nm;
     newtet2.marked = nm;
 
+    newtet1.newest_vertex = oldtet.newest_vertex;
+
 #ifdef DEBUG
     *testout << "newtet1,before = " << newtet1 << endl;
     *testout << "newtet2,before = " << newtet2 << endl;
@@ -1264,6 +1271,7 @@ namespace netgen
       {
 	if (i == oldtet.tetedge1)
 	  {
+            newtet2.newest_vertex = i;
 	    newtet2.pnums[i] = newp;
 	    newtet2.faceedges[i] = oldtet.faceedges[i];  // inherited face
 	    newtet2.faceedges[vis1] = i;        // cut faces
@@ -1460,11 +1468,12 @@ namespace netgen
     newtri1.pnums[pe2] = newp;
     newtri1.pgeominfo[pe2] = newpgi;
     newtri1.markededge = pe2;
+    newtri1.newest_vertex = oldtri.newest_vertex;
 
     newtri2.pnums[pe1] = newp;
     newtri2.pgeominfo[pe1] = newpgi;
     newtri2.markededge = pe1;
-
+    newtri2.newest_vertex = pe1;
 
     newtri1.surfid = oldtri.surfid;
     newtri2.surfid = oldtri.surfid;
@@ -2792,6 +2801,20 @@ namespace netgen
     int np = mesh.GetNV();
     mesh.SetNP(np);
 
+
+#ifdef PARALLEL
+    if (mesh.GetCommunicator().Size() > 1)
+      {
+        mesh.GetParallelTopology().IdentifyVerticesAfterRefinement();
+        mesh.GetCommunicator().Barrier();
+        mesh.GetParallelTopology().EnumeratePointsGlobally();
+      }
+#endif
+
+
+
+
+    
     // int ne = mesh.GetNE();
     // int nse = mesh.GetNSE();
     // int i, j, l;
@@ -3698,6 +3721,7 @@ namespace netgen
             el.SetOrder (tet.order);
             for (int j = 0; j < 4; j++)
               el[j] = tet.pnums[j];
+            el.NewestVertex() = tet.newest_vertex;
             mesh.SetVolumeElement (ElementIndex(i), el);
           }
        });
@@ -3794,6 +3818,7 @@ namespace netgen
                 el[j] = trig.pnums[j];
                 el.GeomInfoPi(j+1) = trig.pgeominfo[j];
               }
+            el.NewestVertex() = trig.newest_vertex;
             mesh.SetSurfaceElement (SurfaceElementIndex(i), el);
           }
        });

libsrc/meshing/boundarylayer.cpp

@@ -919,6 +919,7 @@ void BoundaryLayerTool ::CreateNewFaceDescriptors() {
         FaceDescriptor new_fd(-1, isIn ? new_mat_nrs[i] : fd.DomainIn(),
                               isIn ? fd.DomainOut() : new_mat_nrs[i], -1);
         new_fd.SetBCProperty(new_si);
+        new_fd.SetSurfColour(fd.SurfColour());
         mesh.AddFaceDescriptor(new_fd);
         si_map[i] = new_si;
         moved_surfaces.SetBit(i);

libsrc/meshing/meshclass.cpp

@@ -901,6 +901,13 @@ namespace netgen
             outfile << " " << type;
           }
         outfile << "\n";
+        outfile << "identificationnames\n";
+        outfile << ident -> GetMaxNr() << "\n";
+        for (i = 1; i <= ident -> GetMaxNr(); i++)
+          {
+            string name = ident -> GetName(i);
+            outfile << ident->GetName(i) << "\n";
+          }
       }
 
     int cntmat = 0;
@@ -1208,6 +1215,8 @@ namespace netgen
 
     bool endmesh = false;
 
+    bool has_facedescriptors = false;
+    
 
     while (infile.good() && !endmesh)
       {
@@ -1226,6 +1235,7 @@ namespace netgen
 
         if (strcmp (str, "facedescriptors") == 0)
           {
+            has_facedescriptors = true;
             int nfd;
             infile >> nfd;
             for([[maybe_unused]] auto i : Range(nfd))
@@ -1448,6 +1458,17 @@ namespace netgen
                 ident -> SetType(i,Identifications::ID_TYPE(type));
               }
           }
+        if (strcmp (str, "identificationnames") == 0)
+          {
+            infile >> n;
+            PrintMessage (3, n, " identificationnames");
+            for (i = 1; i <= n; i++)
+              {
+                string name;
+                infile >> name;
+                ident -> SetName(i,name);
+              }
+          }
 
         if (strcmp (str, "materials") == 0)
           {
@@ -1610,7 +1631,11 @@ namespace netgen
 
                  surfnr--;
 
-                 if(surfnr > 0) 
+                 if(has_facedescriptors)
+                 {
+                    GetFaceDescriptor(i).SetSurfColour(surfcolour);
+                 }
+                 else if(surfnr > 0)
                  {
                     for(int facedesc = 1; facedesc <= cnt_facedesc; facedesc++)
                     {
@@ -1637,7 +1662,14 @@ namespace netgen
                     double transp;
                     infile >> surfnr >> transp;
                     surfnr--;
-                    if(surfnr > 0)
+                    if(has_facedescriptors)
+                    {
+                       auto& fd = GetFaceDescriptor(index);
+                       auto scol = fd.SurfColour();
+                       scol[3] = transp;
+                       fd.SetSurfColour(scol);
+                    }
+                    else if(surfnr > 0)
                       {
                         for(int facedesc = 1; facedesc <= cnt_facedesc; facedesc++)
                           {
@@ -1796,6 +1828,33 @@ namespace netgen
             archive & copy_el1d;
           }
 
+
+        // sending 0D elements
+        auto copy_el0d  (pointelements);
+        for (auto & el : copy_el0d)
+          {
+            auto & pi = el.pnum;
+            if (pi != PointIndex(PointIndex::INVALID))
+              pi = globnum[pi];
+          }
+        
+        if (comm.Rank() > 0)
+          comm.Send(copy_el0d, 0, 200);
+        else
+          {
+            Array<Element0d> el0di;
+            for (int j = 1; j < comm.Size(); j++)
+              {
+                comm.Recv(el0di, j, 200);
+                for (auto & el : el0di)
+                  copy_el0d += el;
+              }
+            archive & copy_el0d;
+          }
+
+
+
+        
         if (comm.Rank() == 0)
           {
             archive & facedecoding;
@@ -1825,6 +1884,7 @@ namespace netgen
     archive & surfelements;
     archive & volelements;
     archive & segments;
+    archive & pointelements;
     archive & facedecoding;
     archive & materials & bcnames & cd2names & cd3names;
     archive & numvertices;
@@ -6772,14 +6832,14 @@ namespace netgen
   //   }
   // #endif
 
-  int Mesh::IdentifyPeriodicBoundaries(const string &s1,
+  int Mesh::IdentifyPeriodicBoundaries(const string& id_name,
-                                       const string &s2,
+                                       const string &s1,
                                        const Transformation<3> &mapping,
                                        double pointTolerance)
   {
-    auto nr = ident->GetMaxNr() + 1;
+    auto nr = ident->GetNr(id_name);
     ident->SetType(nr, Identifications::PERIODIC);
-    double lami[4];
+    // double lami[4];
     set<int> identified_points;
     if(pointTolerance < 0.)
       {
@@ -6787,44 +6847,39 @@ namespace netgen
         GetBox(pmin, pmax);
         pointTolerance = 1e-8 * (pmax-pmin).Length();
       }
-    for(const auto& se : surfelements)
+    size_t nse = GetDimension() == 3 ? surfelements.Size() : segments.Size();
+    for(auto sei : Range(nse))
       {
-        if(GetBCName(se.index-1) != s1)
+        auto name = GetDimension() == 3 ? GetBCName(surfelements[sei].index-1) :
+          GetBCName(segments[sei].edgenr-1);
+        if(name != s1)
           continue;
 
-        for(const auto& pi : se.PNums())
+        const auto& pnums = GetDimension() == 3 ? surfelements[sei].PNums() :
+          segments[sei].PNums();
+        for(const auto& pi : pnums)
           {
             if(identified_points.find(pi) != identified_points.end())
               continue;
             auto pt = (*this)[pi];
             auto mapped_pt = mapping(pt);
-            auto other_nr = GetElementOfPoint(mapped_pt, lami, true);
+            bool found = false;
-            int index = -1;
+            for(auto other_pi : Range(points))
-            if(other_nr != 0)
               {
-                auto other_el = VolumeElement(other_nr);
+                if((mapped_pt - (*this)[other_pi]).Length() < pointTolerance)
-                for(auto i : Range(other_el.PNums().Size()))
-                  if((mapped_pt - (*this)[other_el.PNums()[i]]).Length() < pointTolerance)
                   {
-                      index = i;
+                    identified_points.insert(pi);
+                    ident->Add(pi, other_pi, nr);
+                    found = true;
                     break;
                   }
-                if(index == -1)
+              }
+            if(!found)
               {
                 cout << "point coordinates = " << pt << endl;
                 cout << "mapped coordinates = " << mapped_pt << endl;
                 throw Exception("Did not find mapped point with nr " + ToString(pi) + ", are you sure your mesh is periodic?");
               }
-                auto other_pi = other_el.PNums()[index];
-                identified_points.insert(pi);
-                ident->Add(pi, other_pi, nr);
-              }
-            else
-              {
-                cout << "point coordinates = " << pt << endl;
-                cout << "mapped coordinates = " << mapped_pt << endl;
-                throw Exception("Mapped point with nr " + ToString(pi) + " is outside of mesh, are you sure your mesh is periodic?");
-              }
           }
       }
     return nr;
@@ -7233,7 +7288,7 @@ namespace netgen
 
     for(auto dom : Range(ndomains))
     {
-      if(regex_match(mesh.GetMaterial(dom), regex_domains))
+      if(regex_match(mesh.GetMaterial(dom+1), regex_domains))
         keep_domain.SetBit(dom);
     }
 
@@ -7244,7 +7299,7 @@ namespace netgen
         keep_face.SetBit(fd.BCProperty());
     }
 
-    auto filter_elements = [&mesh, &keep_point](auto & elements, auto & keep_region)
+    auto filter_elements = [&keep_point](auto & elements, auto & keep_region)
     {
       for(auto & el : elements)
       {

libsrc/meshing/meshclass.hpp

@@ -776,8 +776,8 @@ namespace netgen
     { return facedecoding[i-1]; }      
     // { return facedecoding.Elem(i); }
 
-    int IdentifyPeriodicBoundaries(const string& s1,
+    int IdentifyPeriodicBoundaries(const string& id_name,
-                                   const string& s2,
+                                   const string& s1,
                                    const Transformation<3>& mapping,
                                    double pointTolerance);
 

libsrc/meshing/meshfunc.cpp

@@ -685,7 +685,7 @@ namespace netgen
     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))
+    for([[maybe_unused]] auto i : Range(mp.optsteps3d))
       {
         auto [total_badness, max_badness, bad_els] = optmesh.UpdateBadness();
         if(bad_els==0) break;

libsrc/meshing/meshing3.cpp

@@ -372,7 +372,8 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 	      onlytri = 0;
 	  
 	  if (onlytri && groupfaces.Size() <= 20 + 2*stat.qualclass &&
-	      FindInnerPoint (grouppoints, groupfaces, inp))
+	      FindInnerPoint (grouppoints, groupfaces, inp) &&
+              !adfront->PointInsideGroup(grouppindex, groupfaces))
 	    {
 	      (*testout) << "inner point found" << endl;
 

libsrc/meshing/meshtype.cpp

@@ -144,7 +144,40 @@ namespace netgen
 #endif
 
 
+#ifdef PARALLEL
+  NG_MPI_Datatype Element0d :: MyGetMPIType()
+  {
+    static NG_MPI_Datatype type = NG_MPI_DATATYPE_NULL;
+    static NG_MPI_Datatype htype = NG_MPI_DATATYPE_NULL;
+    if (type == NG_MPI_DATATYPE_NULL)
+      {
+	Element0d hel;
+	int blocklen[] = { 1, 1 };
+	NG_MPI_Aint displ[] =
+          { (char*)&hel.pnum - (char*)&hel,
+            (char*)&hel.index - (char*)&hel,
+          };
+	NG_MPI_Datatype types[] = {
+          GetMPIType(hel.pnum), GetMPIType(hel.index)
+        };
+	NG_MPI_Type_create_struct (2, blocklen, displ, types, &htype);
+	NG_MPI_Type_commit ( &htype );
+	NG_MPI_Aint lb, ext;
+	NG_MPI_Type_get_extent (htype, &lb, &ext);
+	// *testout << "lb = " << lb << endl;
+	// *testout << "ext = " << ext << endl;
+	ext = sizeof (Element0d);
+	NG_MPI_Type_create_resized (htype, lb, ext, &type);
+	NG_MPI_Type_commit ( &type );
+      }
+    return type;
+  }
+#endif
 
+ void Element0d :: DoArchive (Archive & ar)
+ {
+   ar & pnum & index;
+ }
 
   Segment :: Segment() 
     : is_curved(false)
@@ -2693,6 +2726,8 @@ namespace netgen
           for (auto & t : type)
             ar & (unsigned char&)(t);
         }
+      if (ar.GetVersion("netgen") > "v6.2.2404-66")
+        ar & names;
     }    
 
   

libsrc/meshing/meshtype.hpp

@@ -425,6 +425,7 @@ namespace netgen
     // control whether it is visible or not
     bool visible:1;  // element visible
     bool is_curved;   // element is (high order) curved
+    int8_t newest_vertex = -1; // from refinement via bisection
     /// order for hp-FEM
     unsigned int orderx:6;
     unsigned int ordery:6;
@@ -562,6 +563,9 @@ namespace netgen
     ///
     const PointGeomInfo & GeomInfoPiMod (int i) const { return geominfo[(i-1) % np]; }
 
+    auto & NewestVertex() { return newest_vertex; }
+    auto NewestVertex() const { return newest_vertex; }
+
     void DoArchive (Archive & ar)
     {
       short _np, _typ;
@@ -731,6 +735,7 @@ namespace netgen
     ELEMENT_TYPE typ;
     /// number of points (4..tet, 5..pyramid, 6..prism, 8..hex, 10..quad tet, 12..quad prism)
     int8_t np;
+    int8_t newest_vertex = -1; // from refinement via bisection
     
     /// sub-domain index
     int index;
@@ -856,6 +861,9 @@ namespace netgen
     ///
     const PointIndex & PNumMod (int i) const { return pnum[(i-1) % np]; }
 
+    auto & NewestVertex() { return newest_vertex; }
+    auto NewestVertex() const { return newest_vertex; }
+
     void DoArchive (Archive & ar)
     {
       short _np, _typ;
@@ -1156,6 +1164,12 @@ namespace netgen
     Element0d () = default;
     Element0d (PointIndex _pnum, int _index)
       : pnum(_pnum), index(_index) { ; }
+
+#ifdef PARALLEL
+    static NG_MPI_Datatype MyGetMPIType();
+#endif
+    
+    void DoArchive (Archive & ar);
   };
 
   ostream & operator<<(ostream  & s, const Element0d & el);
@@ -1629,6 +1643,19 @@ namespace netgen
          names.Append(name);
       return names.Pos(name)+1;
     }
+    string GetName(int nr) const
+    {
+      if (nr <= names.Size())
+        return names[nr - 1];
+      else
+        return "";
+    }
+    void SetName(int nr, string name)
+    {
+      while(names.Size() < nr)
+        names.Append("");
+      names[nr-1] = name;
+    }
 
     /// remove secondorder
     void SetMaxPointNr (int maxpnum);
@@ -1666,6 +1693,9 @@ namespace ngcore
   template <> struct MPI_typetrait<netgen::Segment> {
     static NG_MPI_Datatype MPIType ()  { return netgen::Segment::MyGetMPIType(); }
   };
+  template <> struct MPI_typetrait<netgen::Element0d> {
+    static NG_MPI_Datatype MPIType ()  { return netgen::Element0d::MyGetMPIType(); }
+  };
 
 }
 #endif

libsrc/meshing/paralleltop.cpp

@@ -292,7 +292,7 @@ namespace netgen
 
 
 
-
+  /*
 
   void ParallelMeshTopology :: UpdateCoarseGridGlobal ()
   {
@@ -387,6 +387,7 @@ namespace netgen
     
     is_updated = true;
   }
+  */
   
   
   void ParallelMeshTopology :: IdentifyVerticesAfterRefinement()

libsrc/meshing/paralleltop.hpp

@@ -36,8 +36,8 @@ namespace netgen
 
     
     void UpdateCoarseGrid();
-    [[deprecated("should not need it anymore")]]                    
+    // [[deprecated("should not need it anymore")]]                    
-    void UpdateCoarseGridGlobal();
+    // void UpdateCoarseGridGlobal();
     void IdentifyVerticesAfterRefinement();
     void EnumeratePointsGlobally ();
         

libsrc/meshing/python_mesh.cpp

@@ -1,4 +1,3 @@
-#include "pybind11/pytypes.h"
 #ifdef NG_PYTHON
 
 #include <regex>
@@ -180,6 +179,34 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
 
   py::implicitly_convertible<py::tuple, Vec<3>>();
 
+  py::class_<Mat<3,3>>(m, "Mat33")
+    .def(py::init([](py::tuple m)
+    {
+      if(m.size() != 9)
+        throw std::length_error("Invalid dimension of input array!");
+      Mat<3,3> mat;
+      for(int i = 0; i < 3; i++)
+        for(int j = 0; j < 3; j++)
+          mat(i,j) = m[i*3+j].cast<double>();
+      return mat;
+    }))
+    .def("__getitem__", [](Mat<3,3>& mat, py::tuple index)
+    {
+      if(index.size() != 2)
+        throw std::length_error("Invalid dimension of input array!");
+      return mat(index[0].cast<int>(), index[1].cast<int>());
+    })
+    .def("__setitem__", [](Mat<3,3>& mat, py::tuple index, double val)
+    {
+      if(index.size() != 2)
+        throw std::length_error("Invalid dimension of input array!");
+      mat(index[0].cast<int>(), index[1].cast<int>()) = val;
+    })
+    .def("__str__", &ToString<Mat<3,3>>)
+    ;
+
+  py::implicitly_convertible<py::tuple, Mat<3,3>>();
+
   m.def ("Vec", FunctionPointer
            ([] (double x, double y, double z) { return global_trafo(Vec<3>(x,y,z)); }));
   m.def("Vec", [](py::array_t<double> np_array)
@@ -204,6 +231,11 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     .def("__mul__", [](Transformation<3> a, Transformation<3> b)->Transformation<3>
          { Transformation<3> res; res.Combine(a,b); return res; })
     .def("__call__", [] (Transformation<3> trafo, Point<3> p) { return trafo(p); })
+    .def_property("mat", &Transformation<3>::GetMatrix,
+                  [](Transformation<3>& self, const Mat<3,3>& mat)
+                  {
+                    self.GetMatrix() = mat;
+                  })
     ;
 
   m.def ("GetTransformation", [] () { return global_trafo; });
@@ -536,14 +568,24 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
                                      li.append (py::cast(self.surfnr2));
                                      return li;
                                    }))
-    .def_property_readonly("index", FunctionPointer([](const Segment &self) -> size_t
+    .def_property("index",
+                  [](const Segment &self) -> size_t
                   {
                     return self.si;
-		  }))
+                  },
-    .def_property_readonly("edgenr", FunctionPointer([](const Segment & self) -> size_t
+                  [](Segment& self, int index)
+                  {
+                    self.si = index;
+                  })
+    .def_property("edgenr",
+                  [](const Segment & self) -> size_t
                   {
                     return self.edgenr;
-						     }))
+                  },
+                  [](Segment& self, int edgenr)
+                  {
+                    self.edgenr = edgenr;
+                  })
     .def_property("singular",
                   [](const Segment & seg) { return seg.singedge_left; },
                   [](Segment & seg, double sing) { seg.singedge_left = sing; seg.singedge_right=sing; })
@@ -1215,7 +1257,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
            py::arg("identnr"),
            py::arg("type")=Identifications::PERIODIC)
     .def("IdentifyPeriodicBoundaries", &Mesh::IdentifyPeriodicBoundaries,
-         py::arg("face1"), py::arg("face2"), py::arg("mapping"), py::arg("point_tolerance") = -1.)
+         py::arg("identification_name"), py::arg("face1"), py::arg("mapping"),
+py::arg("point_tolerance") = -1.)
     .def("GetNrIdentifications", [](Mesh& self)
                                  {
                                    return self.GetIdentifications().GetMaxNr();
@@ -1234,15 +1277,12 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
            {
              MeshingParameters mp;
              if(pars) mp = *pars;
-             {
-               py::gil_scoped_acquire acquire;
              CreateMPfromKwargs(mp, kwargs);
-             }
+             py::gil_scoped_release gil_release;
              MeshVolume (mp, self);
              OptimizeVolume (mp, self);
            }, py::arg("mp")=nullptr,
-          meshingparameter_description.c_str(),
+          meshingparameter_description.c_str())
-          py::call_guard<py::gil_scoped_release>())
 
     .def ("OptimizeVolumeMesh", [](Mesh & self, MeshingParameters* pars)
           {

libsrc/occ/occ_utils.cpp

@@ -61,7 +61,7 @@ namespace netgen
     Standard_Integer BuildTriangulation( const TopoDS_Shape & shape )
     {
        BRepTools::Clean (shape);
-       double deflection = 0.01;
+       // double deflection = 0.01;
 
        // https://dev.opencascade.org/doc/overview/html/occt_user_guides__mesh.html
        // from Standard_Boolean meshing_imeshtools_parameters()

libsrc/occ/occ_utils.hpp

@@ -65,15 +65,14 @@ namespace netgen
 
     DLL_HEADER Box<3> GetBoundingBox( const TopoDS_Shape & shape );
 
-    class OCCIdentification
+    struct OCCIdentification
     {
-    public:
       TopoDS_Shape from;
       TopoDS_Shape to;
-      Transformation<3> trafo;
+      optional<Transformation<3>> trafo = nullopt;
       string name;
       Identifications::ID_TYPE type;
-      bool opposite_direction;
+      bool opposite_direction = false;
     };
 
     Standard_Integer BuildTriangulation( const TopoDS_Shape & shape );

libsrc/occ/occgeom.cpp

@@ -14,6 +14,7 @@
 #include "occgeom.hpp"
 #include "Partition_Spliter.hxx"
 
+#include <BinTools.hxx>
 #include <BOPAlgo_Builder.hxx>
 #include <BRepBndLib.hxx>
 #include <BRepBuilderAPI_Copy.hxx>
@@ -1637,9 +1638,13 @@ namespace netgen
       Standard_Boolean result = BRepTools::Read(occgeo->shape, filename.string().c_str(), aBuilder);
 
       if(!result)
+      {
+        result = BinTools::Read(occgeo->shape, filename.string().c_str());
+        if (!result)
           {
             delete occgeo;
-         return NULL;
+            throw Exception("Could not read BREP file " + filename.string());
+          }
       }
 
       occgeo->changed = 1;
@@ -1767,7 +1772,18 @@ namespace netgen
                   id_from = shape_map.FindIndex(id.from)-1;
                   id_to = shape_map.FindIndex(id.to)-1;
                 }
-                ar & id_from & id_to & id.trafo & id.name;
+                ar & id_from & id_to;
+
+                // trafo is now optional -> special treatment necessary for backward compatibility
+                if(ar.Output() || netgen_version >= "v6.2.2403-34-g571cbbe4")
+                  ar & id.trafo;
+                else
+                {
+                  auto trafo = Transformation<3>();
+                  ar & trafo;
+                  id.trafo = trafo;
+                }
+                ar & id.name;
                 if(ar.Input())
                 {
                     id.from = shape_list[id_from];
@@ -2365,10 +2381,12 @@ namespace netgen
               Array<Handle(StepRepr_RepresentationItem)> items;
               items.Append(MakeReal(ident.from == shape ? 1 : 0));
               items.Append(to);
-              auto & m = ident.trafo.GetMatrix();
+              Transformation<3> trafo;
+              if(ident.trafo) trafo = *ident.trafo;
+              auto & m = trafo.GetMatrix();
               for(auto i : Range(9))
                   items.Append(MakeReal(m(i)));
-              auto & v = ident.trafo.GetVector();
+              auto & v = trafo.GetVector();
               for(auto i : Range(3))
                   items.Append(MakeReal(v(i)));
               items.Append(MakeInt(ident.type));
@@ -2407,12 +2425,15 @@ namespace netgen
                   ident.to = shape_origin;
                 }
 
-              auto & m = ident.trafo.GetMatrix();
+              Transformation<3> trafo;
+              auto & m = trafo.GetMatrix();
               for(auto i : Range(9))
                   m(i) = ReadReal(id_item->ItemElementValue(3+i));
-              auto & v = ident.trafo.GetVector();
+              auto & v = trafo.GetVector();
               for(auto i : Range(3))
                   v(i) = ReadReal(id_item->ItemElementValue(12+i));
+              if(FlatVector(9, &trafo.GetMatrix()(0,0)).L2Norm() != .0 && trafo.GetVector().Length2() != .0)
+                  ident.trafo = trafo;
               ident.type = Identifications::ID_TYPE(ReadInt(id_item->ItemElementValue(15)));
               result.push_back(ident);
           }

libsrc/occ/occgeom.hpp

@@ -516,11 +516,13 @@ namespace netgen
                       if(from.IsSame(from_mapped) && to.IsSame(to_mapped))
                         continue;
   
-                      Transformation<3> trafo_mapped = ident.trafo;
+                      if(!ident.trafo) continue;
+                      Transformation<3> trafo_mapped = *ident.trafo;
+
                       if(trafo)
                       {
                           Transformation<3> trafo_temp;
-                          trafo_temp.Combine(ident.trafo, trafo_inv);
+                          trafo_temp.Combine(*ident.trafo, trafo_inv);
                           trafo_mapped.Combine(*trafo, trafo_temp);
                       }
   

libsrc/occ/python_occ.cpp

@@ -86,7 +86,7 @@ DLL_HEADER void ExportNgOCC(py::module &m)
     try {
       if(p) std::rethrow_exception(p);
     } catch (const Standard_Failure& e) {
-      exc((string(e.DynamicType()->Name()) + ": " + e.GetMessageString()).c_str());
+      py::set_error(PyExc_RuntimeError, (string(e.DynamicType()->Name()) + ": " + e.GetMessageString()).c_str());
     }
   });
   
@@ -243,8 +243,6 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                          {
                            MeshingParameters mp;
                            OCCParameters occparam;
-                           {
-                             py::gil_scoped_acquire aq;
                            if(pars)
                              {
                                auto mp_kwargs = CreateDictFromFlags(pars->geometrySpecificParameters);
@@ -253,7 +251,7 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                              }
                            CreateOCCParametersFromKwargs(occparam, kwargs);
                            CreateMPfromKwargs(mp, kwargs);
-                           }
+                           py::gil_scoped_release gil_release;
                            geo->SetOCCParameters(occparam);
                            if(!mesh)
                              mesh = make_shared<Mesh>();
@@ -279,7 +277,7 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                              }
                            return mesh;
                          }, py::arg("mp") = nullptr, py::arg("comm")=NgMPI_Comm{},
-         py::arg("mesh")=nullptr, py::call_guard<py::gil_scoped_release>(),
+         py::arg("mesh")=nullptr,
          (meshingparameter_description + occparameter_description).c_str())
     .def_property_readonly("shape", [](const OCCGeometry & self) { return self.GetShape(); })
     ;

libsrc/occ/python_occ_basic.cpp

@@ -25,7 +25,7 @@ DLL_HEADER void ExportNgOCCBasic(py::module &m)
     .def(py::init([] (py::tuple pnt)
                   {
                     if (py::len(pnt) != 3)
-                      throw Exception("need 3-tuple to create gp_Pnt");
+                      throw std::length_error("need 3-tuple to create gp_Pnt");
                     
                     return gp_Pnt(py::cast<double>(pnt[0]),
                                   py::cast<double>(pnt[1]),

libsrc/occ/python_occ_shapes.cpp

@@ -780,7 +780,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
            PropagateProperties(builder, shape, occ2ng(trafo));
            return builder.Shape();
          }, py::arg("axes"),
-         "copy shape, and mirror over plane defined by 'axes'")
+         "copy shape, and mirror over XY - plane defined by 'axes'")
     
     .def("Mirror", [] (const TopoDS_Shape & shape, const gp_Ax1 & ax)
          {
@@ -790,7 +790,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
            PropagateProperties(builder, shape, occ2ng(trafo));
            return builder.Shape();
          }, py::arg("axes"),
-         "copy shape, and mirror around axis 'axis'")
+         "copy shape, and rotate by 180 deg around axis 'axis'")
     
     .def("Scale", [](const TopoDS_Shape & shape, const gp_Pnt p, double s)
          {
@@ -1150,6 +1150,59 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
          py::arg("removeIntersectingEdges") = false,
          "makes shell-like solid from faces")
 
+    .def("Offset", [](const TopoDS_Shape & shape, 
+                      double offset, double tol, bool intersection,
+                      string joinT, bool removeIntEdges, optional<string> identification_name) {
+           BRepOffsetAPI_MakeOffsetShape maker;
+           GeomAbs_JoinType joinType;
+           if(joinT == "arc")
+             joinType = GeomAbs_Arc;
+           else if(joinT == "intersection")
+             joinType = GeomAbs_Intersection;
+           else if(joinT == "tangent")
+            joinType = GeomAbs_Tangent;
+           else
+             throw Exception("Only joinTypes 'arc', 'intersection' and 'tangent' exist!");
+           
+           maker.PerformByJoin(shape, offset, tol,
+                               BRepOffset_Skin, intersection,
+                               false, joinType, removeIntEdges);
+
+           // PropagateProperties (maker, shape);
+           for (auto typ : { TopAbs_FACE,  TopAbs_EDGE, TopAbs_VERTEX })
+             for (TopExp_Explorer e(shape, typ); e.More(); e.Next())
+               {
+                 auto s = e.Current();
+                 auto prop = OCCGeometry::GetProperties(s);
+                 for (auto mods : maker.Generated(s))
+                   {
+                     if(OCCGeometry::HaveProperties(s))
+                       {
+                         auto & new_props = OCCGeometry::GetProperties(mods);
+                         new_props.Merge(prop);
+                         if (prop.name) new_props.name = string("offset_")+(*prop.name);
+                       }
+                     if(identification_name)
+                       {
+                         OCCIdentification ident;
+                         ident.from = s;
+                         ident.to = mods;
+                         ident.name = *identification_name;
+                         ident.type = Identifications::CLOSESURFACES;
+                         OCCGeometry::GetIdentifications(s).push_back(ident);
+                       }
+                   }
+               }
+           
+           return maker.Shape();
+       }, py::arg("offset"), py::arg("tol"),
+         py::arg("intersection") = false,py::arg("joinType")="arc",
+         py::arg("removeIntersectingEdges") = false,
+         py::arg("identification_name") = nullopt,
+         "makes shell-like solid from faces")
+
+
+    
     .def("MakeTriangulation", [](const TopoDS_Shape & shape)
          {
            BuildTriangulation(shape);
@@ -1209,7 +1262,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
          })
     .def("_webgui_data", [](const TopoDS_Shape & shape)
          {
-           auto status = BuildTriangulation(shape);
+           [[maybe_unused]] auto status = BuildTriangulation(shape);
            // cout << "status = " << aStatus << endl;
            
            std::vector<double> p[3];
@@ -1479,6 +1532,24 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
               throw NgException("error in wire builder: "+errstr.str());
             }
         }))
+    .def("Offset", [](const TopoDS_Wire & wire, const TopoDS_Face & face, double dist,
+                      string joinT, bool openresult)
+    {
+      GeomAbs_JoinType joinType;
+      if(joinT == "arc")
+        joinType = GeomAbs_Arc;
+      else if(joinT == "intersection")
+        joinType = GeomAbs_Intersection;
+      else if(joinT == "tangent")
+        joinType = GeomAbs_Tangent;
+      else
+        throw Exception("Only joinTypes 'arc', 'tangent', and 'intersection' exist!");
+      BRepOffsetAPI_MakeOffset builder(face, joinType, openresult);
+      builder.AddWire(wire);
+      builder.Perform(dist);
+      auto shape = builder.Shape();    
+      return shape;
+    })
     ;
 
   py::class_<TopoDS_Face, TopoDS_Shape> (m, "Face")

libsrc/stlgeom/python_stl.cpp

@@ -188,7 +188,6 @@ NGCORE_API_EXPORT void ExportSTL(py::module & m)
                          {
                            MeshingParameters mp;
                            STLParameters stlparam;
-                           { py::gil_scoped_acquire aq;
                            if(pars)
                            {
                              auto mp_flags = pars->geometrySpecificParameters;
@@ -198,7 +197,7 @@ NGCORE_API_EXPORT void ExportSTL(py::module & m)
                            }
                            CreateSTLParametersFromKwargs(stlparam, kwargs);
                            CreateMPfromKwargs(mp, kwargs); // this will throw if any kwargs are not passed
-                           }
+                           py::gil_scoped_release gil_release;
                            if(!mesh)
                              {
                                mesh = make_shared<Mesh>();
@@ -215,7 +214,6 @@ NGCORE_API_EXPORT void ExportSTL(py::module & m)
 
                            return mesh;
                          }, py::arg("mp") = nullptr, py::arg("mesh") = nullptr,
-      py::call_guard<py::gil_scoped_release>(),
          (meshingparameter_description + stlparameter_description).c_str())
     .def("Draw", FunctionPointer
          ([] (shared_ptr<STLGeometry> self)

libsrc/visualization/mvdraw.hpp

@@ -145,16 +145,16 @@ namespace netgen
     int filledtimestamp = -1;
     int linetimestamp = -1;
     int edgetimestamp = -1;
-    int pointnumbertimestamp = -1;
+    // int pointnumbertimestamp = -1;
 
     int tettimestamp = -1;
     int prismtimestamp = -1;
     int pyramidtimestamp = -1;
     int hextimestamp = -1;
 
-    int badeltimestamp = -1;
+    // int badeltimestamp = -1;
-    int identifiedtimestamp = -1;
+    // int identifiedtimestamp = -1;
-    int domainsurftimestamp = -1;
+    // int domainsurftimestamp = -1;
 
     struct {
       unsigned texture = -1;

ng/ngappinit.cpp

@@ -13,6 +13,10 @@
 
 #include "../libsrc/interface/writeuser.hpp"
 
+#ifdef NETGEN_PYTHON
+#include <core/python_ngcore.hpp>
+#endif
+
 namespace netgen
 {
   DLL_HEADER extern Flags parameters;
@@ -251,6 +255,10 @@ int main(int argc, char ** argv)
       // start event-loop
       Tk_MainLoop();
       Tcl_DeleteInterp (myinterp); 
+#ifdef NETGEN_PYTHON
+      py::gil_scoped_acquire ensure_gil;
+#endif
+
       Tcl_Exit(0);
     }
 

ng/ngpkg.cpp

@@ -2078,7 +2078,6 @@ namespace netgen
       return TCL_ERROR;
     const char * filename = Tcl_GetString(argv[2]);
 
-    int len = strlen(filename);
     int w = Togl_PixelScale(togl)*Togl_Width (togl);
     int h = Togl_PixelScale(togl)*Togl_Height (togl);
 
@@ -2088,6 +2087,7 @@ namespace netgen
     glReadPixels (0, 0, w, h, GL_RGB, GL_UNSIGNED_BYTE, &buffer[0]);
 
 #ifdef JPEGLIB
+    int len = strlen(filename);
     if (strcmp ("jpg", filename+len-3) == 0)
       {
         cout << "Snapshot to file '" << filename << "'" << endl;
@@ -2130,19 +2130,15 @@ namespace netgen
 #endif // JPEGLIB
     {
         string command;
-        string filename2;
+        std::filesystem::path filepath(filename);
 
-        filename2 = filename;
+        bool need_conversion = filepath.extension() != ".ppm";
+        if (need_conversion)
+          filepath += ".ppm";
 
-        if(filename2.substr(len-3) != ".ppm")
+        cout << IM(3) << "Snapshot to file '" << filepath.string() << endl;
-            filename2 += ".ppm";
 
-        cout << "Snapshot to file '" << filename << endl;
+        ofstream outfile(filepath);
-
-        // int w = Togl_Width (togl);
-        // int h = Togl_Height (togl);
-
-        ofstream outfile(filename2);
         outfile << "P6" << endl
           << "# CREATOR: Netgen" << endl
           << w << " " << h << endl
@@ -2153,12 +2149,10 @@ namespace netgen
                     outfile.put (buffer[k+3*j+3*w*(h-i-1)]);
         outfile << flush;
 
-        if (filename2 == string(filename))
+        if (need_conversion)
-            return TCL_OK;
-        else
         {
           // convert image file (Unix/Linux only):
-          command = string("convert -quality 100 ") + filename2 + " " + filename;
+          command = string("convert -quality 100 ") + filepath.string() + " " + filename;
 
           int err = system(command.c_str());
           if (err != 0)
@@ -2167,18 +2161,12 @@ namespace netgen
               return TCL_ERROR;
           }
 
-          command  = string("rm ") + filename2;
+          std::filesystem::remove(filepath);
-          err = system(command.c_str());
+        }
 
-          if (err != 0)
-          {
-              Tcl_SetResult (Togl_Interp(togl), (char*)"Cannot delete temporary file", TCL_VOLATILE);
-              return TCL_ERROR;
-          }
         return TCL_OK;
     }
   }
-  }
 
 
 

python/__init__.py

@@ -1,10 +1,65 @@
 import os
 import sys
+from pathlib import Path
 
 from . import config
 _netgen_bin_dir=os.path.realpath(os.path.join(os.path.dirname(__file__),'..',config.NETGEN_PYTHON_RPATH_BIN))
 _netgen_lib_dir=os.path.realpath(os.path.join(os.path.dirname(__file__),'..',config.NETGEN_PYTHON_RPATH))
 
+def load_occ_libs():
+    try:
+        try:
+            import importlib.metadata as metadata
+        except ImportError:
+            import importlib_metadata as metadata
+        import ctypes
+        metadata.metadata('netgen-occt')
+        lib_names = [
+            "TKOffset",
+            "TKFillet",
+            "TKDEIGES",
+            "TKBool",
+            "TKDESTEP",
+            "TKXSBase",
+            "TKDESTL",
+            "TKXCAF",
+            "TKVCAF",
+            "TKCAF",
+            "TKBO",
+            "TKPrim",
+            "TKLCAF",
+            "TKCDF",
+            "TKV3d",
+            "TKHLR",
+            "TKMesh",
+            "TKService",
+            "TKShHealing",
+            "TKTopAlgo",
+            "TKGeomAlgo",
+            "TKBRep",
+            "TKGeomBase",
+            "TKG3d",
+            "TKG2d",
+            "TKMath",
+            "TKDE",
+            "TKernel",
+        ]
+        lib_names.reverse()
+        lib_paths = {}
+        for f in metadata.files('netgen-occt'):
+            if f.match('*libTK*') or f.match("*.dll"):
+                p = f.locate()
+                name = p.name.split('.')[0].lower().replace("lib","")
+                lib_paths[name] = str(p)
+        for lib_name in lib_names:
+            p = lib_paths[lib_name.lower()]
+            ctypes.CDLL(p, mode=ctypes.RTLD_GLOBAL)
+
+    except metadata.PackageNotFoundError:
+        pass
+
+load_occ_libs()
+
 __diagnostics_template = """
 Netgen diagnostics:
     sys.platform:          {sys.platform}

python/webgui.py

@@ -333,6 +333,15 @@ class WebGLScene(base):
             else:
                 d["objects"].append(obj._GetWebguiData())
 
+        if 'center' in kwargs:
+            center = list(kwargs['center'])
+            if len(center) == 2:
+                center.append(0.)
+            d["mesh_center"] = center
+
+        if 'radius' in kwargs:
+            d["mesh_radius"] = kwargs['radius']
+
         return d
 
 
@@ -389,6 +398,7 @@ def Draw(obj, *args, show=True, **kwargs):
 
             html = scene.GenerateHTML()
             display(HTML(html))
+            return
         else:
             import webgui_jupyter_widgets as wjw
             from packaging.version import parse
@@ -400,8 +410,6 @@ def Draw(obj, *args, show=True, **kwargs):
                 scene.Draw(
                     kwargs_with_defaults["width"], kwargs_with_defaults["height"]
                 )
-            return scene
-    else:
     if "filename" in kwargs_with_defaults:
         scene.GenerateHTML(filename=kwargs_with_defaults["filename"])
     return scene

rules/CMakeLists.txt

@@ -8,6 +8,9 @@ if(EMSCRIPTEN)
   set(rules_command ${CMAKE_BINARY_DIR}/makerls)
 else(EMSCRIPTEN)
   add_executable(makerls rules/makerlsfile.cpp)
+  if(USE_CCACHE)
+    set_target_properties(makerls PROPERTIES RULE_LAUNCH_COMPILE "")
+  endif(USE_CCACHE)
   set(rules_command makerls)
 endif()
 

rules/tetrules.rls

@@ -139,7 +139,7 @@ endrule
 
 rule "Tetrahedron Vis a Vis Point (1)"
 
-quality 100
+quality 20
 
 mappoints
 (0, 0, 0);

setup.py

@@ -1,17 +1,24 @@
 import glob
 import os.path
+import os
 import sys
 import pathlib
 import sysconfig
+import importlib.metadata
 
 from skbuild import setup
 import skbuild.cmaker
 from subprocess import check_output
 
-setup_requires = ['pybind11-stubgen==2.5']
+setup_requires = ['pybind11-stubgen>=2.5', 'netgen-occt-devel']
 
 pyprefix = pathlib.Path(sys.prefix).as_posix()
 
+def find_occt_dir():
+    for f in importlib.metadata.files("netgen-occt-devel"):
+        if f.match("OpenCASCADEConfig.cmake"):
+            return f.locate().parent.resolve().absolute().as_posix()
+
 def install_filter(cmake_manifest):
     print(cmake_manifest)
     return cmake_manifest
@@ -28,14 +35,15 @@ def _patched_parse_manifests(self):
 # patch the parse_manifests function to point to the actual netgen cmake project within the superbuild
 skbuild.cmaker.CMaker._parse_manifests = _patched_parse_manifests
 
-git_version = check_output(['git', 'describe', '--tags']).decode('utf-8').strip()
+def is_dev_build():
-version = git_version[1:].split('-')
+    if 'NG_NO_DEV_PIP_VERSION' in os.environ:
-if len(version)>2:
+        return False
-    version = version[:2]
+    if 'CI_COMMIT_REF_NAME' in os.environ and os.environ['CI_COMMIT_REF_NAME'] == 'release':
-if len(version)>1:
+        return False
-    version = '.post'.join(version) + '.dev'
+    return True
-else:
+
-    version = version[0]
+git_version = check_output([sys.executable, os.path.join('tests', 'utils.py'), '--get-git-version']).decode('utf-8').strip()
+version = check_output([sys.executable, os.path.join('tests', 'utils.py'), '--get-version']).decode('utf-8').strip()
 
 py_install_dir = os.path.relpath(sysconfig.get_path('platlib'), sysconfig.get_path('data')).replace('\\','/')
 
@@ -44,6 +52,7 @@ arch = None
 cmake_args = [
         f'-DNETGEN_VERSION_GIT={git_version}',
         f'-DNETGEN_VERSION_PYTHON={version}',
+        f'-DOpenCascade_DIR={find_occt_dir()}',
     ]
 
 if 'NETGEN_ARCH' in os.environ and os.environ['NETGEN_ARCH'] == 'avx2':
@@ -129,7 +138,8 @@ cmake_args += [
         '-DUSE_GUI=ON',
         '-DUSE_NATIVE_ARCH=OFF',
         '-DBUILD_ZLIB=ON',
-        '-DBUILD_OCC=ON',
+        '-DZLIB_USE_STATIC_LIBS=ON',
+        '-DBUILD_OCC=OFF',
         '-DUSE_OCC=ON',
         '-DBUILD_FOR_CONDA=ON',
         f'-DNETGEN_PYTHON_PACKAGE_NAME={name}',
@@ -146,6 +156,7 @@ setup(
     license="LGPL2.1",
     packages=packages,
     #package_dir={'netgen': 'python'},
+    install_requires=[f"netgen-occt=={importlib.metadata.version('netgen-occt-devel')}"],
     tests_require=['pytest'],
     #include_package_data=True,
     cmake_process_manifest_hook=install_filter,

tests/build_pip.ps1

@@ -10,6 +10,12 @@ $env:NETGEN_ARCH = 'avx2'
 
 $pydir=$args[0]
 & $pydir\python.exe --version
+& $pydir\python.exe -m pip install packaging
+& $pydir\python.exe tests\utils.py --check-pip
+if ($LASTEXITCODE -ne 0) {
+    exit 0
+}
 & $pydir\python.exe -m pip install scikit-build wheel numpy twine pybind11-stubgen
+& $pydir\python.exe -m pip install --upgrade netgen-occt==7.8.1 netgen-occt-devel==7.8.1
 & $pydir\python setup.py bdist_wheel -G"Visual Studio 16 2019"
 & $pydir\python -m twine upload dist\*.whl

tests/build_pip.sh

@@ -9,7 +9,7 @@ dpkg-deb -R openmpi-dev.deb /opt/openmpi
 mv /opt/openmpi/usr/lib/x86_64-linux-gnu/openmpi/include /opt/openmpi/include
 
 
-curl http://ftp.de.debian.org/debian/pool/main/m/mpich/libmpich-dev_4.2.0-5.1_amd64.deb -o mpich.deb
+curl http://ftp.de.debian.org/debian/pool/main/m/mpich/libmpich-dev_4.2.1-2_amd64.deb -o mpich.deb
 dpkg-deb -R mpich.deb /opt/mpich
 mv /opt/mpich/usr/lib/x86_64-linux-gnu/mpich/include /opt/mpich/include
 
@@ -19,23 +19,27 @@ export NETGEN_CCACHE=1
 
 /opt/python/cp39-cp39/bin/python tests/fix_auditwheel_policy.py
 
-for pyversion in 38 39 310 311 312
+for pyversion in 312 311 310 39 38
 do
     export PYDIR="/opt/python/cp${pyversion}-cp${pyversion}/bin"
     echo $PYDIR
-    $PYDIR/pip install -U pytest-check numpy wheel scikit-build pybind11-stubgen
+    $PYDIR/pip install requests packaging
+    $PYDIR/python3 ./tests/utils.py --check-pip || continue
+    $PYDIR/pip install -U pytest-check numpy wheel scikit-build pybind11-stubgen netgen-occt==7.8.1 netgen-occt-devel==7.8.1
     $PYDIR/pip install -i https://pypi.anaconda.org/mpi4py/simple/ --pre mpi4py
 
     rm -rf _skbuild
     NETGEN_ARCH=avx2 $PYDIR/pip wheel .
-    auditwheel repair netgen_mesher*-cp${pyversion}-*.whl
+    mkdir -p wheelhouse
-    rm netgen_mesher-*.whl
+    #auditwheel repair netgen_mesher*-cp${pyversion}-*.whl
+    rename linux_x86_64 manylinux_2_17_x86_64.manylinux2014_x86_64 netgen_mesher*-cp${pyversion}-*.whl
+    mv netgen_mesher*-cp${pyversion}-*.whl wheelhouse/
 
     $PYDIR/pip install wheelhouse/netgen_mesher*-cp${pyversion}-*.whl
     $PYDIR/python3 -c 'import netgen'
+    $PYDIR/pip install -U twine
+    $PYDIR/twine upload --skip-existing wheelhouse/netgen_mesher*-cp${pyversion}*manylinux*.whl
     #cd ../tests/pytest
     #$PYDIR/python3 -m pytest
 done
 
-$PYDIR/pip install -U twine
-$PYDIR/twine upload wheelhouse/*manylinux*.whl

tests/build_pip_mac.sh

@@ -7,7 +7,10 @@ export PATH=$PYDIR:/Applications/CMake.app/Contents/bin:$PATH
 export NETGEN_CCACHE=1
 
 $PYDIR/python3 --version
-$PYDIR/pip3 install --user numpy twine scikit-build wheel pybind11-stubgen
+$PYDIR/python3 -m pip install packaging
+$PYDIR/python3 tests/utils.py --check-pip || exit 0
+$PYDIR/python3 -m pip install --user numpy twine scikit-build wheel pybind11-stubgen
+$PYDIR/python3 -m pip install --user -U netgen-occt==7.8.1 netgen-occt-devel==7.8.1
 
 export CMAKE_OSX_ARCHITECTURES='arm64;x86_64'
 export NETGEN_ARCH='avx2'

tests/pytest/test_array.py

@@ -1,3 +1,4 @@
+import netgen
 from pyngcore import *
 from numpy import sort, array
 

tests/pytest/test_bitarray.py

@@ -1,3 +1,4 @@
+import netgen
 from pyngcore import BitArray
 
 def test_bitarray():

tests/utils.py

@@ -0,0 +1,123 @@
+import argparse
+import os
+import requests
+import sys
+from subprocess import check_output
+from packaging import tags
+from packaging.utils import parse_wheel_filename
+
+
+_sys_tags = None
+
+
+def _is_wheel_compatible(wheel_filename: str):
+    global _sys_tags
+    try:
+        if _sys_tags is None:
+            _sys_tags = set(tags.sys_tags())
+
+        for tag in parse_wheel_filename(wheel_filename)[-1]:
+            if tag in _sys_tags:
+                return True
+
+        return False
+    except Exception as e:
+        print(f"Error parsing wheel file: {e}")
+        return False
+
+
+def is_package_available(package_name, version):
+    url = f"https://pypi.org/pypi/{package_name}/{version}/json"
+
+    try:
+        response = requests.get(url)
+        if response.status_code != 200:
+            return False
+
+        data = response.json()
+
+        for file_info in data["urls"]:
+            name = file_info.get("filename", "")
+            if _is_wheel_compatible(name):
+                return True
+
+        return False
+
+    except requests.RequestException as e:
+        print(f"Error checking package: {e}")
+        return False
+
+
+def is_dev_build():
+    if "NG_NO_DEV_PIP_VERSION" in os.environ:
+        return False
+    if (
+        "CI_COMMIT_REF_NAME" in os.environ
+        and os.environ["CI_COMMIT_REF_NAME"] == "release"
+    ):
+        return False
+    return True
+
+
+def get_git_version(cwd):
+    return check_output(["git", "describe", "--tags"], cwd=cwd).decode("utf-8").strip()
+
+
+def get_version(cwd):
+    git_version = get_git_version(cwd)
+
+    version = git_version[1:].split("-")
+    if len(version) > 2:
+        version = version[:2]
+    if len(version) > 1:
+        version = ".post".join(version)
+        if is_dev_build():
+            version += ".dev0"
+    else:
+        version = version[0]
+
+    return version
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Netgen pip building utilities")
+    parser.add_argument(
+        "--check-pip",
+        action="store_true",
+        help="Check if package is on pypi already, fails with exit code 1 if available",
+    )
+    parser.add_argument(
+        "--get-git-version",
+        action="store_true",
+        help="Generate the current package git version string",
+    )
+    parser.add_argument(
+        "--get-version",
+        action="store_true",
+        help="Generate the current package version using git",
+    )
+    parser.add_argument("--dir", type=str, default=".", help="CWD to run git commands")
+    parser.add_argument(
+        "--package",
+        type=str,
+        default="netgen-mesher",
+        help="Package name to check on pypi",
+    )
+
+    args = parser.parse_args()
+
+    if args.get_git_version:
+        print(get_git_version(args.dir))
+    elif args.get_version:
+        print(get_version(args.dir))
+    elif args.check_pip:
+        version = get_version(args.dir)
+        if is_package_available(args.package, version):
+            print(f"{args.package}=={version} is already on pypi")
+            sys.exit(1)
+    else:
+        print("no action")
+
+
+if __name__ == "__main__":
+    main()