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

libsrc/core/archive.hpp

@@ -42,6 +42,32 @@ namespace ngcore
     operator T&() { return val; }
   };
 
+  // Helper to detect shared_from_this
+  template <typename T>
+  class has_shared_from_this2
+    {
+    private:
+      // typedef T* T_ptr;
+      template <typename C> static std::true_type test(decltype(((C*)nullptr)->shared_from_this()));
+      template <typename C> static std::false_type test(...);
+      
+    public:
+      // If the test returns true_type, then T has shared_from_this
+      static constexpr bool value = decltype(test<T>(0))::value;
+  };
+  
+  
+
+  
+  template <typename T, typename = void>
+  class has_shallow_archive : public std::false_type {};
+  
+  template <typename T>
+  class has_shallow_archive<T, std::void_t<decltype(T::shallow_archive)>>
+    : public std::is_same<decltype(T::shallow_archive), std::true_type> {};
+  
+
+  
 #ifdef NETGEN_PYTHON
   pybind11::object CastAnyToPy(const std::any& a);
 #endif // NETGEN_PYTHON
@@ -486,6 +512,13 @@ namespace ngcore
     template <typename T>
     Archive& operator & (std::shared_ptr<T>& ptr)
     {
+      if constexpr(has_shallow_archive<T>::value)
+        if (shallow_to_python)
+          {
+            Shallow (ptr);
+            return *this;
+          }
+          
       if(Output())
         {
           // save -2 for nullptr

libsrc/core/exception.hpp

@@ -92,10 +92,10 @@ namespace ngcore
       throw ngcore::RangeException(__FILE__ ":" NETGEN_CORE_NGEXEPTION_STR(__LINE__) "\t", int(value), int(min), int(max_plus_one)); }
 #define NETGEN_CHECK_SHAPE(a,b) \
   { if(a.Shape() != b.Shape()) \
-      throw ngcore::Exception(__FILE__": shape don't match"); }
+      throw ngcore::Exception(__FILE__ ":" NETGEN_CORE_NGEXEPTION_STR(__LINE__) "\t: shape don't match"); }
 #define NETGEN_CHECK_SAME(a,b) \
   { if(a != b) \
-      throw ngcore::Exception(__FILE__": not the same, a="+ToString(a) + ", b="+ToString(b)); }
+      throw ngcore::Exception(__FILE__ ":" NETGEN_CORE_NGEXEPTION_STR(__LINE__) "\t: not the same, a="+ToString(a) + ", b="+ToString(b)); }
 #define NETGEN_NOEXCEPT 
 #else // defined(NETGEN_ENABLE_CHECK_RANGE) && !defined(__CUDA_ARCH__)
 #define NETGEN_CHECK_RANGE(value, min, max)

libsrc/core/register_archive.hpp

@@ -73,8 +73,8 @@ namespace ngcore {
       };
 #ifdef NETGEN_PYTHON
     info.anyToPyCaster = [](const std::any &a) {
-      if constexpr(has_shared_from_this<T>::value) {
+      if constexpr(has_shared_from_this2<T>::value) {
-        std::shared_ptr<T> val = std::any_cast<std::shared_ptr<T>>(&a);
+        std::shared_ptr<T> val = std::any_cast<std::shared_ptr<T>>(a);
         return pybind11::cast(val);
       } else {
         const T* val = std::any_cast<T>(&a);

libsrc/core/signal.hpp

@@ -2,6 +2,7 @@
 #define NGCORE_SIGNALS_HPP
 
 #include <list>
+#include <map>
 #include <functional>
 
 namespace ngcore
@@ -43,6 +44,39 @@ namespace ngcore
     }
     inline bool GetEmitting() const { return is_emitting; }
   };
+
+
+
+
+  class SimpleSignal
+  {
+  private:
+    // std::map<void*,std::function<void()>> funcs;
+    std::list<std::pair<void*,std::function<void()>>> funcs;
+  public:
+    SimpleSignal() = default;
+
+    template<typename FUNC>
+    void Connect(void* var, FUNC f)
+    {
+      // funcs[var] = f;
+      funcs.push_back ( { var, f } );
+    }
+
+    void Remove(void* var)
+    {
+      // funcs.erase(var);
+      funcs.remove_if([&] (auto var_f) { return var_f.first==var; });
+    }
+
+    inline void Emit()
+    {
+      for (auto [key,f] : funcs)
+        f();
+    }
+  };
+
+  
 } // namespace ngcore
 
 #endif // NGCORE_SIGNALS_HPP

