Merge branch 'master' into refactor_geo_meshing

.gitlab-ci.yml

@@ -56,6 +56,9 @@ build_win:
         -G Ninja
         -DCMAKE_INSTALL_PREFIX=%INSTALL_DIR%
         -DUSE_OCC=ON
+        -DOCC_LIBRARY=C:/install_opencascade_7.4.0_static/win64/vc14/lib/TKernel.lib
+        -DOCC_INCLUDE_DIR=C:/install_opencascade_7.4.0_static/inc
+        -DOCC_LINK_FREETYPE=ON
         -DUSE_CCACHE=ON
         -DENABLE_UNIT_TESTS=ON
         -DCMAKE_BUILD_TYPE=Release
@@ -239,6 +242,10 @@ build_mac:
         -DENABLE_UNIT_TESTS=ON
         -DCMAKE_OSX_DEPLOYMENT_TARGET=10.12
         -DCMAKE_OSX_SYSROOT=/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk
+        -DUSE_OCC=ON
+        -DOCC_LIBRARY=/usr/local/opt/opencascade-7.4.0/lib/libTKernel.a
+        -DOCC_INCLUDE_DIR=/usr/local/opt/opencascade-7.4.0/include/opencascade
+        -DOCC_LINK_FREETYPE=ON
     - make -j5 install
 
 test_mac:

cmake/cmake_modules/FindOpenCasCade.cmake

@@ -80,6 +80,14 @@ set(OCC_LIBRARY_NAMES
     TKXSBase
 )
 
+if(OCC_LINK_FREETYPE)
+    set(OCC_LIBRARY_NAMES ${OCC_LIBRARY_NAMES} freetype)
+endif(OCC_LINK_FREETYPE)
+
+if(OCC_VERSION_STRING VERSION_GREATER_EQUAL "7.3.0")
+    set(OCC_LIBRARY_NAMES ${OCC_LIBRARY_NAMES} TKVCAF)
+endif()
+
 foreach( libname ${OCC_LIBRARY_NAMES} )
     find_library( ${libname} ${libname} ${OCC_LIBRARY_DIR} )
     set(OCC_LIBRARIES ${OCC_LIBRARIES} ${${libname}})

libsrc/core/ngcore_api.hpp

@@ -74,6 +74,22 @@ inline void operator delete[]( void* ptr, std::align_val_t al ) noexcept
     delete[] (char*)ptr;
 }
 
+inline void operator delete  ( void* ptr, std::size_t sz, std::align_val_t al ) noexcept
+{
+  if (int(al) > __STDCPP_DEFAULT_NEW_ALIGNMENT__)
+     _mm_free(ptr);
+  else
+    delete (char*)ptr;
+}
+
+inline void operator delete[]( void* ptr, std::size_t sz, std::align_val_t al ) noexcept
+{
+  if (int(al) > __STDCPP_DEFAULT_NEW_ALIGNMENT__)
+     _mm_free(ptr);
+  else
+    delete[] (char*)ptr;
+}
+
 #endif // __MAC_OS_X_VERSION_MIN_REQUIRED
 #endif // __MAC_OS_X_VERSION_MIN_REQUIRED < 101300
 

libsrc/core/taskmanager.cpp

@@ -35,11 +35,11 @@ namespace ngcore
   int TaskManager :: num_threads = 1;
 
   
-#ifndef __clang__      
+  // #ifndef __clang__      
   thread_local int TaskManager :: thread_id = 0;
-#else
-  __thread int TaskManager :: thread_id;
-#endif
+  // #else
+  // __thread int TaskManager :: thread_id;
+  // #endif
   
   const function<void(TaskInfo&)> * TaskManager::func;
   const function<void()> * TaskManager::startup_function = nullptr;

libsrc/core/taskmanager.hpp

@@ -72,11 +72,11 @@ namespace ngcore
 
 
 
-#ifndef __clang__      
+    // #ifndef __clang__      
     static thread_local int thread_id;
-#else
-    static __thread int thread_id;
-#endif
+    // #else
+    // static __thread int thread_id;
+    // #endif
     
     NGCORE_API static bool use_paje_trace;
   public:

