From 7aae5369c455c04933e0fa7af87511c09352e915 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Thu, 13 Mar 2025 18:39:21 +0100 Subject: [PATCH 1/2] move all searchtrees to use elementindex --- libsrc/interface/nginterface.cpp | 6 +- libsrc/interface/nginterface_v2.cpp | 27 ++- libsrc/interface/writefluent.cpp | 11 +- libsrc/meshing/fieldlines.cpp | 30 +++- libsrc/meshing/meshclass.cpp | 264 ++++++++++++++-------------- libsrc/meshing/meshclass.hpp | 77 ++++---- libsrc/meshing/meshfunc.cpp | 7 +- libsrc/visualization/vssolution.cpp | 10 +- 8 files changed, 222 insertions(+), 210 deletions(-) diff --git a/libsrc/interface/nginterface.cpp b/libsrc/interface/nginterface.cpp index 3db2fa1c..9186a5ef 100644 --- a/libsrc/interface/nginterface.cpp +++ b/libsrc/interface/nginterface.cpp @@ -651,14 +651,14 @@ int Ng_FindElementOfPoint (double * p, double * lami, int build_searchtree, { Point3d p3d(p[0], p[1], p[2]); ind = - mesh->GetElementOfPoint(p3d, lami, dummy, build_searchtree != 0); + mesh->GetElementOfPoint(p3d, lami, dummy, build_searchtree != 0) + 1; } else { double lam3[3]; Point3d p2d(p[0], p[1], 0); ind = - mesh->GetElementOfPoint(p2d, lam3, dummy, build_searchtree != 0); + mesh->GetSurfaceElementOfPoint(p2d, lam3, dummy, build_searchtree != 0) + 1; if (ind > 0) { @@ -697,7 +697,7 @@ int Ng_FindSurfaceElementOfPoint (double * p, double * lami, int build_searchtre { Point3d p3d(p[0], p[1], p[2]); ind = - mesh->GetSurfaceElementOfPoint(p3d, lami, dummy, build_searchtree != 0); + mesh->GetSurfaceElementOfPoint(p3d, lami, dummy, build_searchtree != 0) + 1; } else { diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp index d238f401..2c95505f 100644 --- a/libsrc/interface/nginterface_v2.cpp +++ b/libsrc/interface/nginterface_v2.cpp @@ -1043,13 +1043,26 @@ namespace netgen int * const indices, int numind) const { - if(build_searchtree) - mesh->BuildElementSearchTree(2); Point<3> pp(p[0], p[1], 0.); if(mesh->GetDimension() == 3) pp[2] = p[2]; - NgArray ind(numind, indices); - return Find2dElement(*mesh, pp, lami, &ind, mesh->GetElementSearchTree(2)); + FlatArray ind(numind, indices); + double lam3[3]; + auto elnr = mesh->GetSurfaceElementOfPoint(pp, lam3, ind, build_searchtree); + if(elnr.IsValid()) + { + if((*mesh)[elnr].GetType() == QUAD || (*mesh)[elnr].GetType() == TRIG6) + { + lami[0] = lam3[0]; + lami[1] = lam3[1]; + } + else + { + lami[0] = 1-lam3[0]-lam3[1]; + lami[1] = lam3[0]; + } + } + return elnr; } @@ -1060,11 +1073,9 @@ namespace netgen int * const indices, int numind) const { - if(build_searchtree) - mesh->BuildElementSearchTree(3); Point<3> pp(p[0], p[1], p[2]); - NgArray ind(numind, indices); - return Find3dElement(*mesh, pp, lami, &ind, mesh->GetElementSearchTree(3)); + FlatArray ind(numind, indices); + return mesh->GetElementOfPoint(pp, lami, ind, build_searchtree); } void Ngx_Mesh :: Curve (int order) diff --git a/libsrc/interface/writefluent.cpp b/libsrc/interface/writefluent.cpp index 78141c7b..6cf64b72 100644 --- a/libsrc/interface/writefluent.cpp +++ b/libsrc/interface/writefluent.cpp @@ -66,7 +66,7 @@ void WriteFluentFormat (const Mesh & mesh, int i2, j2; NgArray surfaceelp; NgArray surfaceeli; - NgArray locels; + Array locels; //no cells=no tets //no faces=2*tets @@ -117,12 +117,9 @@ void WriteFluentFormat (const Mesh & mesh, int eli2 = 0; int stopsig = 0; - for (i2 = 1; i2 <= nel; i2++) + for (auto locind : locels) { - locind = locels.Get(i2); - //cout << " locind=" << locind << endl; - - Element el2 = mesh.VolumeElement(locind); + Element el2 = mesh[locind]; //if (inverttets) // el2.Invert(); @@ -130,7 +127,7 @@ void WriteFluentFormat (const Mesh & mesh, { el2.GetFace(j2, face2); - if (face2.HasFace(face)) {eli2 = locind; stopsig = 1; break;} + if (face2.HasFace(face)) {eli2 = locind+1; stopsig = 1; break;} } if (stopsig) break; } diff --git a/libsrc/meshing/fieldlines.cpp b/libsrc/meshing/fieldlines.cpp index d2ee2e21..35dfa0fe 100644 --- a/libsrc/meshing/fieldlines.cpp +++ b/libsrc/meshing/fieldlines.cpp @@ -9,6 +9,25 @@ namespace netgen { + inline int GetVolElement(const Mesh& mesh, const Point<3>& p, + double* lami) + { + if(mesh.GetDimension() == 3) + { + auto ei = mesh.GetElementOfPoint(p, lami, true); + if(!ei.IsValid()) + return -1; + return ei; + } + else + { + auto ei = mesh.GetSurfaceElementOfPoint(p, lami, true); + if(!ei.IsValid()) + return -1; + return ei; + } + } + RKStepper :: ~RKStepper() { delete a; @@ -190,9 +209,9 @@ namespace netgen for(int i=0; i & p, - double * lami, - const NgArray * const indices, - BoxTree<3> * searchtree, - const bool allowindex) + SurfaceElementIndex + Find2dElement (const Mesh& mesh, + const netgen::Point<3> & p, + double * lami, + std::optional> indices, + BoxTree<3, SurfaceElementIndex> * searchtree, + bool allowindex) { double vlam[3]; - int velement = 0; + cout << "find2delement " << p << endl; + ElementIndex velement = ElementIndex::INVALID; if(mesh.GetNE()) { if(searchtree) const_cast(mesh).BuildElementSearchTree(3); - velement = Find3dElement(mesh, p,vlam,NULL,searchtree ? mesh.GetElementSearchTree(3) : nullptr,allowindex); + velement = Find3dElement(mesh, p,vlam, nullopt,searchtree ? mesh.GetElementSearchTree() : nullptr,allowindex); + cout << "found volume element = " << velement << endl; } //(*testout) << "p " << p << endl; //(*testout) << "velement " << velement << endl; // first try to find a volume element containing p and project to face - if(velement!=0) + if(velement.IsValid()) { auto & topology = mesh.GetTopology(); - // NgArray faces; - // topology.GetElementFaces(velement,faces); - auto faces = Array (topology.GetFaces(ElementIndex(velement-1))); + const auto & fnrs = topology.GetFaces(velement); + auto faces = ArrayMem(); + for(auto face : fnrs) + faces.Append(topology.GetFace2SurfaceElement(face)); - //(*testout) << "faces " << faces << endl; - - for(int i=0; iIsSurfaceElementCurved(element-1)) + if (surfelements[ei].GetType() ==TRIG6 || curvedelems->IsSurfaceElementCurved(ei)) { // netgen::Point<2> lam(1./3,1./3); netgen::Point<2> lam(sol.X(), sol.Y()); - if(SurfaceElement(element).GetType() != TRIG6) + if(surfelements[ei].GetType() != TRIG6) { lam[0] = 1-sol.X()-sol.Y(); lam[1] = sol.X(); @@ -5800,7 +5811,7 @@ namespace netgen const int maxits = 30; while(delta > 1e-16 && iCalcSurfaceTransformation(lam,element-1,x,Jac); + curvedelems->CalcSurfaceTransformation(lam,ei,x,Jac); rhs = p-x; Jac.Solve(rhs,deltalam); @@ -5819,7 +5830,7 @@ namespace netgen sol.X() = lam(0); sol.Y() = lam(1); - if (SurfaceElement(element).GetType() !=TRIG6 ) + if (surfelements[ei].GetType() !=TRIG6 ) { sol.Z() = sol.X(); sol.X() = sol.Y(); @@ -5851,7 +5862,7 @@ namespace netgen bool Mesh :: PointContainedIn3DElement(const Point3d & p, double lami[3], - const int element) const + ElementIndex ei) const { //bool oldresult = PointContainedIn3DElementOld(p,lami,element); //(*testout) << "old result: " << oldresult @@ -5862,7 +5873,7 @@ namespace netgen const double eps = 1.e-4; - const Element & el = VolumeElement(element); + const Element & el = volelements[ei]; netgen::Point<3> lam = 0.0; @@ -5897,7 +5908,7 @@ namespace netgen const int maxits = 30; while(delta > 1e-16 && iCalcElementTransformation(lam,element-1,x,Jac); + curvedelems->CalcElementTransformation(lam,ei,x,Jac); rhs = p-x; Jac.Solve(rhs,deltalam); @@ -6019,91 +6030,72 @@ namespace netgen } - int Mesh :: GetElementOfPoint (const netgen::Point<3> & p, - double lami[3], - bool build_searchtree, - const int index, - const bool allowindex) const + ElementIndex Mesh :: GetElementOfPoint (const netgen::Point<3> & p, + double* lami, + bool build_searchtree, + int index, + bool allowindex) const { if(index != -1) { - NgArray dummy(1); + Array dummy(1); dummy[0] = index; - return GetElementOfPoint(p,lami,&dummy,build_searchtree,allowindex); + return GetElementOfPoint(p,lami,dummy,build_searchtree,allowindex); } else - return GetElementOfPoint(p,lami,NULL,build_searchtree,allowindex); + return GetElementOfPoint(p,lami,nullopt,build_searchtree,allowindex); } - int Mesh :: GetElementOfPoint (const netgen::Point<3> & p, - double lami[3], - const NgArray * const indices, - bool build_searchtree, - const bool allowindex) const + ElementIndex Mesh :: GetElementOfPoint (const netgen::Point<3> & p, + double* lami, + std::optional> indices, + bool build_searchtree, + bool allowindex) const { - if ( (dimension == 2 && !GetNSE()) || - (dimension == 3 && !GetNE() && !GetNSE()) ) - return -1; - - - if (dimension == 2 || (dimension==3 && !GetNE() && GetNSE())) - { - if (build_searchtree) - const_cast(*this).BuildElementSearchTree (2); - return Find2dElement(*this, p, lami, indices, elementsearchtree[2].get(), allowindex); - } - if (build_searchtree) const_cast(*this).BuildElementSearchTree (3); - return Find3dElement(*this, p, lami, indices, elementsearchtree[3].get(), allowindex); + return Find3dElement(*this, p, lami, indices, elementsearchtree_vol.get(), allowindex); } - int Mesh :: GetSurfaceElementOfPoint (const netgen::Point<3> & p, - double lami[3], - bool build_searchtree, - const int index, - const bool allowindex) const + SurfaceElementIndex Mesh :: + GetSurfaceElementOfPoint (const netgen::Point<3> & p, + double* lami, + bool build_searchtree, + int index, + bool allowindex) const { - if(index != -1) + if(index != -1) { - NgArray dummy(1); + Array dummy(1); dummy[0] = index; - return GetSurfaceElementOfPoint(p,lami,&dummy,build_searchtree,allowindex); + return GetSurfaceElementOfPoint(p,lami,dummy,build_searchtree,allowindex); } else - return GetSurfaceElementOfPoint(p,lami,NULL,build_searchtree,allowindex); + return GetSurfaceElementOfPoint(p,lami,nullopt,build_searchtree,allowindex); } - - - - int Mesh :: GetSurfaceElementOfPoint (const netgen::Point<3> & p, - double lami[3], - const NgArray * const indices, - bool build_searchtree, - const bool allowindex) const + SurfaceElementIndex Mesh :: + GetSurfaceElementOfPoint (const netgen::Point<3> & p, + double* lami, + std::optional> indices, + bool build_searchtree, + bool allowindex) const { - if (dimension == 2) - return Find1dElement(*this, p, lami, indices, nullptr, allowindex); - else - { - if (build_searchtree) - const_cast(*this).BuildElementSearchTree(2); - return Find2dElement(*this, p, lami, indices, elementsearchtree[2].get(), allowindex); - } - return 0; + if (build_searchtree) + const_cast(*this).BuildElementSearchTree(2); + return Find2dElement(*this, p, lami, indices, elementsearchtree_surf.get(), allowindex); } void Mesh::GetIntersectingVolEls(const Point3d& p1, const Point3d& p2, - NgArray & locels) const + Array & locels) const { - elementsearchtree[3]->GetIntersecting (p1, p2, locels); + elementsearchtree_vol->GetIntersecting (p1, p2, locels); } void Mesh :: SplitIntoParts() diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index 6e52f629..6188debe 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -148,7 +148,8 @@ namespace netgen int numvertices; /// geometric search tree for interval intersection search - unique_ptr> elementsearchtree[4] = {nullptr,nullptr,nullptr,nullptr}; + unique_ptr> elementsearchtree_vol; + unique_ptr> elementsearchtree_surf; /// time stamp for tree mutable size_t elementsearchtreets[4]; @@ -200,11 +201,11 @@ namespace netgen DLL_HEADER bool PointContainedIn2DElement(const Point3d & p, double lami[3], - const int element, + SurfaceElementIndex element, bool consider3D = false) const; DLL_HEADER bool PointContainedIn3DElement(const Point3d & p, double lami[3], - const int element) const; + ElementIndex element) const; DLL_HEADER bool PointContainedIn3DElementOld(const Point3d & p, double lami[3], const int element) const; @@ -664,38 +665,47 @@ namespace netgen /// build box-search tree DLL_HEADER void BuildElementSearchTree (int dim); - BoxTree<3>* GetElementSearchTree (int dim) const + BoxTree<3, ElementIndex>* GetElementSearchTree () const { - return elementsearchtree[dim].get(); + return elementsearchtree_vol.get(); + } + + BoxTree<3, SurfaceElementIndex>* GetSurfaceElementSearchTree () const + { + return elementsearchtree_surf.get(); } void SetPointSearchStartElement(const int el) const {ps_startelement = el;} /// gives element of point, barycentric coordinates - DLL_HEADER int GetElementOfPoint (const netgen::Point<3> & p, - double * lami, - bool build_searchtree = 0, - const int index = -1, - const bool allowindex = true) const; - DLL_HEADER int GetElementOfPoint (const netgen::Point<3> & p, - double * lami, - const NgArray * const indices, - bool build_searchtree = 0, - const bool allowindex = true) const; - DLL_HEADER int GetSurfaceElementOfPoint (const netgen::Point<3> & p, - double * lami, - bool build_searchtree = 0, - const int index = -1, - const bool allowindex = true) const; - DLL_HEADER int GetSurfaceElementOfPoint (const netgen::Point<3> & p, - double * lami, - const NgArray * const indices, - bool build_searchtree = 0, - const bool allowindex = true) const; + DLL_HEADER ElementIndex + GetElementOfPoint (const netgen::Point<3> & p, + double * lami, + bool build_searchtree = false, + int index = -1, + bool allowindex = true) const; + DLL_HEADER ElementIndex + GetElementOfPoint (const netgen::Point<3> & p, + double * lami, + std::optional> indices, + bool build_searchtree = 0, + bool allowindex = true) const; + DLL_HEADER SurfaceElementIndex + GetSurfaceElementOfPoint (const netgen::Point<3> & p, + double * lami, + bool build_searchtree = false, + int index = -1, + bool allowindex = true) const; + DLL_HEADER SurfaceElementIndex + GetSurfaceElementOfPoint (const netgen::Point<3> & p, + double * lami, + std::optional> indices, + bool build_searchtree = false, + bool allowindex = true) const; /// give list of vol elements which are int the box(p1,p2) void GetIntersectingVolEls(const Point3d& p1, const Point3d& p2, - NgArray & locels) const; + Array & locels) const; /// int AddFaceDescriptor(const FaceDescriptor& fd) @@ -1040,21 +1050,6 @@ namespace netgen } DLL_HEADER void AddFacesBetweenDomains(Mesh & mesh); - - int Find2dElement (const Mesh& mesh, - const netgen::Point<3> & p, - double * lami, - const NgArray * const indices, - BoxTree<3> * searchtree, - const bool allowindex = true); - - int Find3dElement (const Mesh& mesh, - const netgen::Point<3> & p, - double * lami, - const NgArray * const indices, - BoxTree<3> * searchtree, - const bool allowindex = true); - } #endif // NETGEN_MESHCLASS_HPP diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index b0f6d6db..5dee0579 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -808,17 +808,16 @@ namespace netgen double lam[3]; ElementIndex ei_start = mesh.GetElementOfPoint(p_new, lam, false, domain); - if(ei_start == 0) { + if(!ei_start.IsValid()) { PrintMessage(1, "Could not find volume element with new point"); return; } - ei_start -= 1; if(mesh[ei_start].IsDeleted()) return; double max_inside = -1.; - ElementIndex ei_max_inside = -1; + ElementIndex ei_max_inside = ElementIndex::INVALID; // search for adjacent volume element, where the new point is "most inside", // i.e. the minimal barycentric coordinate is maximal @@ -828,7 +827,7 @@ namespace netgen if(mesh[ei1].IsDeleted()) return; - if(!mesh.PointContainedIn3DElement(p_new, lam, ei1+1)) + if(!mesh.PointContainedIn3DElement(p_new, lam, ei1)) continue; double inside = min(min(lam[0], lam[1]), min(lam[2], 1.0-lam[0]-lam[1])); diff --git a/libsrc/visualization/vssolution.cpp b/libsrc/visualization/vssolution.cpp index 7560aea6..47459941 100644 --- a/libsrc/visualization/vssolution.cpp +++ b/libsrc/visualization/vssolution.cpp @@ -4337,7 +4337,7 @@ namespace netgen } double lami[3]; - int elnr = mesh->GetElementOfPoint (hp, lami,0,cindex,allowindex)-1; + int elnr = mesh->GetElementOfPoint (hp, lami,0,cindex,allowindex); if (elnr != -1) { @@ -4858,18 +4858,18 @@ namespace netgen if(sol.iscomplex && rcomponent != 0) { rcomponent = 2 * ((rcomponent-1)/2) + 1; - GetValue(&sol, el3d-1, lami[0], lami[1], lami[2], rcomponent+1, + GetValue(&sol, el3d, lami[0], lami[1], lami[2], rcomponent+1, imag); comp = (scalcomp-1)/2 + 1; } - GetValue(&sol, el3d-1, lami[0], lami[1], lami[2], rcomponent, val); + GetValue(&sol, el3d, lami[0], lami[1], lami[2], rcomponent, val); printScalValue(sol, comp, val, imag, sol.iscomplex && comp > 0); } if(vecfunction!=-1 && soldata[vecfunction]->draw_volume) { auto & sol = *soldata[vecfunction]; ArrayMem values(sol.components); - GetValues(&sol, el3d-1, lami[0], lami[1], lami[2], &values[0]); + GetValues(&sol, el3d, lami[0], lami[1], lami[2], &values[0]); printVecValue(sol, values); } return; @@ -4880,7 +4880,7 @@ namespace netgen double lami[3] = {0.0, 0.0, 0.0}; // Check if unprojected Point is close to surface element (eps of 1e-3 due to z-Buffer accuracy) bool found_2del = false; - if(selelement>0 && mesh->PointContainedIn2DElement(p, lami, selelement, false && fabs(lami[2])<1e-3)) + if(selelement>0 && mesh->PointContainedIn2DElement(p, lami, selelement-1, false && fabs(lami[2])<1e-3)) { // Found it, use coordinates of point projected to surface element mesh->GetCurvedElements().CalcSurfaceTransformation({1.0-lami[0]-lami[1], lami[0]}, selelement-1, p); From b8d722d6a8347218681407d87114f4836e0a3714 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Thu, 13 Mar 2025 18:41:38 +0100 Subject: [PATCH 2/2] remove debug output --- libsrc/meshing/meshclass.cpp | 19 ++----------------- 1 file changed, 2 insertions(+), 17 deletions(-) diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 7e9209d3..cecce6cd 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -83,7 +83,6 @@ namespace netgen bool allowindex) { double vlam[3]; - cout << "find2delement " << p << endl; ElementIndex velement = ElementIndex::INVALID; if(mesh.GetNE()) @@ -91,7 +90,6 @@ namespace netgen if(searchtree) const_cast(mesh).BuildElementSearchTree(3); velement = Find3dElement(mesh, p,vlam, nullopt,searchtree ? mesh.GetElementSearchTree() : nullptr,allowindex); - cout << "found volume element = " << velement << endl; } //(*testout) << "p " << p << endl; @@ -106,8 +104,6 @@ namespace netgen for(auto face : fnrs) faces.Append(topology.GetFace2SurfaceElement(face)); - cout << "faces = " << faces << endl; - for(int i=0; i