libsrc/interface/nginterface_v2.cpp

@@ -1234,7 +1234,7 @@ namespace netgen
   {
     NgLock meshlock (mesh->MajorMutex(), true);
     Refinement & ref = const_cast<Refinement&> (mesh->GetGeometry()->GetRefinement());
-    ::netgen::HPRefinement (*mesh, &ref, SPLIT_ALFELD, 1, 0.5, true, true);
+    ::netgen::HPRefinement (*mesh, &ref, SPLIT_ALFELD, 1, 1.0/3.0, true, true);
   }
 
   

libsrc/meshing/hpref_trig.hpp

@@ -775,18 +775,7 @@ HPRef_Struct reftrig_3singedges =
   reftrig_3singedges_newels
 };
 
-
+// HP_TRIG_ALFELD
-
-
-
-
-
-
-
-
-
-
-// HP_TRIG_3SINGEDGES
 int reftrig_Alfeld_splitedges[][3] =
 {
   { 0, 0, 0 }
@@ -818,3 +807,42 @@ HPRef_Struct reftrig_Alfeld =
   reftrig_Alfeld_newels
 };
 
+
+// HP_TRIG_POWELL
+int reftrig_Powell_splitedges[][3] =
+{
+  { 1, 2, 4 },
+  { 2, 3, 5 },
+  { 3, 1, 6 },
+  { 0, 0, 0 },
+};
+int reftrig_Powell_splitfaces[][4] =
+{
+  { 1, 2, 3, 7 },
+  { 0, 0, 0, 0 }
+};
+
+HPREF_ELEMENT_TYPE reftrig_Powell_newelstypes[] =
+{
+  HP_TRIG, HP_TRIG, HP_TRIG, HP_TRIG, HP_TRIG, HP_TRIG,
+  HP_NONE, 
+};
+int reftrig_Powell_newels[][8] =
+{
+  { 1, 4, 7 }, 
+  { 4, 2, 7 }, 
+  { 2, 5, 7 },
+  { 5, 3, 7 },
+  { 3, 6, 7 },
+  { 6, 1, 7 },
+};
+HPRef_Struct reftrig_Powell =
+{
+  HP_TRIG,
+  reftrig_Powell_splitedges, 
+  reftrig_Powell_splitfaces,
+  0,
+  reftrig_Powell_newelstypes,
+  reftrig_Powell_newels
+};
+

libsrc/meshing/hprefinement.cpp

@@ -185,6 +185,8 @@ namespace netgen
  
       case HP_TRIG_ALFELD:
         hps = &reftrig_Alfeld; break;
+      case HP_TRIG_POWELL:
+        hps = &reftrig_Powell; break;
 
         
       case HP_QUAD:
@@ -674,9 +676,7 @@ namespace netgen
     INDEX_2_HASHTABLE<int> newpts(elements.Size()+1);
     INDEX_3_HASHTABLE<int> newfacepts(elements.Size()+1);
 
-    // prepare new points  
+    double fac2 = max(0.001,min(1.0/3,fac1)); // factor for face points
-    
-    fac1 = max(0.001,min(0.33,fac1));
     PrintMessage(3, " in HP-REFINEMENT with fac1 ", fac1); 
     *testout << " in HP-REFINEMENT with fac1 " << fac1 <<  endl; 
    
@@ -726,8 +726,8 @@ namespace netgen
 		{
 		  Point<3> np; 
 		  	for( int l=0;l<3;l++)
-			  np(l) = (1-2*fac1)*mesh.Point(i3.I1())(l) 
+			  np(l) = (1-2*fac2)*mesh.Point(i3.I1())(l) 
-			    + fac1*mesh.Point(i3.I2())(l)  + fac1*mesh.Point(i3.I3())(l);  
+			    + fac2*mesh.Point(i3.I2())(l)  + fac2*mesh.Point(i3.I3())(l);  
 		  int npi = mesh.AddPoint (np);
 		  newfacepts.Set (i3, npi);
 		}
@@ -815,9 +815,9 @@ namespace netgen
 
 	      for (int l = 0; l < 3; l++)
 		newparam[hprs->splitfaces[j][3]-1][l] =
