From b7b168e265b4a381abbd008b597b8f32eaa7c4ce Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 31 Dec 2024 12:16:53 +0100 Subject: [PATCH 01/19] base Index template --- libsrc/core/array.hpp | 4 + libsrc/interface/writegmsh2.cpp | 32 ++++---- libsrc/meshing/bisect.cpp | 2 - libsrc/meshing/clusters.cpp | 8 +- libsrc/meshing/improve3.cpp | 4 +- libsrc/meshing/meshclass.cpp | 6 +- libsrc/meshing/meshing3.cpp | 2 - libsrc/meshing/meshtype.hpp | 135 +++++++++++++++++++++++++++++++- libsrc/meshing/parallelmesh.cpp | 13 ++- libsrc/meshing/refine.cpp | 15 ++-- libsrc/meshing/topology.hpp | 7 ++ libsrc/occ/occmeshsurf.cpp | 19 +++-- libsrc/occ/occmeshsurf.hpp | 17 +++- libsrc/occ/occpkg.cpp | 18 +++-- nglib/nglib.cpp | 6 +- 15 files changed, 227 insertions(+), 61 deletions(-) diff --git a/libsrc/core/array.hpp b/libsrc/core/array.hpp index 00859e55..35d1ab30 100644 --- a/libsrc/core/array.hpp +++ b/libsrc/core/array.hpp @@ -216,6 +216,10 @@ namespace ngcore template constexpr T IndexBASE () { return T(0); } + template + constexpr T IndexBASE (T ind) { return IndexBASE(); } + + class IndexFromEnd { diff --git a/libsrc/interface/writegmsh2.cpp b/libsrc/interface/writegmsh2.cpp index 845556bd..74843195 100644 --- a/libsrc/interface/writegmsh2.cpp +++ b/libsrc/interface/writegmsh2.cpp @@ -58,7 +58,7 @@ namespace netgen int np = mesh.GetNP(); /// number of points in mesh int ne = mesh.GetNE(); /// number of 3D elements in mesh int nse = mesh.GetNSE(); /// number of surface elements (BC) - int i, j, k, l; + // int i, j, k, l; /* @@ -66,8 +66,8 @@ namespace netgen */ if ((ne > 0) - && (mesh.VolumeElement(1).GetNP() <= 10) - && (mesh.SurfaceElement(1).GetNP() <= 6)) + && (mesh.VolumeElements().First().GetNP() <= 10) + && (mesh.SurfaceElements().First().GetNP() <= 6)) { cout << "Write GMSH v2.xx Format \n"; cout << "The GMSH v2.xx export is currently available for elements upto 2nd Order\n" << endl; @@ -86,7 +86,7 @@ namespace netgen outfile << "$Nodes\n"; outfile << np << "\n"; - for (i = 1; i <= np; i++) + for (int i = 1; i <= np; i++) { const Point3d & p = mesh.Point(i); outfile << i << " "; /// node number @@ -101,13 +101,13 @@ namespace netgen outfile << "$Elements\n"; outfile << ne + nse << "\n"; //// number of elements + number of surfaces BC - for (i = 1; i <= nse; i++) + for (auto sei : Range(mesh.SurfaceElements())) { int elType = 0; - Element2d el = mesh.SurfaceElement(i); + Element2d el = mesh[sei]; // .SurfaceElement(i); if(invertsurf) el.Invert(); - + if(el.GetNP() == 3) elType = GMSH_TRIG; //// GMSH Type for a 3 node triangle if(el.GetNP() == 6) elType = GMSH_TRIG6; //// GMSH Type for a 6 node triangle if(elType == 0) @@ -116,7 +116,7 @@ namespace netgen return; } - outfile << i; + outfile << sei-IndexBASE(sei)+1; outfile << " "; outfile << elType; outfile << " "; @@ -125,7 +125,7 @@ namespace netgen outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " "; /// that means that physical entity = elementary entity (arbitrary approach) outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " "; - for (j = 1; j <= el.GetNP(); j++) + for (int j = 1; j <= el.GetNP(); j++) { outfile << " "; outfile << el.PNum(triGmsh[j]); @@ -133,12 +133,12 @@ namespace netgen outfile << "\n"; } - - for (i = 1; i <= ne; i++) + for (ElementIndex ei : Range(mesh.VolumeElements())) { + int i = ei-IndexBASE(ei)+1; int elType = 0; - Element el = mesh.VolumeElement(i); + Element el = mesh[ei]; if (inverttets) el.Invert(); if(el.GetNP() == 4) elType = GMSH_TET; //// GMSH Element type for 4 node tetrahedron @@ -160,7 +160,7 @@ namespace netgen outfile << " "; outfile << 100000 + el.GetIndex(); /// volume number outfile << " "; - for (j = 1; j <= el.GetNP(); j++) + for (int j = 1; j <= el.GetNP(); j++) { outfile << " "; outfile << el.PNum(tetGmsh[j]); @@ -193,7 +193,7 @@ namespace netgen outfile << "$Nodes\n"; outfile << np << "\n"; - for (i = 1; i <= np; i++) + for (int i = 1; i <= np; i++) { const Point3d & p = mesh.Point(i); outfile << i << " "; /// node number @@ -207,7 +207,7 @@ namespace netgen outfile << "$Elements\n"; outfile << nse << "\n"; - for (k = 1; k <= nse; k++) + for (int k = 1; k <= nse; k++) { int elType = 0; @@ -232,7 +232,7 @@ namespace netgen outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " "; /// that means that physical entity = elementary entity (arbitrary approach) outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " "; - for (l = 1; l <= el.GetNP(); l++) + for (int l = 1; l <= el.GetNP(); l++) { outfile << " "; if((elType == GMSH_TRIG) || (elType == GMSH_TRIG6)) diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp index 8043cb87..2870bb1d 100644 --- a/libsrc/meshing/bisect.cpp +++ b/libsrc/meshing/bisect.cpp @@ -5,10 +5,8 @@ #include "meshing.hpp" // quickfix for parallel - #define noDEBUG - namespace netgen { class MarkedTet diff --git a/libsrc/meshing/clusters.cpp b/libsrc/meshing/clusters.cpp index 2cba410c..c71bc733 100644 --- a/libsrc/meshing/clusters.cpp +++ b/libsrc/meshing/clusters.cpp @@ -143,19 +143,21 @@ namespace netgen cluster_reps.Elem(nnums[j]) = nnums[j]; } */ + + ngcore::ParallelForRange (mesh.SurfaceElements().Range(), [&] (auto myrange) { NgArrayMem nnums; // , ednums; - for (int i_ : myrange) + for (SurfaceElementIndex i_ : myrange) { int i = i_+1; - const Element2d & el = mesh.SurfaceElement(i); + const Element2d & el = mesh[i_]; // .SurfaceElement(i); ELEMENT_TYPE typ = el.GetType(); // top.GetSurfaceElementEdges (i, ednums); - auto ednums = top.GetEdges (SurfaceElementIndex(i_)); + auto ednums = top.GetEdges (i_); // cout << "ednums = " << ednums << endl; int fanum = top.GetSurfaceElementFace (i); diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp index a0811112..a021f8bc 100644 --- a/libsrc/meshing/improve3.cpp +++ b/libsrc/meshing/improve3.cpp @@ -357,7 +357,7 @@ void MeshOptimize3d :: CombineImprove () { static Timer t("MeshOptimize3d::CombineImprove"); RegionTimer reg(t); static Timer topt("Optimize"); - static Timer tsearch("Search"); + static Timer tsearch("Search-combine"); static Timer tbuild_elements_table("Build elements table"); mesh.BuildBoundaryEdges(false); @@ -613,7 +613,7 @@ void MeshOptimize3d :: SplitImprove () { static Timer t("MeshOptimize3d::SplitImprove"); RegionTimer reg(t); static Timer topt("Optimize"); - static Timer tsearch("Search"); + static Timer tsearch("Search-split"); // int np = mesh.GetNP(); int ne = mesh.GetNE(); diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index be26bf80..16ee4a0a 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -3051,6 +3051,7 @@ namespace netgen // bool buggy = false; // ofstream bout("buggy.out"); + for (int i = 1; i <= GetNSE(); i++) { const Element2d & el = SurfaceElement(i); @@ -3060,10 +3061,10 @@ namespace netgen { for (int j = 1; j <= el.GetNP(); j++) { - INDEX_3 seg (el.PNumMod(j), el.PNumMod(j+1), el.GetIndex()); + PointIndices<3> seg (el.PNumMod(j), el.PNumMod(j+1), el.GetIndex()); // int data; - if (seg.I1() < PointIndex::BASE || seg.I2() < PointIndex::BASE) + if (!seg.I1().IsValid() || !seg.I2().IsValid()) cerr << "seg = " << seg << endl; if (faceht.Used(seg)) @@ -7299,6 +7300,7 @@ namespace netgen { int eli0, eli1; GetTopology().GetSurface2VolumeElement(sei+1, eli0, eli1); + // auto [ei0,ei1] = GetTopology().GetSurface2VolumeElement(sei); // the way to go auto & sel = (*this)[sei]; int face = sel.GetIndex(); int domin = VolumeElement(eli0).GetIndex(); diff --git a/libsrc/meshing/meshing3.cpp b/libsrc/meshing/meshing3.cpp index 57e19d49..0ed40bb1 100644 --- a/libsrc/meshing/meshing3.cpp +++ b/libsrc/meshing/meshing3.cpp @@ -9,8 +9,6 @@ double minwithoutother; - - MeshingStat3d :: MeshingStat3d () { cntsucc = cnttrials = cntelem = qualclass = 0; diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 843ba9ca..1fa2ed47 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -27,7 +27,6 @@ namespace netgen */ - enum ELEMENT_TYPE : unsigned char { SEGMENT = 1, SEGMENT3 = 2, TRIG = 10, QUAD=11, TRIG6 = 12, QUAD6 = 13, QUAD8 = 14, @@ -149,8 +148,107 @@ namespace netgen } + template + class Index + { + public: + + T i; + static constexpr int BASE = Base; + + public: + class t_invalid { public: constexpr t_invalid() = default; }; + static constexpr t_invalid INVALID{}; + + constexpr Index () = default; + constexpr Index (const Index&) = default; + constexpr Index (Index &&) = default; + Index & operator= (const Index&) = default; + Index & operator= (Index&&) = default; + + // private: + constexpr Index (int ai) : i(ai) + { +#ifdef DEBUG + if (ai < Base) + cout << "illegal Index, use Index::INVALID instead" << endl; +#endif + } + // friend constexpr netgen::TIndex ngcore::IndexBASE (); + // template friend class PointIndices; + + /* + friend auto operator+ (Index, int) -> TIndex; + friend TIndex operator+ (Index, size_t); + friend TIndex operator+ (int, Index); + friend TIndex operator+ (size_t, Index); + friend constexpr TIndex operator- (Index, int); + friend int operator- (Index, Index); + friend bool operator< (Index a, Index b); + friend bool operator> (Index a, Index b); + friend bool operator>= (Index a, Index b); + friend bool operator<= (Index a, Index b); + friend bool operator== (Index a, Index b); + friend bool operator!= (Index a, Index b); + */ + + public: + constexpr Index (t_invalid inv) : i(long(BASE)-1) { ; } + // TIndex & operator= (const TIndex &ai) { i = ai.i; return *this; } + // private: + constexpr operator const int& () const { return i; } + explicit constexpr operator int& () { return i; } + public: + constexpr operator TIndex() const { return TIndex(i); } + constexpr operator TIndex&() { return static_cast(*this); } + TIndex operator++ (int) { TIndex hi(*this); i++; return hi; } + TIndex operator-- (int) { TIndex hi(*this); i--; return hi; } + TIndex & operator++ () { i++; return *this; } + TIndex operator-- () { i--; return *this; } + TIndex operator+= (int add) { i += add; return *this; } + void Invalidate() { i = long(TIndex::BASE)-1; } + bool IsValid() const { return i+1 != TIndex::BASE; } + // operator bool() const { return IsValid(); } + + void DoArchive (Archive & ar) { ar & i; } + }; + + + /* + template + auto operator+ (Index pi, int i) -> TIndex { return TIndex(pi.i+i); } + template + inline TIndex operator+ (Index pi, size_t i) { return TIndex(pi.i+i); } + template + inline TIndex operator+ (int i, Index pi) { return TIndex(pi.i+i); } + template + inline TIndex operator+ (size_t i, Index pi) { return TIndex(pi.i+i); } + + template + constexpr inline auto operator- (Index pi, int i) -> TIndex { return TIndex(pi.i-i); } + + template + inline int operator- (Index pa, Index pb) { return pa.i-pb.i; } + template + inline bool operator< (Index a, Index b) { return a.i-b.i < 0; } + template + inline bool operator> (Index a, Index b) { return a.i-b.i > 0; } + template + inline bool operator>= (Index a, Index b) { return a.i-b.i >= 0; } + template + inline bool operator<= (Index a, Index b) { return a.i-b.i <= 0; } + template + inline bool operator== (Index a, Index b) { return a.i == b.i; } + template + inline bool operator!= (Index a, Index b) { return a.i != b.i; } + */ + + + + + /* class PointIndex { @@ -196,7 +294,9 @@ namespace netgen // #define BASE0 - + + + /* class PointIndex { int i; @@ -260,6 +360,14 @@ namespace netgen void DoArchive (Archive & ar) { ar & i; } }; + */ + + class PointIndex : public Index + { + public: + using Index::Index; + }; + constexpr inline PointIndex operator+ (PointIndex pi, int i) { return PointIndex(pi.i+i); } constexpr inline PointIndex operator+ (PointIndex pi, size_t i) { return PointIndex(pi.i+i); } constexpr inline PointIndex operator+ (int i, PointIndex pi) { return PointIndex(pi.i+i); } @@ -501,8 +609,8 @@ namespace netgen inline bool operator> (ElementIndex ei1, ElementIndex ei2) { return int(ei1) > int(ei2); }; inline bool operator>= (ElementIndex ei1, ElementIndex ei2) { return int(ei1) >= int(ei2); }; inline bool operator<= (ElementIndex ei1, ElementIndex ei2) { return int(ei1) <= int(ei2); }; - // these should not be needed: + // these should not be needed: inline bool operator== (ElementIndex ei1, int ei2) { return int(ei1) == int(ei2); }; inline bool operator< (size_t s, ElementIndex ei2) { return int(s) < int(ei2); }; inline bool operator< (ElementIndex ei1, size_t s) { return int(ei1) < int(s); }; // should not need @@ -530,6 +638,27 @@ namespace netgen SurfaceElementIndex & operator+= (int inc) { i+=inc; return *this; } void DoArchive (Archive & ar) { ar & i; } }; + + inline SurfaceElementIndex operator+ (SurfaceElementIndex ei, int i) { return SurfaceElementIndex { int(ei) + i }; } + inline SurfaceElementIndex operator+ (size_t s, SurfaceElementIndex ei) { return SurfaceElementIndex(int(ei) + s); } + inline SurfaceElementIndex operator+ (SurfaceElementIndex ei, size_t s) { return SurfaceElementIndex(int(ei) + s); } + inline bool operator== (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) == int(ei2); }; + inline bool operator!= (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) != int(ei2); }; + inline bool operator< (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) < int(ei2); }; + inline bool operator> (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) > int(ei2); }; + inline bool operator>= (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) >= int(ei2); }; + inline bool operator<= (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) <= int(ei2); }; + + // these should not be needed: + inline bool operator== (SurfaceElementIndex ei1, int ei2) { return int(ei1) == int(ei2); }; + inline bool operator== (int ei2, SurfaceElementIndex ei1) { return int(ei1) == int(ei2); }; + inline bool operator!= (SurfaceElementIndex ei1, int ei2) { return int(ei1) != int(ei2); }; + inline bool operator< (size_t s, SurfaceElementIndex ei2) { return int(s) < int(ei2); }; + inline bool operator< (SurfaceElementIndex ei1, size_t s) { return int(ei1) < int(s); }; // should not need + inline bool operator< (SurfaceElementIndex ei1, int s) { return int(ei1) < int(s); }; // should not need + inline bool operator>= (size_t s, SurfaceElementIndex ei2) { return int(s) >= int(ei2); }; + inline bool operator>= (SurfaceElementIndex ei1, int s) { return int(ei1) >= int(s); }; + inline void SetInvalid (SurfaceElementIndex & id) { id = -1; } inline bool IsInvalid (SurfaceElementIndex & id) { return id == -1; } diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp index 4d7b9a1d..7745f33a 100644 --- a/libsrc/meshing/parallelmesh.cpp +++ b/libsrc/meshing/parallelmesh.cpp @@ -592,18 +592,17 @@ namespace netgen } DynamicTable elementarrays(elarraysize); - - for (int ei = 1; ei <= GetNE(); ei++) + + for (ElementIndex ei : VolumeElements().Range()) { const Element & el = VolumeElement (ei); - // int dest = el.GetPartition(); - int dest = vol_partition[ei-1]; + int dest = vol_partition[ei]; - elementarrays.Add (dest, ei); + elementarrays.Add (dest, int(ei+1)); elementarrays.Add (dest, el.GetIndex()); elementarrays.Add (dest, el.GetNP()); - for (int i = 0; i < el.GetNP(); i++) - elementarrays.Add (dest, el[i]); + for (PointIndex pi : el.PNums()) + elementarrays.Add (dest, pi); } tbuildelementtable.Stop(); diff --git a/libsrc/meshing/refine.cpp b/libsrc/meshing/refine.cpp index 9c5f5556..dbbfec10 100644 --- a/libsrc/meshing/refine.cpp +++ b/libsrc/meshing/refine.cpp @@ -865,9 +865,10 @@ namespace netgen for (int k = 1; k <= 3; k++) { fhelp.Clear(); - for (int i = 1; i <= mesh.GetNE(); i++) + // for (int i = 1; i <= mesh.GetNE(); i++) + for (const Element & el : mesh.VolumeElements()) { - const Element & el = mesh.VolumeElement(i); + // const Element & el = mesh.VolumeElement(i); int freeel = 0; for (int j = 1; j <= el.GetNP(); j++) if (free.Test(el.PNum(j))) @@ -897,19 +898,19 @@ namespace netgen wrongels = 0; - for (int i = 1; i <= mesh.GetNE(); i++) + for (ElementIndex ei : mesh.VolumeElements().Range()) { - if (mesh.VolumeElement(i).Volume(mesh.Points()) < 0) + if (mesh.VolumeElement(ei).Volume(mesh.Points()) < 0) { wrongels++; - mesh.VolumeElement(i).Flags().badel = 1; + mesh.VolumeElement(ei).Flags().badel = 1; (*testout) << "wrong el: "; for (int j = 1; j <= 4; j++) - (*testout) << mesh.VolumeElement(i).PNum(j) << " "; + (*testout) << mesh.VolumeElement(ei).PNum(j) << " "; (*testout) << endl; } else - mesh.VolumeElement(i).Flags().badel = 0; + mesh.VolumeElement(ei).Flags().badel = 0; } cout << "wrongels = " << wrongels << endl; } diff --git a/libsrc/meshing/topology.hpp b/libsrc/meshing/topology.hpp index fcb1b967..782d1ba6 100644 --- a/libsrc/meshing/topology.hpp +++ b/libsrc/meshing/topology.hpp @@ -185,6 +185,13 @@ public: elnr2 = surf2volelement.Get(selnr)[1]; } + std::array GetSurface2VolumeElement (SurfaceElementIndex sei) + { + return { ElementIndex( surf2volelement.Get(sei+1)[0] - 1), + ElementIndex( surf2volelement.Get(sei+1)[1] - 1) }; + } + + int GetFace2SurfaceElement (int fnr) const { return face2surfel[fnr-1]; } SegmentIndex GetSegmentOfEdge(int edgenr) const { return edge2segment[edgenr-1]; } diff --git a/libsrc/occ/occmeshsurf.cpp b/libsrc/occ/occmeshsurf.cpp index 37f23cfc..16a0e427 100644 --- a/libsrc/occ/occmeshsurf.cpp +++ b/libsrc/occ/occmeshsurf.cpp @@ -366,8 +366,11 @@ namespace netgen void OCCSurface :: Project (Point<3> & ap, PointGeomInfo & gi) { - // static Timer t("OccSurface::Project"); RegionTimer reg(t); - // static Timer t2("OccSurface::Project actual"); + static Timer t("OccSurface::Project"); RegionTimer reg(t); + static Timer tanal("OccSurface::Project analysis"); + static Timer ttol("OccSurface::Project approximation"); + + static Timer t2("OccSurface::Project actual"); // try Newton's method ... @@ -472,14 +475,18 @@ namespace netgen */ // double u,v; - // JS : shouldn't we move these 2 lines to the constructor ? + // JS : shouldn't we move these 2 lines to the constructor ? + // tanal.Start(); Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface ); + // ShapeAnalysis_Surface su( occface ); + // tanal.Stop(); + ttol.Start(); auto toltool = BRep_Tool::Tolerance( topods_face ); - + ttol.Stop(); // gp_Pnt2d suval = su->ValueOfUV ( pnt, toltool); - // t2.Start(); + t2.Start(); gp_Pnt2d suval = su->NextValueOfUV (gp_Pnt2d(u,v), pnt, toltool); - // t2.Stop(); + t2.Stop(); suval.Coord( u, v); pnt = occface->Value( u, v ); diff --git a/libsrc/occ/occmeshsurf.hpp b/libsrc/occ/occmeshsurf.hpp index b4c4b06e..c045da69 100644 --- a/libsrc/occ/occmeshsurf.hpp +++ b/libsrc/occ/occmeshsurf.hpp @@ -5,10 +5,13 @@ #include "occgeom.hpp" #include "mydefs.hpp" - +#include #include #include #include +#include +#include + #define PARAMETERSPACE -1 #define PLANESPACE 1 @@ -30,7 +33,8 @@ public: Handle(Geom_Surface) occface; TopAbs_Orientation orient; int projecttype; - + ShapeAnalysis_Surface su; + Standard_Real toltool; protected: Point<3> p1; Point<3> p2; @@ -60,10 +64,15 @@ protected: public: OCCSurface (const TopoDS_Face & aface, int aprojecttype) + : topods_face(aface), + occface(BRep_Tool::Surface(topods_face)), + su( occface ), + toltool(BRep_Tool::Tolerance(topods_face)) + { static Timer t("occurface ctor"); RegionTimer r(t); topods_face = aface; - occface = BRep_Tool::Surface(topods_face); + // occface = BRep_Tool::Surface(topods_face); orient = topods_face.Orientation(); projecttype = aprojecttype; ShapeAnalysis::GetFaceUVBounds (topods_face, umin, umax, vmin, vmax); @@ -72,6 +81,8 @@ public: umax += fabs(umax-umin)/100.0; vmax += fabs(vmax-vmin)/100.0; // projecttype = PLANESPACE; + + // su = ShapeAnalysis_Surface( occface ); /* TopExp_Explorer exp1; exp1.Init (topods_face, TopAbs_WIRE); diff --git a/libsrc/occ/occpkg.cpp b/libsrc/occ/occpkg.cpp index ecc11af8..57b1b78e 100644 --- a/libsrc/occ/occpkg.cpp +++ b/libsrc/occ/occpkg.cpp @@ -817,21 +817,29 @@ namespace netgen if(strcmp(argv[1], "showall") == 0) { + /* for(int i = 1; i <= mesh->GetNSE(); i++) { - mesh->SurfaceElement(i).Visible(1); + mesh->SurfaceElement(i).Visible(1); } - - mesh->SetNextTimeStamp(); + */ + for (auto & el : mesh->SurfaceElements()) + el.Visible(1); + + mesh->SetNextTimeStamp(); } if(strcmp(argv[1], "hideall") == 0) { + /* for(int i = 1; i <= mesh->GetNSE(); i++) { - mesh->SurfaceElement(i).Visible(0); + mesh->SurfaceElement(i).Visible(0); } - + */ + for (auto & el : mesh->SurfaceElements()) + el.Visible(0); + mesh->SetNextTimeStamp(); } diff --git a/nglib/nglib.cpp b/nglib/nglib.cpp index 11a43ddb..a73fb324 100644 --- a/nglib/nglib.cpp +++ b/nglib/nglib.cpp @@ -270,7 +270,7 @@ namespace nglib NGLIB_API Ng_Surface_Element_Type Ng_GetSurfaceElement (Ng_Mesh * mesh, int num, int * pi) { - const Element2d & el = ((Mesh*)mesh)->SurfaceElement(num); + const Element2d & el = ((Mesh*)mesh)->SurfaceElement(SurfaceElementIndex(num-1)); for (int i = 1; i <= el.GetNP(); i++) pi[i-1] = el.PNum(i); Ng_Surface_Element_Type et; @@ -301,7 +301,7 @@ namespace nglib NGLIB_API Ng_Volume_Element_Type Ng_GetVolumeElement (Ng_Mesh * mesh, int num, int * pi) { - const Element & el = ((Mesh*)mesh)->VolumeElement(num); + const Element & el = ((Mesh*)mesh)->VolumeElement(ElementIndex(num-1)); for (int i = 1; i <= el.GetNP(); i++) pi[i-1] = el.PNum(i); Ng_Volume_Element_Type et; @@ -439,7 +439,7 @@ namespace nglib NGLIB_API Ng_Surface_Element_Type Ng_GetElement_2D (Ng_Mesh * mesh, int num, int * pi, int * matnum) { - const Element2d & el = ((Mesh*)mesh)->SurfaceElement(num); + const Element2d & el = ((Mesh*)mesh)->SurfaceElement(SurfaceElementIndex(num-1)); for (int i = 1; i <= el.GetNP(); i++) pi[i-1] = el.PNum(i); From 9ab086f81941e611711371b8f9e848e1e68aed1e Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 31 Dec 2024 12:54:25 +0100 Subject: [PATCH 02/19] unified Index class --- libsrc/meshing/meshtype.hpp | 258 +++++---------------------------- libsrc/meshing/python_mesh.cpp | 2 +- 2 files changed, 38 insertions(+), 222 deletions(-) diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 1fa2ed47..04921418 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -198,16 +198,18 @@ namespace netgen constexpr Index (t_invalid inv) : i(long(BASE)-1) { ; } // TIndex & operator= (const TIndex &ai) { i = ai.i; return *this; } // private: - constexpr operator const int& () const { return i; } - explicit constexpr operator int& () { return i; } + constexpr operator T () const { return i; } + explicit constexpr operator T& () { return i; } public: constexpr operator TIndex() const { return TIndex(i); } constexpr operator TIndex&() { return static_cast(*this); } + TIndex operator++ (int) { TIndex hi(*this); i++; return hi; } TIndex operator-- (int) { TIndex hi(*this); i--; return hi; } TIndex & operator++ () { i++; return *this; } TIndex operator-- () { i--; return *this; } - TIndex operator+= (int add) { i += add; return *this; } + constexpr TIndex & operator+= (T add) & { i += add; return *this; } + constexpr TIndex operator+= (T add) && { i += add; return *this; } void Invalidate() { i = long(TIndex::BASE)-1; } bool IsValid() const { return i+1 != TIndex::BASE; } // operator bool() const { return IsValid(); } @@ -216,19 +218,17 @@ namespace netgen }; - /* template - auto operator+ (Index pi, int i) -> TIndex { return TIndex(pi.i+i); } + // constexpr auto operator+ (Index ind, int i) { return TIndex(T(ind)+i); } + constexpr auto operator+ (Index ind, int i) { return TIndex(ind)+=i; } template - inline TIndex operator+ (Index pi, size_t i) { return TIndex(pi.i+i); } + constexpr TIndex operator+ (Index pi, size_t i) { return TIndex(pi.i+i); } template - inline TIndex operator+ (int i, Index pi) { return TIndex(pi.i+i); } + constexpr TIndex operator+ (int i, Index pi) { return TIndex(pi.i+i); } template inline TIndex operator+ (size_t i, Index pi) { return TIndex(pi.i+i); } - template constexpr inline auto operator- (Index pi, int i) -> TIndex { return TIndex(pi.i-i); } - template inline int operator- (Index pa, Index pb) { return pa.i-pb.i; } template @@ -243,143 +243,15 @@ namespace netgen inline bool operator== (Index a, Index b) { return a.i == b.i; } template inline bool operator!= (Index a, Index b) { return a.i != b.i; } - */ - - - - /* - class PointIndex - { - int i; - public: - class t_invalid { public: constexpr t_invalid() = default; }; - static constexpr t_invalid INVALID{}; - - PointIndex () = default; - PointIndex (const PointIndex&) = default; - PointIndex (PointIndex &&) = default; - PointIndex & operator= (const PointIndex&) = default; - PointIndex & operator= (PointIndex&&) = default; - - constexpr PointIndex (int ai) : i(ai) - { -#ifdef DEBUG - if (ai < PointIndex::BASE) - cout << "illegal PointIndex, use PointIndex::INVALID instead" << endl; - // throw Exception("illegal PointIndex, use PointIndex::INVALID instead"); -#endif - } - constexpr PointIndex (t_invalid inv) : i(PointIndex::BASE-1) { ; } - // PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; } - constexpr operator const int& () const { return i; } - explicit constexpr operator int& () { return i; } - PointIndex operator++ (int) { PointIndex hi(*this); i++; return hi; } - PointIndex operator-- (int) { PointIndex hi(*this); i--; return hi; } - PointIndex & operator++ () { i++; return *this; } - PointIndex operator-- () { i--; return *this; } - PointIndex operator+= (int add) { i += add; return *this; } - void Invalidate() { i = PointIndex::BASE-1; } - bool IsValid() const { return i != PointIndex::BASE-1; } -#ifdef BASE0 - static constexpr size_t BASE = 0; -#else - static constexpr size_t BASE = 1; -#endif - - void DoArchive (Archive & ar) { ar & i; } - } - */ - - - // #define BASE0 - - - /* - class PointIndex - { - int i; - public: - class t_invalid { public: constexpr t_invalid() = default; }; - static constexpr t_invalid INVALID{}; - - PointIndex () = default; - PointIndex (const PointIndex&) = default; - PointIndex (PointIndex &&) = default; - PointIndex & operator= (const PointIndex&) = default; - PointIndex & operator= (PointIndex&&) = default; - - // private: - constexpr PointIndex (int ai) : i(ai) - { -#ifdef DEBUG - if (ai < PointIndex::BASE) - cout << "illegal PointIndex, use PointIndex::INVALID instead" << endl; - // throw Exception("illegal PointIndex, use PointIndex::INVALID instead"); -#endif - } - - friend constexpr netgen::PointIndex ngcore::IndexBASE (); - template friend class PointIndices; - - friend constexpr PointIndex operator+ (PointIndex, int); - friend constexpr PointIndex operator+ (PointIndex, size_t); - friend constexpr PointIndex operator+ (int, PointIndex); - friend constexpr PointIndex operator+ (size_t, PointIndex); - friend constexpr PointIndex operator- (PointIndex, int); - friend constexpr int operator- (PointIndex, PointIndex); - friend bool operator< (PointIndex a, PointIndex b); - friend bool operator> (PointIndex a, PointIndex b); - friend bool operator>= (PointIndex a, PointIndex b); - friend bool operator<= (PointIndex a, PointIndex b); - friend constexpr bool operator== (PointIndex a, PointIndex b); - friend constexpr bool operator!= (PointIndex a, PointIndex b); - - public: - constexpr PointIndex (t_invalid inv) : i(long(PointIndex::BASE)-1) { ; } - // PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; } - // private: - constexpr operator const int& () const { return i; } - explicit constexpr operator int& () { return i; } - public: - PointIndex operator++ (int) { PointIndex hi(*this); i++; return hi; } - PointIndex operator-- (int) { PointIndex hi(*this); i--; return hi; } - PointIndex & operator++ () { i++; return *this; } - PointIndex operator-- () { i--; return *this; } - PointIndex operator+= (int add) { i += add; return *this; } - void Invalidate() { i = long(PointIndex::BASE)-1; } - bool IsValid() const { return i+1 != PointIndex::BASE; } - // operator bool() const { return IsValid(); } -#ifdef BASE0 - static constexpr size_t BASE = 0; -#else - static constexpr size_t BASE = 1; -#endif - - void DoArchive (Archive & ar) { ar & i; } - }; - - */ class PointIndex : public Index { public: using Index::Index; }; - - constexpr inline PointIndex operator+ (PointIndex pi, int i) { return PointIndex(pi.i+i); } - constexpr inline PointIndex operator+ (PointIndex pi, size_t i) { return PointIndex(pi.i+i); } - constexpr inline PointIndex operator+ (int i, PointIndex pi) { return PointIndex(pi.i+i); } - constexpr inline PointIndex operator+ (size_t i, PointIndex pi) { return PointIndex(pi.i+i); } - constexpr inline PointIndex operator- (PointIndex pi, int i) { return PointIndex(pi.i-i); } - constexpr inline int operator- (PointIndex pa, PointIndex pb) { return pa.i-pb.i; } - inline bool operator< (PointIndex a, PointIndex b) { return a.i-b.i < 0; } - inline bool operator> (PointIndex a, PointIndex b) { return a.i-b.i > 0; } - inline bool operator>= (PointIndex a, PointIndex b) { return a.i-b.i >= 0; } - inline bool operator<= (PointIndex a, PointIndex b) { return a.i-b.i <= 0; } - inline constexpr bool operator== (PointIndex a, PointIndex b) { return a.i == b.i; } - inline constexpr bool operator!= (PointIndex a, PointIndex b) { return a.i != b.i; } + } namespace ngcore @@ -396,16 +268,15 @@ namespace netgen { // int i; ist >> i; pi = PointIndex(i); return ist; int i; ist >> i; - // pi = PointIndex(i); pi = IndexBASE()+i-1; return ist; } inline ostream & operator<< (ostream & ost, const PointIndex & pi) { + // return (ost << int(pi)); int intpi = pi - IndexBASE() + 1; return (ost << intpi); - // return (ost << int(pi)); } template class PointIndices; @@ -571,25 +442,12 @@ namespace std namespace netgen { - class ElementIndex + class ElementIndex : public Index { - int i; public: - ElementIndex () = default; - constexpr /* explicit */ ElementIndex (int ai) : i(ai) { ; } - ElementIndex & operator= (const ElementIndex & ai) { i = ai.i; return *this; } - ElementIndex & operator= (int ai) { i = ai; return *this; } - - constexpr /* explicit */ operator int () const { return i; } - ElementIndex operator++ (int) { return ElementIndex(i++); } - ElementIndex operator-- (int) { return ElementIndex(i--); } - ElementIndex & operator++ () { ++i; return *this; } - ElementIndex & operator-- () { --i; return *this; } - - int operator- (ElementIndex ei2) const { return i-ei2.i; } - friend constexpr ElementIndex ngcore::IndexBASE(); + using Index::Index; }; - + inline istream & operator>> (istream & ist, ElementIndex & pi) { int i; ist >> i; pi = i; return ist; @@ -600,64 +458,31 @@ namespace netgen return (ost << ei-IndexBASE()); } - inline ElementIndex operator+ (ElementIndex ei, int i) { return ElementIndex { int(ei) + i }; } - inline ElementIndex operator+ (size_t s, ElementIndex ei) { return ElementIndex(int(ei) + s); } - inline ElementIndex operator+ (ElementIndex ei, size_t s) { return ElementIndex(int(ei) + s); } - inline bool operator== (ElementIndex ei1, ElementIndex ei2) { return int(ei1) == int(ei2); }; - inline bool operator!= (ElementIndex ei1, ElementIndex ei2) { return int(ei1) != int(ei2); }; - inline bool operator< (ElementIndex ei1, ElementIndex ei2) { return int(ei1) < int(ei2); }; - inline bool operator> (ElementIndex ei1, ElementIndex ei2) { return int(ei1) > int(ei2); }; - inline bool operator>= (ElementIndex ei1, ElementIndex ei2) { return int(ei1) >= int(ei2); }; - inline bool operator<= (ElementIndex ei1, ElementIndex ei2) { return int(ei1) <= int(ei2); }; - - // these should not be needed: - inline bool operator== (ElementIndex ei1, int ei2) { return int(ei1) == int(ei2); }; - inline bool operator< (size_t s, ElementIndex ei2) { return int(s) < int(ei2); }; - inline bool operator< (ElementIndex ei1, size_t s) { return int(ei1) < int(s); }; // should not need - inline bool operator< (ElementIndex ei1, int s) { return int(ei1) < int(s); }; // should not need - inline bool operator>= (size_t s, ElementIndex ei2) { return int(s) >= int(ei2); }; - - class SurfaceElementIndex + // these should not be needed soon + inline bool operator== (Index ei1, int ei2) { return int(ei1) == int(ei2); }; + inline bool operator< (size_t s, Index ei2) { return int(s) < int(ei2); }; + inline bool operator< ( Index ei1, size_t s) { return int(ei1) < int(s); }; // should not need + inline bool operator< ( Index ei1, int s) { return int(ei1) < int(s); }; // should not need + inline bool operator>= (size_t s, Index ei2) { return int(s) >= int(ei2); }; + + + class SurfaceElementIndex : public Index { - int i; public: - SurfaceElementIndex () = default; - constexpr SurfaceElementIndex (int ai) : i(ai) { ; } - /* - SurfaceElementIndex & operator= (const SurfaceElementIndex & ai) - { i = ai.i; return *this; } - */ - SurfaceElementIndex & operator= (const SurfaceElementIndex & ai) = default; - SurfaceElementIndex & operator= (int ai) { i = ai; return *this; } - constexpr operator int () const { return i; } - SurfaceElementIndex operator++ (int) { SurfaceElementIndex hi(*this); i++; return hi; } - SurfaceElementIndex operator-- (int) { SurfaceElementIndex hi(*this); i--; return hi; } - SurfaceElementIndex & operator++ () { ++i; return *this; } - SurfaceElementIndex & operator-- () { --i; return *this; } - SurfaceElementIndex & operator+= (int inc) { i+=inc; return *this; } - void DoArchive (Archive & ar) { ar & i; } + using Index::Index; }; - inline SurfaceElementIndex operator+ (SurfaceElementIndex ei, int i) { return SurfaceElementIndex { int(ei) + i }; } - inline SurfaceElementIndex operator+ (size_t s, SurfaceElementIndex ei) { return SurfaceElementIndex(int(ei) + s); } - inline SurfaceElementIndex operator+ (SurfaceElementIndex ei, size_t s) { return SurfaceElementIndex(int(ei) + s); } - inline bool operator== (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) == int(ei2); }; - inline bool operator!= (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) != int(ei2); }; - inline bool operator< (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) < int(ei2); }; - inline bool operator> (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) > int(ei2); }; - inline bool operator>= (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) >= int(ei2); }; - inline bool operator<= (SurfaceElementIndex ei1, SurfaceElementIndex ei2) { return int(ei1) <= int(ei2); }; - - // these should not be needed: - inline bool operator== (SurfaceElementIndex ei1, int ei2) { return int(ei1) == int(ei2); }; - inline bool operator== (int ei2, SurfaceElementIndex ei1) { return int(ei1) == int(ei2); }; - inline bool operator!= (SurfaceElementIndex ei1, int ei2) { return int(ei1) != int(ei2); }; - inline bool operator< (size_t s, SurfaceElementIndex ei2) { return int(s) < int(ei2); }; - inline bool operator< (SurfaceElementIndex ei1, size_t s) { return int(ei1) < int(s); }; // should not need - inline bool operator< (SurfaceElementIndex ei1, int s) { return int(ei1) < int(s); }; // should not need - inline bool operator>= (size_t s, SurfaceElementIndex ei2) { return int(s) >= int(ei2); }; - inline bool operator>= (SurfaceElementIndex ei1, int s) { return int(ei1) >= int(s); }; + + // these should not be needed soon + inline bool operator== (Index ei1, int ei2) { return int(ei1) == int(ei2); }; + inline bool operator== (int ei2, Index ei1) { return int(ei1) == int(ei2); }; + inline bool operator!= (Index ei1, int ei2) { return int(ei1) != int(ei2); }; + inline bool operator< (size_t s, Index ei2) { return int(s) < int(ei2); }; + inline bool operator< (Index ei1, size_t s) { return int(ei1) < int(s); }; // should not need + inline bool operator< (Index ei1, int s) { return int(ei1) < int(s); }; // should not need + inline bool operator>= (size_t s, Index ei2) { return int(s) >= int(ei2); }; + inline bool operator>= (Index ei1, int s) { return int(ei1) >= int(s); }; inline void SetInvalid (SurfaceElementIndex & id) { id = -1; } @@ -673,20 +498,11 @@ namespace netgen return (ost << int(pi)); } - class SegmentIndex + + class SegmentIndex : public Index { - int i; public: - SegmentIndex () = default; - constexpr SegmentIndex (int ai) : i(ai) { ; } - SegmentIndex & operator= (const SegmentIndex & ai) - { i = ai.i; return *this; } - SegmentIndex & operator= (int ai) { i = ai; return *this; } - constexpr operator int () const { return i; } - SegmentIndex& operator++ () { ++i; return *this; } - SegmentIndex& operator-- () { --i; return *this; } - SegmentIndex operator++ (int) { return i++; } - SegmentIndex operator-- (int) { return i--; } + using Index::Index; }; inline void SetInvalid (SegmentIndex & id) { id = -1; } diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 3faa7d51..1a133165 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -268,7 +268,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def(py::init()) .def("__repr__", &ToString) .def("__str__", &ToString) - .def_property_readonly("nr", &PointIndex::operator const int&) + .def_property_readonly("nr", &PointIndex::operator int) .def("__eq__" , FunctionPointer( [](PointIndex &self, PointIndex &other) { return static_cast(self)==static_cast(other); }) ) .def("__hash__" , FunctionPointer( [](PointIndex &self ) { return static_cast(self); }) ) From 626507f8fb4cc4ff0b246eb6c17c78b5e06df417 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 31 Dec 2024 13:03:46 +0100 Subject: [PATCH 03/19] missing constexpr --- libsrc/meshing/meshtype.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 04921418..289d064a 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -204,12 +204,12 @@ namespace netgen constexpr operator TIndex() const { return TIndex(i); } constexpr operator TIndex&() { return static_cast(*this); } - TIndex operator++ (int) { TIndex hi(*this); i++; return hi; } + TIndex operator++ (int) { TIndex hi{*this}; i++; return hi; } TIndex operator-- (int) { TIndex hi(*this); i--; return hi; } TIndex & operator++ () { i++; return *this; } TIndex operator-- () { i--; return *this; } constexpr TIndex & operator+= (T add) & { i += add; return *this; } - constexpr TIndex operator+= (T add) && { i += add; return *this; } + constexpr TIndex operator+= (T add) && { i += add; return TIndex(*this); } void Invalidate() { i = long(TIndex::BASE)-1; } bool IsValid() const { return i+1 != TIndex::BASE; } // operator bool() const { return IsValid(); } @@ -230,7 +230,7 @@ namespace netgen template constexpr inline auto operator- (Index pi, int i) -> TIndex { return TIndex(pi.i-i); } template - inline int operator- (Index pa, Index pb) { return pa.i-pb.i; } + constexpr inline int operator- (Index pa, Index pb) { return pa.i-pb.i; } template inline bool operator< (Index a, Index b) { return a.i-b.i < 0; } template From 7fac77d28e7031d76fec987d202b4e428d960558 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 31 Dec 2024 13:04:57 +0100 Subject: [PATCH 04/19] calling ctor --- libsrc/meshing/meshtype.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 289d064a..2662a5a2 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -220,7 +220,7 @@ namespace netgen template // constexpr auto operator+ (Index ind, int i) { return TIndex(T(ind)+i); } - constexpr auto operator+ (Index ind, int i) { return TIndex(ind)+=i; } + constexpr auto operator+ (Index ind, int i) { return TIndex{ind}+=i; } template constexpr TIndex operator+ (Index pi, size_t i) { return TIndex(pi.i+i); } template From 990fb0657cc8ef9a14b957c02199b79c159aeb23 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 31 Dec 2024 13:09:33 +0100 Subject: [PATCH 05/19] calling ctor --- libsrc/meshing/meshtype.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 2662a5a2..6be61848 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -220,7 +220,7 @@ namespace netgen template // constexpr auto operator+ (Index ind, int i) { return TIndex(T(ind)+i); } - constexpr auto operator+ (Index ind, int i) { return TIndex{ind}+=i; } + constexpr auto operator+ (Index ind, int i) { return TIndex( Index(ind) +=i ); } template constexpr TIndex operator+ (Index pi, size_t i) { return TIndex(pi.i+i); } template From 3c273bf5376d3eda8914ddd3386034cb9ceeece8 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 31 Dec 2024 13:13:57 +0100 Subject: [PATCH 06/19] calling ctor --- libsrc/meshing/meshtype.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 6be61848..08b3c87f 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -202,14 +202,14 @@ namespace netgen explicit constexpr operator T& () { return i; } public: constexpr operator TIndex() const { return TIndex(i); } - constexpr operator TIndex&() { return static_cast(*this); } + operator TIndex&() { return static_cast(*this); } TIndex operator++ (int) { TIndex hi{*this}; i++; return hi; } TIndex operator-- (int) { TIndex hi(*this); i--; return hi; } TIndex & operator++ () { i++; return *this; } TIndex operator-- () { i--; return *this; } constexpr TIndex & operator+= (T add) & { i += add; return *this; } - constexpr TIndex operator+= (T add) && { i += add; return TIndex(*this); } + constexpr TIndex operator+= (T add) && { i += add; return TIndex{*this}; } void Invalidate() { i = long(TIndex::BASE)-1; } bool IsValid() const { return i+1 != TIndex::BASE; } // operator bool() const { return IsValid(); } From 9bc9ee8e7dfe070e8bf0fee2f483d854ac0f9227 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 31 Dec 2024 13:18:03 +0100 Subject: [PATCH 07/19] remove convert operator --- libsrc/meshing/meshtype.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 08b3c87f..2bbd8659 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -201,7 +201,7 @@ namespace netgen constexpr operator T () const { return i; } explicit constexpr operator T& () { return i; } public: - constexpr operator TIndex() const { return TIndex(i); } + // constexpr operator TIndex() const { return TIndex(i); } operator TIndex&() { return static_cast(*this); } TIndex operator++ (int) { TIndex hi{*this}; i++; return hi; } @@ -250,6 +250,7 @@ namespace netgen { public: using Index::Index; + constexpr PointIndex (Index & ind) : Index(ind) { } }; } From 2b75d091e9b6978831dd436bca640bcff277a264 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 31 Dec 2024 13:31:31 +0100 Subject: [PATCH 08/19] inheriated ctor --- libsrc/meshing/meshtype.hpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 2bbd8659..5c2c740a 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -167,7 +167,7 @@ namespace netgen Index & operator= (Index&&) = default; // private: - constexpr Index (int ai) : i(ai) + constexpr Index (T ai) : i(ai) { #ifdef DEBUG if (ai < Base) @@ -446,7 +446,12 @@ namespace netgen class ElementIndex : public Index { public: - using Index::Index; + using Index::Index; + /* + constexpr ElementIndex () = default; + constexpr ElementIndex (int i) : Index(i) { } + constexpr ElementIndex (Index & ind) : Index(ind) { } + */ }; inline istream & operator>> (istream & ist, ElementIndex & pi) @@ -472,6 +477,7 @@ namespace netgen { public: using Index::Index; + constexpr SurfaceElementIndex (Index & ind) : Index(ind) { } }; From 2fdc293b9aedb26cde8234094e9c8062a56cede5 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 31 Dec 2024 13:38:04 +0100 Subject: [PATCH 09/19] ctor --- libsrc/meshing/meshtype.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 5c2c740a..da5aa339 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -220,7 +220,7 @@ namespace netgen template // constexpr auto operator+ (Index ind, int i) { return TIndex(T(ind)+i); } - constexpr auto operator+ (Index ind, int i) { return TIndex( Index(ind) +=i ); } + constexpr auto operator+ (Index ind, int i) { return TIndex{ Index(ind) +=i }; } template constexpr TIndex operator+ (Index pi, size_t i) { return TIndex(pi.i+i); } template From 104c576caa2f8a7eecf141712f5d03918adb6908 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 31 Dec 2024 13:54:11 +0100 Subject: [PATCH 10/19] ctor --- libsrc/meshing/meshtype.hpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index da5aa339..757073af 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -161,7 +161,7 @@ namespace netgen static constexpr t_invalid INVALID{}; constexpr Index () = default; - constexpr Index (const Index&) = default; + constexpr Index (const Index& i2) = default; constexpr Index (Index &&) = default; Index & operator= (const Index&) = default; Index & operator= (Index&&) = default; @@ -202,14 +202,13 @@ namespace netgen explicit constexpr operator T& () { return i; } public: // constexpr operator TIndex() const { return TIndex(i); } - operator TIndex&() { return static_cast(*this); } + // operator TIndex&() { return static_cast(*this); } TIndex operator++ (int) { TIndex hi{*this}; i++; return hi; } TIndex operator-- (int) { TIndex hi(*this); i--; return hi; } - TIndex & operator++ () { i++; return *this; } - TIndex operator-- () { i--; return *this; } - constexpr TIndex & operator+= (T add) & { i += add; return *this; } - constexpr TIndex operator+= (T add) && { i += add; return TIndex{*this}; } + TIndex operator++ () { i++; return TIndex{*this}; } + TIndex operator-- () { i--; return TIndex{*this}; } + constexpr TIndex operator+= (T add) { i += add; return TIndex{*this}; } void Invalidate() { i = long(TIndex::BASE)-1; } bool IsValid() const { return i+1 != TIndex::BASE; } // operator bool() const { return IsValid(); } @@ -220,7 +219,8 @@ namespace netgen template // constexpr auto operator+ (Index ind, int i) { return TIndex(T(ind)+i); } - constexpr auto operator+ (Index ind, int i) { return TIndex{ Index(ind) +=i }; } + // constexpr auto operator+ (Index ind, int i) { return TIndex{ Index(ind) +=i }; } + constexpr auto operator+ (Index ind, int i) { return ind += i; return TIndex(ind); } template constexpr TIndex operator+ (Index pi, size_t i) { return TIndex(pi.i+i); } template From fe21b0bb8b75015c1caf9cd23fa35b68e35aedbd Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 31 Dec 2024 14:07:27 +0100 Subject: [PATCH 11/19] cleanup --- libsrc/meshing/meshtype.hpp | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 757073af..0f2e780d 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -206,8 +206,8 @@ namespace netgen TIndex operator++ (int) { TIndex hi{*this}; i++; return hi; } TIndex operator-- (int) { TIndex hi(*this); i--; return hi; } - TIndex operator++ () { i++; return TIndex{*this}; } - TIndex operator-- () { i--; return TIndex{*this}; } + TIndex & operator++ () { i++; return static_cast(*this); } + TIndex & operator-- () { i--; return static_cast(*this); } constexpr TIndex operator+= (T add) { i += add; return TIndex{*this}; } void Invalidate() { i = long(TIndex::BASE)-1; } bool IsValid() const { return i+1 != TIndex::BASE; } @@ -218,9 +218,7 @@ namespace netgen template - // constexpr auto operator+ (Index ind, int i) { return TIndex(T(ind)+i); } - // constexpr auto operator+ (Index ind, int i) { return TIndex{ Index(ind) +=i }; } - constexpr auto operator+ (Index ind, int i) { return ind += i; return TIndex(ind); } + constexpr auto operator+ (Index ind, int i) { return TIndex(ind.i+i); } template constexpr TIndex operator+ (Index pi, size_t i) { return TIndex(pi.i+i); } template @@ -250,7 +248,6 @@ namespace netgen { public: using Index::Index; - constexpr PointIndex (Index & ind) : Index(ind) { } }; } @@ -447,11 +444,6 @@ namespace netgen { public: using Index::Index; - /* - constexpr ElementIndex () = default; - constexpr ElementIndex (int i) : Index(i) { } - constexpr ElementIndex (Index & ind) : Index(ind) { } - */ }; inline istream & operator>> (istream & ist, ElementIndex & pi) @@ -477,7 +469,6 @@ namespace netgen { public: using Index::Index; - constexpr SurfaceElementIndex (Index & ind) : Index(ind) { } }; From bcbd390f7d632d582b09c43d6c571de98f4a72a4 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 31 Dec 2024 21:26:04 +0100 Subject: [PATCH 12/19] PointIndex in Identifications --- libsrc/csg/zrefine.cpp | 23 ++++++-- libsrc/include/nginterface_v2_impl.hpp | 4 +- libsrc/interface/nginterface_v2.cpp | 10 +++- libsrc/interface/writeuser.cpp | 39 ++++++++++---- libsrc/meshing/meshclass.cpp | 34 ++++++------ libsrc/meshing/meshclass.hpp | 2 +- libsrc/meshing/meshfunc.cpp | 6 ++- libsrc/meshing/meshtype.cpp | 75 ++++++++++++++++++++------ libsrc/meshing/meshtype.hpp | 12 +++-- libsrc/meshing/python_mesh.cpp | 5 +- libsrc/meshing/secondorder.cpp | 1 + libsrc/meshing/topology.cpp | 15 ++++-- libsrc/visualization/vsmesh.cpp | 45 +++++++++------- 13 files changed, 188 insertions(+), 83 deletions(-) diff --git a/libsrc/csg/zrefine.cpp b/libsrc/csg/zrefine.cpp index 649bb50c..6e154bcb 100644 --- a/libsrc/csg/zrefine.cpp +++ b/libsrc/csg/zrefine.cpp @@ -264,14 +264,21 @@ namespace netgen { auto & identpts = mesh.GetIdentifications().GetIdentifiedPoints (); - + + /* for (int i = 1; i <= identpts.GetNBags(); i++) for (int j = 1; j <= identpts.GetBagSize(i); j++) { INDEX_3 pair; int dummy; identpts.GetData(i, j, pair, dummy); - auto idnr = pair[2]; + */ + for (auto [hash, val] : identpts)\ + { + auto [hash_pts, idnr] = hash; + auto [pi1, pi2] = hash_pts; + // auto idnr = pair[2]; + const CloseSurfaceIdentification * csid = dynamic_cast (geom->identifications.Get(idnr)); @@ -282,17 +289,25 @@ namespace netgen if (first_id.Test (idnr)) { first_id.Clear(idnr); + /* ref_uniform.Append (INDEX_3 (pair.I1(), pair.I2(), csid->RefLevels())); ref_singular.Append (INDEX_3 (pair.I1(), pair.I2(), csid->RefLevels1())); ref_singular.Append (INDEX_3 (pair.I2(), pair.I1(), csid->RefLevels2())); + */ + ref_uniform.Append (INDEX_3 (pi1, pi2, csid->RefLevels())); + ref_singular.Append (INDEX_3 (pi1, pi2, csid->RefLevels1())); + ref_singular.Append (INDEX_3 (pi2, pi1, csid->RefLevels2())); + } } else { //const NgArray & slices = csid->GetSlices(); INDEX_4 i4; - i4[0] = pair.I1(); - i4[1] = pair.I2(); + // i4[0] = pair.I1(); + // i4[1] = pair.I2(); + i4[0] = pi1; + i4[1] = pi2; i4[2] = idnr; i4[3] = csid->GetSlices().Size(); ref_slices.Append (i4); diff --git a/libsrc/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp index 19c41257..d0db9a2b 100644 --- a/libsrc/include/nginterface_v2_impl.hpp +++ b/libsrc/include/nginterface_v2_impl.hpp @@ -339,8 +339,8 @@ NGX_INLINE DLL_HEADER Ng_Buffer Ngx_Mesh :: GetPeriodicVertices(int idnr mesh->GetIdentifications().GetPairs (idnr+1, apairs); for(auto& ind : apairs) { - ind.I1()--; - ind.I2()--; + ind.I1() -= IndexBASE(); + ind.I2() -= IndexBASE(); } typedef int ti2[2]; return { apairs.Size(), (ti2*)(void*)apairs.Release() }; diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp index 9233c5a7..1e3e7af8 100644 --- a/libsrc/interface/nginterface_v2.cpp +++ b/libsrc/interface/nginterface_v2.cpp @@ -1249,9 +1249,15 @@ int Ngx_Mesh::GetElementOrder (int enr) const void Ngx_Mesh::GetElementOrders (int enr, int * ox, int * oy, int * oz) const { if (mesh->GetDimension() == 3) - mesh->VolumeElement(enr).GetOrder(*ox, *oy, *oz); + { + ElementIndex ei = IndexBASE() + enr-1; + mesh->VolumeElement(ei).GetOrder(*ox, *oy, *oz); + } else - mesh->SurfaceElement(enr).GetOrder(*ox, *oy, *oz); + { + SurfaceElementIndex sei = IndexBASE() + enr-1; + mesh->SurfaceElement(sei).GetOrder(*ox, *oy, *oz); + } } void Ngx_Mesh::SetElementOrder (int enr, int order) diff --git a/libsrc/interface/writeuser.cpp b/libsrc/interface/writeuser.cpp index 27370eb0..b4f5f033 100644 --- a/libsrc/interface/writeuser.cpp +++ b/libsrc/interface/writeuser.cpp @@ -94,9 +94,13 @@ void WriteNeutralFormat (const Mesh & mesh, if (mesh.GetDimension() == 3) { outfile << ne << "\n"; + /* for (int i = 1; i <= ne; i++) { Element el = mesh.VolumeElement(i); + */ + for (Element el : mesh.VolumeElements()) + { if (inverttets) el.Invert(); outfile.width(4); @@ -112,9 +116,13 @@ void WriteNeutralFormat (const Mesh & mesh, } outfile << nse << "\n"; + /* for (int i = 1; i <= nse; i++) { Element2d el = mesh.SurfaceElement(i); + */ + for (Element2d el : mesh.SurfaceElements()) + { if (invertsurf) el.Invert(); outfile.width(4); @@ -583,9 +591,13 @@ void WriteFEPPFormat (const Mesh & mesh, outfile << ne << "\n"; + /* for (i = 1; i <= ne; i++) { const Element & el = mesh.VolumeElement(i); + */ + for (const Element & el : mesh.VolumeElements()) + { outfile.width(4); outfile << el.GetIndex() << " "; outfile.width(4); @@ -677,7 +689,6 @@ void WriteEdgeElementFormat (const Mesh & mesh, int nelements = mesh.GetNE(); int nsurfelem = mesh.GetNSE(); int nedges = top->GetNEdges(); - int i, j; int inverttets = mparam.inverttets; int invertsurf = mparam.inverttrigs; @@ -692,7 +703,7 @@ void WriteEdgeElementFormat (const Mesh & mesh, // vertices with coordinates outfile << npoints << "\n"; - for (i = 1; i <= npoints; i++) + for (int i = 1; i <= npoints; i++) { const Point3d & p = mesh.Point(i); @@ -706,16 +717,24 @@ void WriteEdgeElementFormat (const Mesh & mesh, // element - edge - list outfile << nelements << " " << nedges << "\n"; + /* for (i = 1; i <= nelements; i++) { Element el = mesh.VolumeElement(i); + */ + for (ElementIndex ei : Range(mesh.VolumeElements())) + { + int i = ei-IndexBASE(ei)+1; + + Element el = mesh.VolumeElement(ei); + if (inverttets) el.Invert(); outfile.width(4); outfile << el.GetIndex() << " "; outfile.width(8); outfile << el.GetNP(); - for (j = 1; j <= el.GetNP(); j++) + for (int j = 1; j <= el.GetNP(); j++) { outfile << " "; outfile.width(8); @@ -723,11 +742,11 @@ void WriteEdgeElementFormat (const Mesh & mesh, } // top->GetElementEdges(i,edges); - auto eledges = top->GetEdges(ElementIndex(i-1)); + auto eledges = top->GetEdges(ei); outfile << endl << " "; outfile.width(8); outfile << eledges.Size(); - for (j=1; j <= eledges.Size(); j++) + for (int j=1; j <= eledges.Size(); j++) { outfile << " "; outfile.width(8); @@ -738,7 +757,7 @@ void WriteEdgeElementFormat (const Mesh & mesh, // orientation: top->GetElementEdgeOrientations(i,edges); outfile << " "; - for (j=1; j <= edges.Size(); j++) + for (int j=1; j <= edges.Size(); j++) { outfile << " "; outfile.width(8); @@ -749,7 +768,7 @@ void WriteEdgeElementFormat (const Mesh & mesh, // surface element - edge - list (with boundary conditions) outfile << nsurfelem << "\n"; - for (i = 1; i <= nsurfelem; i++) + for (int i = 1; i <= nsurfelem; i++) { Element2d el = mesh.SurfaceElement(i); if (invertsurf) @@ -758,7 +777,7 @@ void WriteEdgeElementFormat (const Mesh & mesh, outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " "; outfile.width(8); outfile << el.GetNP(); - for (j = 1; j <= el.GetNP(); j++) + for (int j = 1; j <= el.GetNP(); j++) { outfile << " "; outfile.width(8); @@ -769,7 +788,7 @@ void WriteEdgeElementFormat (const Mesh & mesh, outfile << endl << " "; outfile.width(8); outfile << edges.Size(); - for (j=1; j <= edges.Size(); j++) + for (int j=1; j <= edges.Size(); j++) { outfile << " "; outfile.width(8); @@ -782,7 +801,7 @@ void WriteEdgeElementFormat (const Mesh & mesh, // int v1, v2; // edge - vertex - list outfile << nedges << "\n"; - for (i=1; i <= nedges; i++) + for (int i=1; i <= nedges; i++) { // top->GetEdgeVertices(i,v1,v2); auto [v1,v2] = top->GetEdgeVertices(i-1); diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 16ee4a0a..4dae9b84 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -6438,22 +6438,24 @@ namespace netgen mapped_points = false; // Add new points - for(auto [p1p2, dummy] : identpts) + for(auto [hash, dummy] : identpts) { - if(p1p2[2] != nr) + auto [hash_pts, hash_nr] = hash; + if(hash_nr != nr) continue; - auto& ipts = inserted_points[{p1p2.I1(), p1p2.I2()}]; - auto p1 = Point(p1p2.I1()); - auto p2 = Point(p1p2.I2()); - ipts.Append(p1p2.I1()); - mapped_points.SetBit(p1p2.I1()); + // auto& ipts = inserted_points[{p1p2.I1(), p1p2.I2()}]; + auto& ipts = inserted_points[ { hash_pts[0], hash_pts[1] }]; + auto p1 = Point(hash_pts.I1()); + auto p2 = Point(hash_pts.I2()); + ipts.Append(hash_pts.I1()); + mapped_points.SetBit(hash_pts.I1()); for(auto slice : slices) { auto np = p1 + slice * (p2-p1); auto npi = AddPoint(np); ipts.Append(npi); } - ipts.Append(p1p2.I2()); + ipts.Append(hash_pts.I2()); } // Split segments @@ -6644,21 +6646,21 @@ namespace netgen void Mesh :: CalcMinMaxAngle (double badellimit, double * retvalues) { - int i, j; int lpi1, lpi2, lpi3, lpi4; double phimax = 0, phimin = 10; double facephimax = 0, facephimin = 10; int illegaltets = 0, negativetets = 0, badtets = 0; - for (i = 1; i <= GetNE(); i++) + // for (int i = 1; i <= GetNE(); i++) + for (ElementIndex ei : Range(VolumeElements())) { int badel = 0; - Element & el = VolumeElement(i); + Element & el = VolumeElement(ei); if (el.GetType() != TET) { - VolumeElement(i).Flags().badel = 0; + VolumeElement(ei).Flags().badel = 0; continue; } @@ -6673,8 +6675,8 @@ namespace netgen { badel = 1; illegaltets++; - (*testout) << "illegal tet: " << i << " "; - for (j = 1; j <= el.GetNP(); j++) + (*testout) << "illegal tet: " << ei << " "; + for (int j = 1; j <= el.GetNP(); j++) (*testout) << el.PNum(j) << " "; (*testout) << endl; } @@ -6713,7 +6715,7 @@ namespace netgen // angles in faces - for (j = 1; j <= 4; j++) + for (int j = 1; j <= 4; j++) { Element2d face(TRIG); el.GetFace (j, face); @@ -6740,7 +6742,7 @@ namespace netgen } - VolumeElement(i).Flags().badel = badel; + VolumeElement(ei).Flags().badel = badel; if (badel) badtets++; } diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index 90619969..f89dddc1 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -377,7 +377,7 @@ namespace netgen const auto & SurfaceElements() const { return surfelements; } auto & SurfaceElements() { return surfelements; } - + DLL_HEADER void RebuildSurfaceElementLists (); DLL_HEADER void GetSurfaceElementsOfFace (int facenr, Array & sei) const; diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 44040fe3..5716f928 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -578,8 +578,10 @@ namespace netgen el.SetIndex(m_.domain); mesh.AddVolumeElement(el); } - for(const auto& [p1p2, dummy] : m.GetIdentifications().GetIdentifiedPoints()) - mesh.GetIdentifications().Add(pmap[p1p2[0]], pmap[p1p2[1]], p1p2[2]); + // for(const auto& [p1p2, dummy] : m.GetIdentifications().GetIdentifiedPoints()) + // mesh.GetIdentifications().Add(pmap[p1p2[0]], pmap[p1p2[1]], p1p2[2]); + for(const auto& [p1p2, dummy] : m.GetIdentifications().GetIdentifiedPoints()) + mesh.GetIdentifications().Add( pmap[ get<0>(p1p2)[0] ], pmap[ get<0>(p1p2)[1]] , get<1>(p1p2) ); for(auto i : Range(m.GetIdentifications().GetMaxNr())) { mesh.GetIdentifications().SetType(i+1, m.GetIdentifications().GetType(i+1)); diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index 24e9f398..2e954f54 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -2708,8 +2708,9 @@ namespace netgen void Identifications :: DoArchive (Archive & ar) { ar & maxidentnr; - ar & identifiedpoints & identifiedpoints_nr; - + // ar & identifiedpoints & identifiedpoints_nr; +#pragma message( "Archive CloseHadhTable missing " __FILE__ ) + ar & idpoints_table; if (ar.Output()) { @@ -2735,11 +2736,11 @@ namespace netgen void Identifications :: Add (PointIndex pi1, PointIndex pi2, int identnr) { // (*testout) << "Identification::Add, pi1 = " << pi1 << ", pi2 = " << pi2 << ", identnr = " << identnr << endl; - INDEX_2 pair (pi1, pi2); + PointIndices<2> pair (pi1, pi2); identifiedpoints.Set (pair, identnr); - INDEX_3 tripl (pi1, pi2, identnr); - identifiedpoints_nr.Set (tripl, 1); + // INDEX_3 tripl (pi1, pi2, identnr); + identifiedpoints_nr.Set ( { { pi1, pi2 }, identnr }, 1); if (identnr > maxidentnr) maxidentnr = identnr; names.SetSize(maxidentnr); @@ -2762,8 +2763,9 @@ namespace netgen bool Identifications :: Get (PointIndex pi1, PointIndex pi2, int nr) const { - INDEX_3 tripl(pi1, pi2, nr); - if (identifiedpoints_nr.Used (tripl)) + // INDEX_3 tripl(pi1, pi2, nr); + // if (identifiedpoints_nr.Used (tripl)) + if (identifiedpoints_nr.Used ( { { pi1, pi1 }, nr } ) ) return 1; else return 0; @@ -2803,23 +2805,30 @@ namespace netgen { cout << "getmap, identnr = " << identnr << endl; + /* for (int i = 1; i <= identifiedpoints_nr.GetNBags(); i++) for (int j = 1; j <= identifiedpoints_nr.GetBagSize(i); j++) + */ + for (auto [hash, val] : identifiedpoints_nr) { + /* INDEX_3 i3; int dummy; identifiedpoints_nr.GetData (i, j, i3, dummy); - - if (i3.I3() == identnr || !identnr) + */ + + auto [hash_pts, hash_nr] = hash; + + if (hash_nr == identnr || !identnr) { /* identmap.Elem(i3.I1()) = i3.I2(); if(symmetric) identmap.Elem(i3.I2()) = i3.I1(); */ - identmap[i3.I1()] = i3.I2(); + identmap[hash_pts.I1()] = hash_pts.I2(); if(symmetric) - identmap[i3.I2()] = i3.I1(); + identmap[hash_pts.I2()] = hash_pts.I1(); } } } @@ -2830,8 +2839,12 @@ namespace netgen Array Identifications :: GetPairs () const { Array pairs; - for(auto [i3, dummy] : identifiedpoints_nr) - pairs.Append(i3); + for(auto [hash, dummy] : identifiedpoints_nr) + // pairs.Append(i3); + { + auto [pts,nr] = hash; + pairs.Append ( { pts[0], pts[1], nr } ); + } return pairs; } @@ -2841,6 +2854,8 @@ namespace netgen identpairs.SetSize(0); if (identnr == 0) + { + /* for (int i = 1; i <= identifiedpoints.GetNBags(); i++) for (int j = 1; j <= identifiedpoints.GetBagSize(i); j++) { @@ -2848,8 +2863,14 @@ namespace netgen int nr; identifiedpoints.GetData (i, j, i2, nr); identpairs.Append (i2); - } + } + */ + for (auto [hash,val] : identifiedpoints) + identpairs.Append (hash); + } else + { + /* for (int i = 1; i <= identifiedpoints_nr.GetNBags(); i++) for (int j = 1; j <= identifiedpoints_nr.GetBagSize(i); j++) { @@ -2859,12 +2880,21 @@ namespace netgen if (i3.I3() == identnr) identpairs.Append (INDEX_2(i3.I1(), i3.I2())); - } + } + */ + for (auto [hash,val] : identifiedpoints_nr) + { + auto [hash_pts, hash_nr] = hash; + if (hash_nr == identnr) + identpairs.Append (hash_pts); + } + } } void Identifications :: SetMaxPointNr (int maxpnum) { + /* for (int i = 1; i <= identifiedpoints.GetNBags(); i++) for (int j = 1; j <= identifiedpoints.GetBagSize(i); j++) { @@ -2878,6 +2908,18 @@ namespace netgen identifiedpoints.SetData (i, j, i2, -1); } } + */ + + // can we get data by reference ? + for (auto [hash,data] : identifiedpoints) + { + if (hash.I1() > IndexBASE()+maxpnum-1 || + hash.I2() > IndexBASE()+maxpnum-1) + { + identifiedpoints[hash] = -1; + } + + } } // Map points in the identifications to new point numbers @@ -2900,7 +2942,8 @@ namespace netgen { ost << "Identifications:" << endl; ost << "pairs: " << endl << identifiedpoints << endl; - ost << "pairs and nr: " << endl << identifiedpoints_nr << endl; + // ost << "pairs and nr: " << endl << identifiedpoints_nr << endl; +#pragma message( "Can't ostream a tuple " __FILE__ ) ost << "table: " << endl << idpoints_table << endl; } diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 0f2e780d..91c5123a 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -244,7 +244,7 @@ namespace netgen - class PointIndex : public Index + class PointIndex : public Index { public: using Index::Index; @@ -1763,13 +1763,15 @@ namespace netgen class Mesh & mesh; /// identify points (thin layers, periodic b.c.) - INDEX_2_HASHTABLE identifiedpoints; + // INDEX_2_HASHTABLE identifiedpoints; + ClosedHashTable, int> identifiedpoints; /// the same, with info about the id-nr - INDEX_3_HASHTABLE identifiedpoints_nr; + // INDEX_3_HASHTABLE identifiedpoints_nr; + ClosedHashTable, int>, int> identifiedpoints_nr; /// sorted by identification nr - TABLE idpoints_table; + TABLE> idpoints_table; NgArray type; @@ -1805,7 +1807,7 @@ namespace netgen // bool HasIdentifiedPoints() const { return identifiedpoints != nullptr; } /// - INDEX_3_HASHTABLE & GetIdentifiedPoints () + auto & GetIdentifiedPoints () { return identifiedpoints_nr; } diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 1a133165..c9c95002 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -1266,7 +1266,10 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) py::list points; for(const auto& pair : self.GetIdentifications().GetIdentifiedPoints()) { - py::tuple pnts = py::make_tuple(pair.first.I1(), pair.first.I2()); + // py::tuple pnts = py::make_tuple(pair.first.I1(), pair.first.I2()); + + auto [pi1, pi2] = get<0> (pair.first); + py::tuple pnts = py::make_tuple(pi1, pi2); points.append(pnts); } return points; diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp index 5527b2a8..56dc7583 100644 --- a/libsrc/meshing/secondorder.cpp +++ b/libsrc/meshing/secondorder.cpp @@ -558,6 +558,7 @@ namespace netgen // (*testout) << "bad els: " << endl; wrongels = 0; for (int i = 1; i <= ne; i++) + { if (!illegalels.Test(i) && mesh.VolumeElement(i). diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp index ce4f9ce7..2af9c76b 100644 --- a/libsrc/meshing/topology.cpp +++ b/libsrc/meshing/topology.cpp @@ -1829,8 +1829,8 @@ namespace netgen void MeshTopology :: GetElementFaces (int elnr, NgArray & elfaces, bool withorientation) const { - int nfa = GetNFaces (mesh->VolumeElement(elnr).GetType()); ElementIndex ei = IndexBASE() +(elnr-1); + int nfa = GetNFaces (mesh->VolumeElement(ei).GetType()); elfaces.SetSize (nfa); @@ -1859,7 +1859,8 @@ namespace netgen void MeshTopology :: GetElementEdgeOrientations (int elnr, NgArray & eorient) const { - int ned = GetNEdges (mesh->VolumeElement(elnr).GetType()); + ElementIndex ei = IndexBASE() +(elnr-1); + int ned = GetNEdges (mesh->VolumeElement(ei).GetType()); eorient.SetSize (ned); for (int i = 1; i <= ned; i++) // eorient.Elem(i) = (edges.Get(elnr)[i-1] > 0) ? 1 : -1; @@ -1869,7 +1870,8 @@ namespace netgen void MeshTopology :: GetElementFaceOrientations (int elnr, NgArray & forient) const { - int nfa = GetNFaces (mesh->VolumeElement(elnr).GetType()); + ElementIndex ei = IndexBASE() +(elnr-1); + int nfa = GetNFaces (mesh->VolumeElement(ei).GetType()); forient.SetSize (nfa); for (int i = 1; i <= nfa; i++) // forient.Elem(i) = faces.Get(elnr)[i-1].forient; @@ -2111,7 +2113,9 @@ namespace netgen int MeshTopology :: GetElementEdgeOrientation (int elnr, int locedgenr) const { - const Element & el = mesh->VolumeElement (elnr); + ElementIndex ei = IndexBASE() +(elnr-1); + + const Element & el = mesh->VolumeElement (ei); const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType()); int k = locedgenr; @@ -2123,7 +2127,8 @@ namespace netgen int MeshTopology :: GetElementFaceOrientation (int elnr, int locfacenr) const { - const Element & el = mesh->VolumeElement (elnr); + ElementIndex ei = IndexBASE() +(elnr-1); + const Element & el = mesh->VolumeElement (ei); const ELEMENT_FACE * elfaces = MeshTopology::GetFaces0 (el.GetType()); diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp index 0e45aa12..f8f6c5d9 100644 --- a/libsrc/visualization/vsmesh.cpp +++ b/libsrc/visualization/vsmesh.cpp @@ -515,20 +515,21 @@ namespace netgen if (vispar.drawelementnumbers) { NgArray v; - for (int i = 1; i <= mesh->GetNE(); i++) + // for (int i = 1; i <= mesh->GetNE(); i++) + for (ElementIndex ei : Range(mesh->VolumeElements())) { // const ELEMENTTYPE & eltype = mesh->ElementType(i); NgArray pnums; Point3d p; - const Element & el = mesh->VolumeElement (i); + const Element & el = mesh->VolumeElement (ei); if ( ! el.PNum(5)) // eltype == TET ) { pnums.SetSize(4); for( int j = 0; j < pnums.Size(); j++) - pnums[j] = mesh->VolumeElement(i).PNum(j+1); + pnums[j] = mesh->VolumeElement(ei).PNum(j+1); const Point3d & p1 = mesh->Point(pnums[0]); @@ -541,7 +542,7 @@ namespace netgen { pnums.SetSize(5); for( int j = 0; j < pnums.Size(); j++) - pnums[j] = mesh->VolumeElement(i).PNum(j+1); + pnums[j] = mesh->VolumeElement(ei).PNum(j+1); const Point3d & p1 = mesh->Point(pnums[0]); @@ -559,7 +560,7 @@ namespace netgen { pnums.SetSize(6); for( int j = 0; j < pnums.Size(); j++) - pnums[j] = mesh->VolumeElement(i).PNum(j+1); + pnums[j] = mesh->VolumeElement(ei).PNum(j+1); const Point3d & p1 = mesh->Point(pnums[0]); const Point3d & p2 = mesh->Point(pnums[1]); @@ -574,7 +575,7 @@ namespace netgen { pnums.SetSize(8); for( int j = 0; j < pnums.Size(); j++) - pnums[j] = mesh->VolumeElement(i).PNum(j+1); + pnums[j] = mesh->VolumeElement(ei).PNum(j+1); const Point3d & p1 = mesh->Point(pnums[0]); const Point3d & p2 = mesh->Point(pnums[1]); @@ -589,7 +590,7 @@ namespace netgen } glRasterPos3d (p.X(), p.Y(), p.Z()); - snprintf (buf, size(buf), "%d", i); + snprintf (buf, size(buf), "%d", ei-IndexBASE(ei)); // glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf); MyOpenGLText (buf); @@ -617,14 +618,15 @@ namespace netgen static float badelcol[] = { 1.0f, 0.0f, 1.0f, 1.0f }; glLineWidth (1.0f); - for (int i = 1; i <= mesh->GetNE(); i++) + //for (int i = 1; i <= mesh->GetNE(); i++) + for (ElementIndex ei : Range(mesh->VolumeElements())) { - if (mesh->VolumeElement(i).Flags().badel || - mesh->VolumeElement(i).Flags().illegal || - (i == vispar.drawelement)) + if (mesh->VolumeElement(ei).Flags().badel || + mesh->VolumeElement(ei).Flags().illegal || + (ei-IndexBASE(ei) == vispar.drawelement)) { // copy to be thread-safe - Element el = mesh->VolumeElement (i); + Element el = mesh->VolumeElement (ei); el.GetSurfaceTriangles (faces); glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, badelcol); @@ -674,9 +676,9 @@ namespace netgen } - for (int i = 1; i <= mesh->GetNE(); i++) + for (ElementIndex ei : Range(mesh->VolumeElements())) { - Element el = mesh->VolumeElement (i); + Element el = mesh->VolumeElement (ei); int hascp = 0; for (int j = 1; j <= el.GetNP(); j++) if (el.PNum(j) == vispar.centerpoint) @@ -684,7 +686,7 @@ namespace netgen if (hascp) { - (*testout) << "draw el " << i << " : "; + (*testout) << "draw el " << ei << " : "; for (int j = 1; j <= el.GetNP(); j++) (*testout) << el.PNum(j) << " "; (*testout) << endl; @@ -860,17 +862,22 @@ namespace netgen { auto & idpts = mesh->GetIdentifications().GetIdentifiedPoints(); - + + /* for (int i = 1; i <= idpts.GetNBags(); i++) for (int j = 1; j <= idpts.GetBagSize(i); j++) { INDEX_3 pts; int dummy; // , val; - idpts.GetData (i, j, pts, dummy); + */ + for (auto [hash, val] : idpts) + { + auto [hash_pts, hash_nr] = hash; + auto [pi1, pi2] = hash_pts; // val = pts[2]; - const Point3d & p1 = mesh->Point(pts.I1()); - const Point3d & p2 = mesh->Point(pts.I2()); + const Point3d & p1 = mesh->Point(pi1); + const Point3d & p2 = mesh->Point(pi2); glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, identifiedcol); From 6d6e297a1f3a221fe977f538eb0029c41eee5dce Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Wed, 1 Jan 2025 12:09:03 +0100 Subject: [PATCH 13/19] PointInd for edge on closed surf --- libsrc/csg/edgeflw.cpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/libsrc/csg/edgeflw.cpp b/libsrc/csg/edgeflw.cpp index fe57471f..6b742ff0 100644 --- a/libsrc/csg/edgeflw.cpp +++ b/libsrc/csg/edgeflw.cpp @@ -1885,12 +1885,10 @@ namespace netgen if (seg1.domin != -1 || seg1.domout != -1) { - mesh.AddPoint (p1, layer, EDGEPOINT); - mesh.AddPoint (p2, layer, EDGEPOINT); - seg1[0] = mesh.GetNP()-1; - seg1[1] = mesh.GetNP(); - seg2[1] = mesh.GetNP()-1; - seg2[0] = mesh.GetNP(); + seg1[0] = mesh.AddPoint (p1, layer, EDGEPOINT); + seg1[1] = mesh.AddPoint (p2, layer, EDGEPOINT); + seg2[0] = seg1[1]; + seg2[1] = seg1[0]; seg1.geominfo[0].trignum = 1; seg1.geominfo[1].trignum = 1; seg2.geominfo[0].trignum = 1; From 3185256ad39e7a04716adb1da76f2ad9cb48a7f3 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Wed, 1 Jan 2025 12:27:44 +0100 Subject: [PATCH 14/19] PointIndex for csg lockedpnts --- libsrc/csg/genmesh.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libsrc/csg/genmesh.cpp b/libsrc/csg/genmesh.cpp index 4aa292cf..c9bba2f0 100644 --- a/libsrc/csg/genmesh.cpp +++ b/libsrc/csg/genmesh.cpp @@ -35,7 +35,7 @@ namespace netgen auto up = geom.GetUserPoint(i); auto pnum = mesh.AddPoint(up); mesh.Points().Last().Singularity (geom.GetUserPointRefFactor(i)); - mesh.AddLockedPoint (PointIndex (i+1)); + mesh.AddLockedPoint (pnum); int index = up.GetIndex(); if (index == -1) index = mesh.AddCD3Name (up.GetName())+1; From aca27f64214e4e9596948710ab4f5681ff7689a7 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Wed, 1 Jan 2025 13:53:01 +0100 Subject: [PATCH 15/19] some more int->PointIndex --- libsrc/csg/edgeflw.cpp | 8 ++++---- libsrc/csg/identify.cpp | 30 +++++++++++++++++++++--------- libsrc/meshing/secondorder.cpp | 4 ++-- 3 files changed, 27 insertions(+), 15 deletions(-) diff --git a/libsrc/csg/edgeflw.cpp b/libsrc/csg/edgeflw.cpp index 6b742ff0..afa8184f 100644 --- a/libsrc/csg/edgeflw.cpp +++ b/libsrc/csg/edgeflw.cpp @@ -1714,12 +1714,12 @@ namespace netgen if (oldseg.seginfo == 0) continue; - int pi1 = oldseg[0]; - int pi2 = oldseg[1]; + PointIndex pi1 = oldseg[0]; + PointIndex pi2 = oldseg[1]; - int npi1 = geometry.identifications.Get(copyedgeidentification) + PointIndex npi1 = geometry.identifications.Get(copyedgeidentification) -> GetIdentifiedPoint (mesh, pi1); - int npi2 = geometry.identifications.Get(copyedgeidentification) + PointIndex npi2 = geometry.identifications.Get(copyedgeidentification) -> GetIdentifiedPoint (mesh, pi2); //(*testout) << "copy edge, pts = " << npi1 << " - " << npi2 << endl; diff --git a/libsrc/csg/identify.cpp b/libsrc/csg/identify.cpp index 8878c695..e87a6679 100644 --- a/libsrc/csg/identify.cpp +++ b/libsrc/csg/identify.cpp @@ -289,14 +289,14 @@ GetIdentifiedPoint (class Mesh & mesh, PointIndex pi) // project to other surface snew->Project (hp); - int newpi = 0; - for (int i = 1; i <= mesh.GetNP(); i++) - if (Dist2 (mesh.Point(i), hp) < 1e-12) + PointIndex newpi(PointIndex::INVALID); + for (PointIndex pi : Range(mesh.Points())) + if (Dist2 (mesh.Point(pi), hp) < 1e-12) { - newpi = i; + newpi = pi; break; } - if (!newpi) + if (!newpi.IsValid()) newpi = mesh.AddPoint (hp); if (snew == s2) @@ -322,6 +322,7 @@ void PeriodicIdentification :: IdentifyPoints (class Mesh & mesh) mesh.GetBox(p1, p2); auto eps = 1e-6 * (p2-p1).Length(); + /* for (int i = 1; i <= mesh.GetNP(); i++) { Point<3> p = mesh.Point(i); @@ -334,13 +335,24 @@ void PeriodicIdentification :: IdentifyPoints (class Mesh & mesh) if (Dist2(mesh.Point(j), pp) < eps) { mesh.GetIdentifications().Add (i, j, nr); - /* - (*testout) << "Identify points(periodic:), nr = " << nr << ": " - << mesh.Point(i) << " - " << mesh.Point(j) << endl; - */ } } } + */ + + for (auto pi : Range(mesh.Points())) + { + Point<3> p = mesh[pi]; + if (s1->PointOnSurface (p)) + { + Point<3> pp = p; + pp = trafo(pp); + s2->Project (pp); + for (PointIndex pj : Range(mesh.Points())) + if (Dist2(mesh[pj], pp) < eps) + mesh.GetIdentifications().Add (pi, pj, nr); + } + } mesh.GetIdentifications().SetType(nr,Identifications::PERIODIC); } diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp index 56dc7583..a29c1271 100644 --- a/libsrc/meshing/secondorder.cpp +++ b/libsrc/meshing/secondorder.cpp @@ -174,8 +174,8 @@ namespace netgen int nnp = newel.GetNP(); for (int j = 0; j < nnp-onp; j++) { - int pi1 = newel[betw[j][0]]; - int pi2 = newel[betw[j][1]]; + PointIndex pi1 = newel[betw[j][0]]; + PointIndex pi2 = newel[betw[j][1]]; INDEX_2 i2 = INDEX_2::Sort (pi1, pi2); From ce5f6d695ce6b3187002178edeb7ad7ac400abba Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Wed, 1 Jan 2025 15:53:23 +0100 Subject: [PATCH 16/19] all tests passing for PointIndex::BASE=0 --- libsrc/include/nginterface_v2_impl.hpp | 2 +- libsrc/meshing/meshtype.hpp | 2 +- libsrc/meshing/topology.cpp | 12 ++++++------ 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/libsrc/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp index d0db9a2b..305c0073 100644 --- a/libsrc/include/nginterface_v2_impl.hpp +++ b/libsrc/include/nginterface_v2_impl.hpp @@ -327,7 +327,7 @@ template <> NGX_INLINE DLL_HEADER const Ng_Node<2> Ngx_Mesh :: GetNode<2> (int n { Ng_Node<2> node; node.vertices.ptr = (const int*)mesh->GetTopology().GetFaceVerticesPtr(nr); - node.vertices.nv = (node.vertices.ptr[3] == 0) ? 3 : 4; + node.vertices.nv = (node.vertices.ptr[3]+1 == PointIndex::BASE) ? 3 : 4; node.surface_el = mesh->GetTopology().GetFace2SurfaceElement (nr+1)-1; return node; } diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 91c5123a..4f988b81 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -244,7 +244,7 @@ namespace netgen - class PointIndex : public Index + class PointIndex : public Index { public: using Index::Index; diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp index 2af9c76b..622512fb 100644 --- a/libsrc/meshing/topology.cpp +++ b/libsrc/meshing/topology.cpp @@ -160,7 +160,7 @@ namespace netgen { // triangle PointIndices<4> face(el[elfaces[j][0]], el[elfaces[j][1]], - el[elfaces[j][2]], 0); + el[elfaces[j][2]], PointIndex(PointIndex::INVALID)); [[maybe_unused]] int facedir = 0; @@ -266,7 +266,7 @@ namespace netgen PointIndices<4> face(el.PNum(elfaces[0][0]), el.PNum(elfaces[0][1]), - el.PNum(elfaces[0][2]),0); + el.PNum(elfaces[0][2]), PointIndex(PointIndex::INVALID)); // facedir = 0; if (face[0] > face[1]) @@ -1161,7 +1161,7 @@ namespace netgen size_t pos; if (vert2face.PositionCreate(face, pos)) { - face2vert[nfa] = { face[0], face[1], face[2], 0 }; // i4; + face2vert[nfa] = { face[0], face[1], face[2], PointIndex::BASE-1 }; // i4; vert2face.SetData (pos, face, nfa); nfa++; } @@ -2136,7 +2136,7 @@ namespace netgen if (elfaces[j][3] < 0) { // triangle INDEX_4 face(el[elfaces[j][0]], el[elfaces[j][1]], - el[elfaces[j][2]], 0); + el[elfaces[j][2]], PointIndex::BASE-1 ); int facedir = 0; if (face.I1() > face.I2()) @@ -2203,7 +2203,7 @@ namespace netgen if (elfaces[j][3] < 0) { // triangle INDEX_4 face(el[elfaces[j][0]], el[elfaces[j][1]], - el[elfaces[j][2]], 0); + el[elfaces[j][2]], PointIndex(PointIndex::INVALID)); int facedir = 0; if (face.I1() > face.I2()) @@ -2272,7 +2272,7 @@ namespace netgen vertices.SetSize(4); for (int i = 0; i < 4; i++) vertices[i] = face2vert[fnr-1][i]; - if (vertices[3] == 0) + if (vertices[3]+1==PointIndex::BASE) vertices.SetSize(3); } From e926071bb2f24de87858394a119ddeb6a59a6534 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Wed, 1 Jan 2025 16:42:11 +0100 Subject: [PATCH 17/19] archiving ngscore::CloseHashTable --- libsrc/core/hashtable.hpp | 6 ++++++ libsrc/meshing/meshtype.cpp | 5 ++--- libsrc/meshing/meshtype.hpp | 10 ++++++++++ libsrc/meshing/secondorder.cpp | 1 - 4 files changed, 18 insertions(+), 4 deletions(-) diff --git a/libsrc/core/hashtable.hpp b/libsrc/core/hashtable.hpp index af27155c..555838a0 100644 --- a/libsrc/core/hashtable.hpp +++ b/libsrc/core/hashtable.hpp @@ -834,6 +834,12 @@ namespace ngcore hash = T_HASH(invalid); used = 0; } + + void DoArchive (Archive & ar) + { + ar & hash & cont; + ar & size & mask & used; + } class Iterator { diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index 2e954f54..839c2bdb 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -2708,10 +2708,9 @@ namespace netgen void Identifications :: DoArchive (Archive & ar) { ar & maxidentnr; - // ar & identifiedpoints & identifiedpoints_nr; -#pragma message( "Archive CloseHadhTable missing " __FILE__ ) - + ar & identifiedpoints & identifiedpoints_nr; ar & idpoints_table; + if (ar.Output()) { size_t s = type.Size(); diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 4f988b81..6f70248f 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -292,6 +292,9 @@ namespace netgen constexpr PointIndex operator[] (int i) const { return PointIndex(INDEX_2::operator[](i)); } PointIndex & operator[] (int i) { return reinterpret_cast(INDEX_2::operator[](i)); } + template + void DoArchive(ARCHIVE& ar) { ar.Do(&I1(), 2); } + PointIndex & I1 () { return (*this)[0]; } PointIndex & I2 () { return (*this)[1]; } PointIndex I1 () const { return (*this)[0]; } @@ -314,6 +317,10 @@ namespace netgen constexpr PointIndices (PointIndex i1, PointIndex i2, PointIndex i3) : INDEX_3(i1,i2,i3) { ; } PointIndex operator[] (int i) const { return PointIndex(INDEX_3::operator[](i)); } PointIndex & operator[] (int i) { return reinterpret_cast(INDEX_3::operator[](i)); } + + template + void DoArchive(ARCHIVE& ar) { ar.Do(&I1(), 3); } + PointIndex & I1 () { return (*this)[0]; } PointIndex & I2 () { return (*this)[1]; } PointIndex & I3 () { return (*this)[2]; } @@ -336,6 +343,9 @@ namespace netgen PointIndex operator[] (int i) const { return PointIndex(INDEX_4::operator[](i)); } PointIndex & operator[] (int i) { return reinterpret_cast(INDEX_4::operator[](i)); } + template + void DoArchive(ARCHIVE& ar) { ar.Do(&I1(), 4); } + PointIndex & I1 () { return (*this)[0]; } PointIndex & I2 () { return (*this)[1]; } PointIndex & I3 () { return (*this)[2]; } diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp index a29c1271..611a61f3 100644 --- a/libsrc/meshing/secondorder.cpp +++ b/libsrc/meshing/secondorder.cpp @@ -424,7 +424,6 @@ namespace netgen { PrintMessage (3, "Validate mesh"); int np = mesh.GetNP(); - int ne = mesh.GetNE(); // int i, j; NgArray parents(np); From 3b3491a597affe1f2cc1da7eb703cda40e4bfc7b Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Thu, 2 Jan 2025 10:17:24 +0100 Subject: [PATCH 18/19] some little steps --- libsrc/general/hashtabl.cpp | 14 +++++++++++++- libsrc/general/hashtabl.hpp | 16 +++++++++++----- libsrc/general/template.hpp | 9 ++++----- libsrc/meshing/adfront2.hpp | 25 +++++++------------------ libsrc/meshing/meshtype.hpp | 6 ++++++ 5 files changed, 41 insertions(+), 29 deletions(-) diff --git a/libsrc/general/hashtabl.cpp b/libsrc/general/hashtabl.cpp index 30f6ed1f..c0be3c1c 100644 --- a/libsrc/general/hashtabl.cpp +++ b/libsrc/general/hashtabl.cpp @@ -292,7 +292,19 @@ namespace netgen - + BASE_INDEX_3_CLOSED_HASHTABLE :: + BASE_INDEX_3_CLOSED_HASHTABLE (size_t size) + : hash(RoundUp2(size)) + { + // cout << "orig size = " << size + // << ", roundup size = " << hash.Size(); + size = hash.Size(); + mask = size-1; + // cout << "mask = " << mask << endl; + invalid = -1; + for (size_t i = 0; i < size; i++) + hash[i].I1() = invalid; + } void BASE_INDEX_3_CLOSED_HASHTABLE :: diff --git a/libsrc/general/hashtabl.hpp b/libsrc/general/hashtabl.hpp index 110e4df4..07c0c7b2 100644 --- a/libsrc/general/hashtabl.hpp +++ b/libsrc/general/hashtabl.hpp @@ -861,9 +861,10 @@ inline ostream & operator<< (ostream & ost, const INDEX_2_CLOSED_HASHTABLE & for (int i = 0; i < ht.Size(); i++) if (ht.UsedPos(i)) { - INDEX_2 hash; - T data; - ht.GetData0 (i, hash, data); + // INDEX_2 hash; + // T data; + // ht.GetData0 (i, hash, data); + auto [hash,data] = ht.GetBoth(i); ost << "hash = " << hash << ", data = " << data << endl; } return ost; @@ -880,7 +881,8 @@ protected: size_t mask; protected: - BASE_INDEX_3_CLOSED_HASHTABLE (size_t size) + BASE_INDEX_3_CLOSED_HASHTABLE (size_t size); + /* : hash(RoundUp2(size)) { // cout << "orig size = " << size @@ -892,6 +894,7 @@ protected: for (size_t i = 0; i < size; i++) hash[i].I1() = invalid; } + */ public: int Size() const @@ -1073,9 +1076,12 @@ inline ostream & operator<< (ostream & ost, const INDEX_3_CLOSED_HASHTABLE & for (int i = 0; i < ht.Size(); i++) if (ht.UsedPos(i)) { + /* INDEX_3 hash; T data; - ht.GetData (i, hash, data); + ht.GetData (i, hash, data); + */ + auto [hash, data] = ht.GetBoth(); ost << "hash = " << hash << ", data = " << data << endl; } return ost; diff --git a/libsrc/general/template.hpp b/libsrc/general/template.hpp index a006a5c2..a3b9457d 100644 --- a/libsrc/general/template.hpp +++ b/libsrc/general/template.hpp @@ -114,8 +114,10 @@ class INDEX_2 public: /// + // protected: INDEX_2 () { } INDEX_2 (const INDEX_2&) = default; +public: INDEX_2 (INDEX_2&&) = default; INDEX_2 & operator= (const INDEX_2&) = default; @@ -157,7 +159,7 @@ public: return INDEX_2 (i1,i2); } - + operator std::array() { return { i[0], i[1] }; } /// INDEX & I1 () { return i[0]; } /// @@ -213,13 +215,10 @@ public: /// constexpr INDEX_3 (INDEX ai1, INDEX ai2, INDEX ai3) : i{ai1, ai2, ai3} { } - // { i[0] = ai1; i[1] = ai2; i[2] = ai3; } - /// + /// constexpr INDEX_3 (const INDEX_3 & in2) : i{in2.i[0], in2.i[1], in2.i[2]} { } - // { i[0] = in2.i[0]; i[1] = in2.i[1]; i[2] = in2.i[2]; } - static INDEX_3 Sort (INDEX_3 i3) { diff --git a/libsrc/meshing/adfront2.hpp b/libsrc/meshing/adfront2.hpp index 497580ca..d1d03472 100644 --- a/libsrc/meshing/adfront2.hpp +++ b/libsrc/meshing/adfront2.hpp @@ -95,7 +95,7 @@ namespace netgen { private: /// Point Indizes - INDEX_2 l; + INDEX_2 l; // want to replace by std::array l; /// quality class int lineclass; /// geometry specific data @@ -109,23 +109,12 @@ namespace netgen /// FrontLine (const INDEX_2 & al) - { - l = al; - lineclass = 1; - } - + : l(al), lineclass(1) { } /// - const INDEX_2 & L () const - { - return l; - } - + const auto & L () const { return l; } /// - int LineClass() const - { - return lineclass; - } + int LineClass() const { return lineclass; } /// void IncrementClass () @@ -141,13 +130,13 @@ namespace netgen /// bool Valid () const { - return l.I1() != -1; + return l[0] != -1; } /// void Invalidate () { - l.I1() = -1; - l.I2() = -1; + l[0] = -1; + l[1] = -1; lineclass = 1000; } diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 6f70248f..1dc2ac98 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -454,6 +454,12 @@ namespace netgen { public: using Index::Index; + /* + private: + operator int() const { return i; } + public: + operator ptrdiff_t() const { return i; } + */ }; inline istream & operator>> (istream & ist, ElementIndex & pi) From 643898c5e28de62a77637169b9ee9877c97999f2 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Thu, 2 Jan 2025 20:51:11 +0100 Subject: [PATCH 19/19] avoid shared ptr copy --- libsrc/meshing/meshclass.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index f89dddc1..4e05dbf2 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -507,7 +507,7 @@ namespace netgen /// LocalH & LocalHFunction (int layer=1) { return * lochfunc[layer-1]; } - shared_ptr GetLocalH(int layer=1) const + shared_ptr & GetLocalH(int layer=1) const { if(lochfunc.Size() == 1) return lochfunc[0];