From 48eb4fed07c9bf72d14dae93963d249785da413f Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Wed, 10 Jan 2024 16:31:05 +0100 Subject: [PATCH 01/14] add tolerance to occ-edge projection --- libsrc/occ/occ_edge.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/libsrc/occ/occ_edge.cpp b/libsrc/occ/occ_edge.cpp index 0c907d78..4805bb17 100644 --- a/libsrc/occ/occ_edge.cpp +++ b/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); From 2d2503bbbb48d7b5a20d66461ffe656265ba40bb Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Sat, 13 Jan 2024 21:15:55 +0100 Subject: [PATCH 02/14] auto-shallow shared_ptr with enable_shared_from_this --- libsrc/core/archive.hpp | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/libsrc/core/archive.hpp b/libsrc/core/archive.hpp index e9f28ef9..52e7ed83 100644 --- a/libsrc/core/archive.hpp +++ b/libsrc/core/archive.hpp @@ -42,6 +42,21 @@ namespace ngcore operator T&() { return val; } }; + // Helper to detect shared_from_this + template + class has_shared_from_this2 + { + private: + // typedef T* T_ptr; + template static std::true_type test(decltype(((C*)nullptr)->shared_from_this())); + template static std::false_type test(...); + + public: + // If the test returns true_type, then T has shared_from_this + static constexpr bool value = decltype(test(0))::value; + }; + + #ifdef NETGEN_PYTHON pybind11::object CastAnyToPy(const std::any& a); #endif // NETGEN_PYTHON @@ -486,6 +501,12 @@ namespace ngcore template Archive& operator & (std::shared_ptr& ptr) { + if (has_shared_from_this2::value && shallow_to_python) + { + Shallow (ptr); + return *this; + } + if(Output()) { // save -2 for nullptr From 1ff8c97b1dafb0750f71baa22c47029d557a6c6e Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Sun, 14 Jan 2024 04:33:55 +0100 Subject: [PATCH 03/14] fix has_shared_from_this any_cast --- libsrc/core/register_archive.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/libsrc/core/register_archive.hpp b/libsrc/core/register_archive.hpp index 7fbc4c87..b7be05d9 100644 --- a/libsrc/core/register_archive.hpp +++ b/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::value) { - std::shared_ptr val = std::any_cast>(&a); + if constexpr(has_shared_from_this2::value) { + std::shared_ptr val = std::any_cast>(a); return pybind11::cast(val); } else { const T* val = std::any_cast(&a); From 6b346926ec638edf06507abd01fb5c1ffeec8ed0 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Sun, 14 Jan 2024 04:52:19 +0100 Subject: [PATCH 04/14] if constexpr --- libsrc/core/archive.hpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/libsrc/core/archive.hpp b/libsrc/core/archive.hpp index 52e7ed83..6db37a07 100644 --- a/libsrc/core/archive.hpp +++ b/libsrc/core/archive.hpp @@ -501,11 +501,12 @@ namespace ngcore template Archive& operator & (std::shared_ptr& ptr) { - if (has_shared_from_this2::value && shallow_to_python) - { - Shallow (ptr); - return *this; - } + if constexpr(has_shared_from_this2::value) + if (shallow_to_python) + { + Shallow (ptr); + return *this; + } if(Output()) { From c0d394ebf55de436fd589e5868e490e9369fd4d1 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Mon, 15 Jan 2024 08:16:14 +0100 Subject: [PATCH 05/14] introduce 'shallow_archive' member --- libsrc/core/archive.hpp | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/libsrc/core/archive.hpp b/libsrc/core/archive.hpp index 6db37a07..d7546060 100644 --- a/libsrc/core/archive.hpp +++ b/libsrc/core/archive.hpp @@ -57,6 +57,17 @@ namespace ngcore }; + + + template + class has_shallow_archive : public std::false_type {}; + + template + class has_shallow_archive> + : public std::is_same {}; + + + #ifdef NETGEN_PYTHON pybind11::object CastAnyToPy(const std::any& a); #endif // NETGEN_PYTHON @@ -501,7 +512,7 @@ namespace ngcore template Archive& operator & (std::shared_ptr& ptr) { - if constexpr(has_shared_from_this2::value) + if constexpr(has_shallow_archive::value) if (shallow_to_python) { Shallow (ptr); From 890f59b8b4c483d4e3cc9988ac69bd43801604ad Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Mon, 15 Jan 2024 11:42:00 +0000 Subject: [PATCH 06/14] Exposed Alfeld splits --- libsrc/meshing/python_mesh.cpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 1a50a250..284b7997 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -1350,7 +1350,15 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) }), py::arg("adaptive")=false, py::call_guard()) .def("ZRefine", &Mesh::ZRefine) - + .def("Split2Tets", &Mesh::Split2Tets) + .def ("SplitAlfeld", FunctionPointer + ([](Mesh & self) + { + NgLock meshlock (self.MajorMutex(), true); + Refinement & ref = const_cast (self.GetGeometry()->GetRefinement()); + ::netgen::HPRefinement (self, &ref, SPLIT_ALFELD, 1, 0.5, true, true); + } + ), py::call_guard()) .def ("SecondOrder", [](Mesh & self) { self.GetGeometry()->GetRefinement().MakeSecondOrder(self); From 696620828fd4b36fba0edb946ed7d12e700c0941 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Mon, 15 Jan 2024 21:42:18 +0100 Subject: [PATCH 07/14] don't restrict refinement parameter in HPRefinement (more user responsibility) --- libsrc/meshing/hprefinement.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp index 1fae823f..066d1942 100644 --- a/libsrc/meshing/hprefinement.cpp +++ b/libsrc/meshing/hprefinement.cpp @@ -676,7 +676,7 @@ namespace netgen // prepare new points - fac1 = max(0.001,min(0.33,fac1)); + // fac1 = max(0.001,min(0.33,fac1)); PrintMessage(3, " in HP-REFINEMENT with fac1 ", fac1); *testout << " in HP-REFINEMENT with fac1 " << fac1 << endl; From 00747fb947601cf1e8749b4808c34551dc45e552 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Mon, 15 Jan 2024 22:43:18 +0000 Subject: [PATCH 08/14] Powell Sabin splits --- libsrc/meshing/hpref_trig.hpp | 52 +++++++++++++++++++++++++-------- libsrc/meshing/hprefinement.cpp | 8 ++++- libsrc/meshing/hprefinement.hpp | 3 +- libsrc/meshing/python_mesh.cpp | 8 +++++ 4 files changed, 57 insertions(+), 14 deletions(-) diff --git a/libsrc/meshing/hpref_trig.hpp b/libsrc/meshing/hpref_trig.hpp index 52f3582e..8c527c48 100644 --- a/libsrc/meshing/hpref_trig.hpp +++ b/libsrc/meshing/hpref_trig.hpp @@ -775,18 +775,7 @@ HPRef_Struct reftrig_3singedges = reftrig_3singedges_newels }; - - - - - - - - - - - -// HP_TRIG_3SINGEDGES +// HP_TRIG_ALFELD 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 +}; + diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp index 1fae823f..bb250aee 100644 --- a/libsrc/meshing/hprefinement.cpp +++ b/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: @@ -1906,6 +1908,8 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & 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 +1973,10 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & 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; diff --git a/libsrc/meshing/hprefinement.hpp b/libsrc/meshing/hprefinement.hpp index 1506cfb4..3cc7ca9e 100644 --- a/libsrc/meshing/hprefinement.hpp +++ b/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); diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 284b7997..886122e6 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -1359,6 +1359,14 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) ::netgen::HPRefinement (self, &ref, SPLIT_ALFELD, 1, 0.5, true, true); } ), py::call_guard()) + .def ("SplitPowellSabin", FunctionPointer + ([](Mesh & self) + { + NgLock meshlock (self.MajorMutex(), true); + Refinement & ref = const_cast (self.GetGeometry()->GetRefinement()); + ::netgen::HPRefinement (self, &ref, SPLIT_POWELL, 1, 0.5, true, true); + } + ), py::call_guard()) .def ("SecondOrder", [](Mesh & self) { self.GetGeometry()->GetRefinement().MakeSecondOrder(self); From 87c4e543ad008bbb7f394de0d0b6ca6b8c43d649 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Mon, 15 Jan 2024 23:54:44 +0000 Subject: [PATCH 09/14] Introduced fac2 to fix issue with face splits --- libsrc/meshing/hprefinement.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp index 2c59be5d..dc7988c5 100644 --- a/libsrc/meshing/hprefinement.cpp +++ b/libsrc/meshing/hprefinement.cpp @@ -677,8 +677,9 @@ namespace netgen INDEX_3_HASHTABLE newfacepts(elements.Size()+1); // prepare new points - // fac1 = max(0.001,min(0.33,fac1)); + double fac2; + fac2 = max(0.001,min(0.33,fac1)); PrintMessage(3, " in HP-REFINEMENT with fac1 ", fac1); *testout << " in HP-REFINEMENT with fac1 " << fac1 << endl; @@ -728,8 +729,8 @@ namespace netgen { Point<3> np; for( int l=0;l<3;l++) - np(l) = (1-2*fac1)*mesh.Point(i3.I1())(l) - + fac1*mesh.Point(i3.I2())(l) + fac1*mesh.Point(i3.I3())(l); + np(l) = (1-2*fac2)*mesh.Point(i3.I1())(l) + + fac2*mesh.Point(i3.I2())(l) + fac2*mesh.Point(i3.I3())(l); int npi = mesh.AddPoint (np); newfacepts.Set (i3, npi); } @@ -817,9 +818,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] + - fac1 * el.param[hprs->splitfaces[j][1]-1][l] + - fac1 * el.param[hprs->splitfaces[j][2]-1][l]; + (1-2*fac2) * el.param[hprs->splitfaces[j][0]-1][l] + + fac2 * el.param[hprs->splitfaces[j][1]-1][l] + + fac2 * el.param[hprs->splitfaces[j][2]-1][l]; j++; } // split elements From 29f0a5d647c9674cbda8a083785a498acfb4cba4 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 16 Jan 2024 10:17:19 +0100 Subject: [PATCH 10/14] simple signal without smart pointers --- libsrc/core/signal.hpp | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/libsrc/core/signal.hpp b/libsrc/core/signal.hpp index 082dd32b..bd9cd6f0 100644 --- a/libsrc/core/signal.hpp +++ b/libsrc/core/signal.hpp @@ -2,6 +2,7 @@ #define NGCORE_SIGNALS_HPP #include +#include #include namespace ngcore @@ -43,6 +44,36 @@ namespace ngcore } inline bool GetEmitting() const { return is_emitting; } }; + + + + + class SimpleSignal + { + private: + std::map> funcs; + public: + SimpleSignal() = default; + + template + void Connect(void* var, FUNC f) + { + funcs[var] = f; + } + + void Remove(void* var) + { + funcs.erase(var); + } + + inline void Emit() + { + for (auto [key,f] : funcs) + f(); + } + }; + + } // namespace ngcore #endif // NGCORE_SIGNALS_HPP From 6c3fcf018896118c12c364273ab5727b1e7642c3 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 16 Jan 2024 12:42:42 +0100 Subject: [PATCH 11/14] Alfeld split uses sub-division factor 1/3 --- libsrc/interface/nginterface_v2.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp index 99485ffc..970070b6 100644 --- a/libsrc/interface/nginterface_v2.cpp +++ b/libsrc/interface/nginterface_v2.cpp @@ -1234,7 +1234,7 @@ namespace netgen { NgLock meshlock (mesh->MajorMutex(), true); Refinement & ref = const_cast (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); } From cb7759cd0b00c3bce8bc6f06e01222c9dc935c7f Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 16 Jan 2024 12:43:23 +0100 Subject: [PATCH 12/14] line number in NETGEN_CHECK_SAME macro --- libsrc/core/exception.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/libsrc/core/exception.hpp b/libsrc/core/exception.hpp index 94885963..cf068aaf 100644 --- a/libsrc/core/exception.hpp +++ b/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) From 5a4b89c1ed19b14f69081851219fab5a69496a67 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 16 Jan 2024 12:54:06 +0100 Subject: [PATCH 13/14] merge hp face-refinement with limiting fac2 --- libsrc/meshing/hprefinement.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp index dc7988c5..0446f25b 100644 --- a/libsrc/meshing/hprefinement.cpp +++ b/libsrc/meshing/hprefinement.cpp @@ -676,10 +676,7 @@ namespace netgen INDEX_2_HASHTABLE newpts(elements.Size()+1); INDEX_3_HASHTABLE newfacepts(elements.Size()+1); - // prepare new points - // fac1 = max(0.001,min(0.33,fac1)); - double fac2; - fac2 = max(0.001,min(0.33,fac1)); + double fac2 = max(0.001,min(1.0/3,fac1)); // factor for face points PrintMessage(3, " in HP-REFINEMENT with fac1 ", fac1); *testout << " in HP-REFINEMENT with fac1 " << fac1 << endl; From eb90c6ed3be207526deaed6b5bc5846d117bb482 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Thu, 18 Jan 2024 19:43:08 +0100 Subject: [PATCH 14/14] use list instead of map to keep order --- libsrc/core/signal.hpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/libsrc/core/signal.hpp b/libsrc/core/signal.hpp index bd9cd6f0..a994f48c 100644 --- a/libsrc/core/signal.hpp +++ b/libsrc/core/signal.hpp @@ -51,19 +51,22 @@ namespace ngcore class SimpleSignal { private: - std::map> funcs; + // std::map> funcs; + std::list>> funcs; public: SimpleSignal() = default; template void Connect(void* var, FUNC f) { - funcs[var] = f; + // funcs[var] = f; + funcs.push_back ( { var, f } ); } void Remove(void* var) { - funcs.erase(var); + // funcs.erase(var); + funcs.remove_if([&] (auto var_f) { return var_f.first==var; }); } inline void Emit()