-		  (1-2*fac1) * el.param[hprs->splitfaces[j][0]-1][l] + 
+		  (1-2*fac2) * el.param[hprs->splitfaces[j][0]-1][l] + 
-		  fac1 * el.param[hprs->splitfaces[j][1]-1][l] + 
+		  fac2 * el.param[hprs->splitfaces[j][1]-1][l] + 
-		  fac1 * el.param[hprs->splitfaces[j][2]-1][l];
+		  fac2 * el.param[hprs->splitfaces[j][2]-1][l];
 	      j++;
 	    }
 	// split elements
@@ -1906,6 +1906,8 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
 
     if (act_ref == 1 && split == SPLIT_ALFELD)
       sing = true;   
+    if (act_ref == 1 && split == SPLIT_POWELL)
+      sing = true;   
     if(sing==0) return(sing); 
 
     int cnt_undef = 0, cnt_nonimplement = 0;
@@ -1969,8 +1971,10 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
               if (split == SPLIT_HP)
                 hpel.type = ClassifyTrig(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, 
                                          faces, face_edges, surf_edges, facepoint, dim, fd);    
-              else
+              else if (split == SPLIT_ALFELD)
                 hpel.type = HP_TRIG_ALFELD;
+              else if (split == SPLIT_POWELL)
+                hpel.type = HP_TRIG_POWELL;
 
 	      dd = 2;
 	      break; 

libsrc/meshing/hprefinement.hpp

@@ -46,6 +46,7 @@ enum HPREF_ELEMENT_TYPE {
   HP_TRIG_3SINGEDGES = 40,
 
   HP_TRIG_ALFELD,
+  HP_TRIG_POWELL,
   
   HP_QUAD = 50,
   HP_QUAD_SINGCORNER,
@@ -347,7 +348,7 @@ public:
 };
 
 
-enum SplittingType { SPLIT_HP, SPLIT_ALFELD };
+enum SplittingType { SPLIT_HP, SPLIT_ALFELD, SPLIT_POWELL};
 
 DLL_HEADER extern void HPRefinement (Mesh & mesh, Refinement * ref, SplittingType split, int levels,
 			  double fac1=0.125, bool setorders=true, bool ref_level = false);

libsrc/meshing/python_mesh.cpp

@@ -1350,7 +1350,23 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
            }), py::arg("adaptive")=false, py::call_guard<py::gil_scoped_release>())
     
     .def("ZRefine", &Mesh::ZRefine)
-
+    .def("Split2Tets", &Mesh::Split2Tets)
+    .def ("SplitAlfeld", FunctionPointer
+          ([](Mesh & self)
+           {
+            NgLock meshlock (self.MajorMutex(), true);
+            Refinement & ref = const_cast<Refinement&> (self.GetGeometry()->GetRefinement());
+            ::netgen::HPRefinement (self, &ref, SPLIT_ALFELD, 1, 0.5, true, true);
+           }
+           ), py::call_guard<py::gil_scoped_release>())
+    .def ("SplitPowellSabin", FunctionPointer
+          ([](Mesh & self)
+           {
+            NgLock meshlock (self.MajorMutex(), true);
+            Refinement & ref = const_cast<Refinement&> (self.GetGeometry()->GetRefinement());
+            ::netgen::HPRefinement (self, &ref, SPLIT_POWELL, 1, 0.5, true, true);
+           }
+           ), py::call_guard<py::gil_scoped_release>())
     .def ("SecondOrder", [](Mesh & self)
           {
             self.GetGeometry()->GetRefinement().MakeSecondOrder(self);

libsrc/occ/occ_edge.cpp

@@ -61,7 +61,12 @@ namespace netgen
     void OCCEdge::ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const
     {
         auto pnt = ng2occ(p);
-        GeomAPI_ProjectPointOnCurve proj(pnt, curve, s0, s1);
+        // extend the projection parameter range, else projection might fail
+        // for an endpoint
+        // see discussion here: https://forum.ngsolve.org/t/how-to-apply-occidentification-correctly/2555
+        // I do not see a better way using occ tolerances?
+        double eps = 1e-7 * (s1-s0);
+        GeomAPI_ProjectPointOnCurve proj(pnt, curve, s0-eps, s1+eps);
         pnt = proj.NearestPoint();
         if(gi)
             gi->dist = (proj.LowerDistanceParameter() - s0)/(s1-s0);