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; + 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); } //(*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))); - - //(*testout) << "faces " << faces << endl; - - for(int i=0; i(); + for(auto face : fnrs) + faces.Append(topology.GetFace2SurfaceElement(face)); for(int i=0; iSize() != 0 && !indices->Contains(sel.GetIndex())) + if(indices && indices->Size() > 0 && !indices->Contains(sel.GetIndex())) continue; - auto & el = mesh.VolumeElement(velement); - + auto & el = mesh[velement]; if (el.GetType() == TET) { double lam4[4] = { vlam[0], vlam[1], vlam[2], 1.0-vlam[0]-vlam[1]-vlam[2] }; @@ -131,7 +125,7 @@ namespace netgen for(auto k : Range(4)) if(sel[j] == el[k]) lami[j-1] = lam4[k]/(1.0-face_lam); - return faces[i]; + return SurfaceElementIndex(faces[i]); } } @@ -143,8 +137,7 @@ namespace netgen // Did't find any matching face of a volume element, search 2d elements directly int ne; - NgArray locels; - // TODO: build search tree for surface elements + Array locels; if (searchtree) { searchtree->GetIntersecting (p, p, locels); @@ -153,37 +146,38 @@ namespace netgen else ne = mesh.GetNSE(); - for (int i = 1; i <= ne; i++) + for (auto i : Range(ne)) { - int ii; + SurfaceElementIndex ii; if (locels.Size()) - ii = locels.Get(i); + ii = locels[i]; else ii = i; - if(indices != NULL && indices->Size() > 0) + if(indices && indices->Size() > 0) { - bool contained = indices->Contains(mesh.SurfaceElement(ii).GetIndex()); + bool contained = indices->Contains(mesh[ii].GetIndex()); if((allowindex && !contained) || (!allowindex && contained)) continue; } - - if(mesh.PointContainedIn2DElement(p,lami,ii)) return ii; - + if(mesh.PointContainedIn2DElement(p,lami,ii)) + return ii; } return 0; } - int Find1dElement (const Mesh& mesh, - const netgen::Point<3> & p, - double * lami, - const NgArray * const indices, - BoxTree<3> * searchtree, - const bool allowindex = true) + SegmentIndex Find1dElement (const Mesh& mesh, + const netgen::Point<3> & p, + double * lami, + std::optional> indices, + BoxTree<3> * searchtree, + const bool allowindex = true) { double vlam[3]; - int velement = Find2dElement(mesh, p, vlam, NULL, searchtree, allowindex); - if(velement == 0) + if(searchtree) + const_cast(mesh).BuildElementSearchTree(2); + auto velement = Find2dElement(mesh, p, vlam, nullopt, searchtree ? mesh.GetSurfaceElementSearchTree() : nullptr, allowindex); + if(!velement.IsValid()) return 0; vlam[2] = 1.-vlam[0] - vlam[1]; @@ -5257,6 +5251,8 @@ namespace netgen void Mesh :: BuildElementSearchTree (int dim) { + if(dim < 2) + return; if (elementsearchtreets[dim] == GetTimeStamp()) return; @@ -5273,10 +5269,10 @@ namespace netgen Point3d pmin, pmax; GetBox(pmin, pmax); Box<3> box(pmin, pmax); - if (dim > 1) - elementsearchtree[dim] = make_unique>(box); + if (dim == 3) + elementsearchtree_vol = make_unique>(box); else - elementsearchtree[dim] = nullptr; // not yet implemented + elementsearchtree_surf = make_unique>(box); if (dim == 3) { @@ -5316,7 +5312,7 @@ namespace netgen } } box.Scale(1.2); - elementsearchtree[dim] -> Insert (box, ei+1); + elementsearchtree_vol -> Insert (box, ei); } } else if (dim == 2) @@ -5341,7 +5337,7 @@ namespace netgen } box.Scale(1.2); } - elementsearchtree[dim] -> Insert (box, ei+1); + elementsearchtree_surf -> Insert (box, ei); } } } @@ -5381,7 +5377,7 @@ namespace netgen bool Mesh :: PointContainedIn2DElement(const Point3d & p, double lami[3], - const int element, + SurfaceElementIndex ei, bool consider3D) const { Vec3d col1, col2, col3; @@ -5392,9 +5388,9 @@ namespace netgen //SZ - if(SurfaceElement(element).GetType()==QUAD) + if(surfelements[ei].GetType()==QUAD) { - const Element2d & el = SurfaceElement(element); + const Element2d & el = surfelements[ei]; const Point3d & p1 = Point(el.PNum(1)); const Point3d & p2 = Point(el.PNum(2)); @@ -5413,7 +5409,7 @@ namespace netgen int i = 0; while(delta > 1e-16 && i < maxits) { - curvedelems->CalcSurfaceTransformation(lam,element-1,x,Jac); + curvedelems->CalcSurfaceTransformation(lam,ei,x,Jac); rhs = p - x; Jac.Solve(rhs,deltalam); lam += deltalam; @@ -5742,7 +5738,7 @@ namespace netgen { // SurfaceElement(element).GetTets (loctets); loctrigs.SetSize(1); - loctrigs.Elem(1) = SurfaceElement(element); + loctrigs.Elem(1) = surfelements[ei]; @@ -5777,11 +5773,11 @@ namespace netgen //(*testout) << "col1 " << col1 << " col2 " << col2 << " col3 " << col3 << " rhs " << rhs << endl; //(*testout) << "sol " << sol << endl; - if (SurfaceElement(element).GetType() ==TRIG6 || curvedelems->IsSurfaceElementCurved(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 +5796,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 +5815,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 +5847,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 +5858,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 +5893,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 +6015,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 73a2145a..75b65c66 100644 --- a/libsrc/visualization/vssolution.cpp +++ b/libsrc/visualization/vssolution.cpp @@ -4341,7 +4341,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) { @@ -4862,18 +4862,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; @@ -4884,7 +4884,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);