libsrc/meshing/curvedelems.cpp

@@ -539,7 +539,7 @@ namespace netgen
       
 
   CurvedElements :: CurvedElements (const Mesh & amesh)
-    : mesh(amesh), geo(*mesh.GetGeometry())
+    : mesh(amesh)
   {
     order = 1;
     rational = 0;
@@ -555,6 +555,7 @@ namespace netgen
   void CurvedElements :: BuildCurvedElements(const Refinement * ref, int aorder,
                                              bool arational)
   {
+    auto & geo = *mesh.GetGeometry();
 
     ishighorder = 0;
     order = 1;

libsrc/meshing/curvedelems.hpp

@@ -17,7 +17,6 @@ class Refinement;
 class CurvedElements
 {
   const Mesh & mesh;
-  const NetgenGeometry& geo;
 
   NgArray<int> edgeorder;
   NgArray<int> faceorder;

libsrc/meshing/improve3.cpp

@@ -527,6 +527,246 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
 
 
 
+double MeshOptimize3d :: SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table<ElementIndex,PointIndex> & elementsonnode, Array<double> &elerrs, NgArray<INDEX_3> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only)
+{
+  double d_badness = 0.0;
+  int cnt = 0;
+
+  ArrayMem<ElementIndex, 20> hasbothpoints;
+
+  if (mesh.BoundaryEdge (pi1, pi2)) return 0.0;
+
+  for (ElementIndex ei : elementsonnode[pi1])
+    {
+      Element & el = mesh[ei];
+
+      if(el.IsDeleted()) return 0.0;
+      if (mesh[ei].GetType() != TET) return 0.0;
+
+      bool has1 = el.PNums().Contains(pi1);
+      bool has2 = el.PNums().Contains(pi2);
+
+      if (has1 && has2)
+          if (!hasbothpoints.Contains (ei))
+              hasbothpoints.Append (ei);
+    }
+
+  if(mp.only3D_domain_nr)
+      for(auto ei : hasbothpoints)
+          if(mp.only3D_domain_nr != mesh[ei].GetIndex())
+              return 0.0;
+
+
+  double bad1 = 0.0;
+  for (ElementIndex ei : hasbothpoints)
+      bad1 += CalcBad (mesh.Points(), mesh[ei], 0);
+
+  bool puretet = 1;
+  for (ElementIndex ei : hasbothpoints)
+      if (mesh[ei].GetType() != TET)
+          puretet = 0;
+  if (!puretet) return 0.0;
+
+  Point3d p1 = mesh[pi1];
+  Point3d p2 = mesh[pi2];
+
+  locfaces.SetSize(0);
+  for (ElementIndex ei : hasbothpoints)
+    {
+      const Element & el = mesh[ei];
+
+      for (int l = 0; l < 4; l++)
+          if (el[l] == pi1 || el[l] == pi2)
+            {
+              INDEX_3 i3;
+              Element2d face(TRIG);
+              el.GetFace (l+1, face);
+              for (int kk = 1; kk <= 3; kk++)
+                  i3.I(kk) = face.PNum(kk);
+              locfaces.Append (i3);
+            }
+    }
+
+  PointFunction1 pf (mesh.Points(), locfaces, mp, -1);
+  OptiParameters par;
+  par.maxit_linsearch = 50;
+  par.maxit_bfgs = 20;
+
+  Point3d pnew = Center (p1, p2);
+  Vector px(3);
+  px(0) = pnew.X();
+  px(1) = pnew.Y();
+  px(2) = pnew.Z();
+
+  if (bad1 > 0.1 * badmax)
+      BFGS (px, pf, par);
+
+  double bad2 = pf.Func (px);
+
+  pnew.X() = px(0);
+  pnew.Y() = px(1);
+  pnew.Z() = px(2);
+
+  mesh[ptmp] = Point<3>(pnew);
+
+  for (int k = 0; k < hasbothpoints.Size(); k++)
+    {
+      Element & oldel = mesh[hasbothpoints[k]];
+      Element newel1 = oldel;
+      Element newel2 = oldel;
+
+      oldel.flags.illegal_valid = 0;
+      newel1.flags.illegal_valid = 0;
+      newel2.flags.illegal_valid = 0;
+
+      for (int l = 0; l < 4; l++)
+        {
+          if (newel1[l] == pi2) newel1[l] = ptmp;
+          if (newel2[l] == pi1) newel2[l] = ptmp;
+        }
+
+      if (!mesh.LegalTet (oldel)) bad1 += 1e6;
+      if (!mesh.LegalTet (newel1)) bad2 += 1e6;
+      if (!mesh.LegalTet (newel2)) bad2 += 1e6;
+    }
+
+  d_badness = bad2-bad1;
+  if(check_only)
+      return d_badness;
+
+  if (d_badness<0.0)
+    {
+      cnt++;
+
+      PointIndex pinew = mesh.AddPoint (pnew);
+
+      for (ElementIndex ei : hasbothpoints)
+        {
+          Element & oldel = mesh[ei];
+          Element newel1 = oldel;
+          Element newel2 = oldel;
+
+          oldel.flags.illegal_valid = 0;
+          oldel.Delete();
+
+          newel1.flags.illegal_valid = 0;
+          newel2.flags.illegal_valid = 0;
+
+          for (int l = 0; l < 4; l++)
+            {
+              if (newel1[l] == pi2) newel1[l] = pinew;
+              if (newel2[l] == pi1) newel2[l] = pinew;
+            }
+
+          mesh.AddVolumeElement (newel1);
+          mesh.AddVolumeElement (newel2);
+        }
+    }
+  return d_badness;
+}
+
+void MeshOptimize3d :: SplitImprove (Mesh & mesh,
+				     OPTIMIZEGOAL goal)
+{
+  static Timer t("MeshOptimize3d::SplitImprove"); RegionTimer reg(t);
+  static Timer topt("Optimize");
+  static Timer tsearch("Search");
+
+  // return SplitImproveSequential(mesh, goal);
+
+  int np = mesh.GetNP();
+  int ne = mesh.GetNE();
+  double bad = 0.0;
+  double badmax = 0.0;
+
+  auto elementsonnode = mesh.CreatePoint2ElementTable();
+
+  Array<double> elerrs(ne);
+
+  const char * savetask = multithread.task;
+  multithread.task = "Split Improve";
+
+  PrintMessage (3, "SplitImprove");
+  (*testout)  << "start SplitImprove" << "\n";
+
+  ParallelFor( mesh.VolumeElements().Range(), [&] (ElementIndex ei) NETGEN_LAMBDA_INLINE
+    {
+      if(mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex())
+          return;
+
+      elerrs[ei] = CalcBad (mesh.Points(), mesh[ei], 0);
+      bad += elerrs[ei];
+      AtomicMax(badmax, elerrs[ei]);
+    });
+
+  if (goal == OPT_QUALITY)
+    {
+      bad = mesh.CalcTotalBad (mp);
+      (*testout) << "Total badness = " << bad << endl;
+    }
+
+  Array<std::tuple<PointIndex,PointIndex>> edges;
+  BuildEdgeList(mesh, elementsonnode, edges);
+
+  // Find edges with improvement
+  Array<std::tuple<double, int>> candidate_edges(edges.Size());
+  std::atomic<int> improvement_counter(0);
+  auto ptmp = mesh.AddPoint( {0,0,0} );
+
+  tsearch.Start();
+  ParallelForRange(Range(edges), [&] (auto myrange)
+  {
+    NgArray<INDEX_3> locfaces;
+
+    for(auto i : myrange)
+    {
+      auto [p0,p1] = edges[i];
+      double d_badness = SplitImproveEdge (mesh, goal, elementsonnode, elerrs, locfaces, badmax, p0, p1, ptmp, true);
+      if(d_badness<0.0)
+      {
+        int index = improvement_counter++;
+        candidate_edges[index] = make_tuple(d_badness, i);
+      }
+    }
+  }, ngcore::TasksPerThread(4));
+  tsearch.Stop();
+
+  auto edges_with_improvement = candidate_edges.Part(0, improvement_counter.load());
+
+  QuickSort(edges_with_improvement);
+  PrintMessage(5, edges.Size(), " edges");
+  PrintMessage(5, edges_with_improvement.Size(), " edges with improvement");
+
+  // Apply actual optimizations
+  topt.Start();
+  int cnt = 0;
+  NgArray<INDEX_3> locfaces;
+  for(auto [d_badness, ei] : edges_with_improvement)
+  {
+      auto [p0,p1] = edges[ei];
+      if (SplitImproveEdge (mesh, goal, elementsonnode, elerrs, locfaces, badmax, p0, p1, ptmp, false) < 0.0)
+        cnt++;
+  }
+  topt.Stop();
+  mesh.Compress();
+  PrintMessage (5, cnt, " splits performed");
+  (*testout) << "Splitt - Improve done" << "\n";
+
+  if (goal == OPT_QUALITY)
+    {
+      bad = mesh.CalcTotalBad (mp);
+      (*testout) << "Total badness = " << bad << endl;
+
+      int cntill = 0;
+      ne = mesh.GetNE();
+      for (ElementIndex ei = 0; ei < ne; ei++)
+        if (!mesh.LegalTet (mesh[ei]))
+          cntill++;
+      //      cout << cntill << " illegal tets" << endl;
+    }
+
+  multithread.task = savetask;
+}
 
 
 /*
@@ -534,7 +774,7 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
   If mesh quality is improved by inserting a node into an inner edge,
   the edge is split into two parts.
 */
-void MeshOptimize3d :: SplitImprove (Mesh & mesh,
+void MeshOptimize3d :: SplitImproveSequential (Mesh & mesh,
 				     OPTIMIZEGOAL goal)
 {
   static Timer t("MeshOptimize3d::SplitImprove"); RegionTimer reg(t);

libsrc/meshing/improve3.hpp

@@ -24,6 +24,9 @@ public:
   void CombineImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
 
   void SplitImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
+  void SplitImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
+  double SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table<ElementIndex,PointIndex> & elementsonnode, Array<double> &elerrs, NgArray<INDEX_3> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only=false);
+  
 
   double SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, const NgBitArray * working_elements, Table<ElementIndex,PointIndex> & elementsonnode, INDEX_3_HASHTABLE<int> & faces, PointIndex pi1, PointIndex pi2, bool check_only=false);
   void SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,

libsrc/occ/python_occ.cpp

@@ -7,6 +7,7 @@
 
 #include <meshing.hpp>
 #include <occgeom.hpp>
+#include <Standard_Version.hxx>
 
 using namespace netgen;
 
@@ -57,6 +58,7 @@ void CreateOCCParametersFromKwargs(OCCParameters& occparam, py::dict kwargs)
 
 DLL_HEADER void ExportNgOCC(py::module &m) 
 {
+  m.attr("occ_version") = OCC_VERSION_COMPLETE;
   py::class_<OCCGeometry, shared_ptr<OCCGeometry>, NetgenGeometry> (m, "OCCGeometry", R"raw_string(Use LoadOCCGeometry to load the geometry from a *.step file.)raw_string")
     .def(py::init<>())
     .def(py::init([] (const string& filename)

tests/dockerfile

@@ -1,5 +1,5 @@
-FROM ubuntu:18.04
+FROM ubuntu:19.10
 ENV DEBIAN_FRONTEND=noninteractive
 MAINTAINER Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>
-RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-pytest python3-numpy python3-tk clang-tidy python3-distutils clang liboce-ocaf-dev
+RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-pytest python3-numpy python3-tk clang-tidy python3-distutils clang libocct-data-exchange-dev
 ADD . /root/src/netgen

tests/pytest/test_tutorials.py

@@ -8,7 +8,7 @@ from pyngcore import TaskManager
 import json
 try:
     import netgen.occ as occ
-    has_occ = True
+    has_occ = occ.occ_version >= "7.4.0"
 except ImportError:
     has_occ = False