ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-03-18 12:37:56 +05:00
Code Issues Packages Projects Releases Wiki Activity

move all searchtrees to use elementindex

Browse Source
This commit is contained in:
Christopher Lackner 2025-03-13 18:39:21 +01:00
parent 787c6043fa
commit 7aae5369c4
8 changed files with 222 additions and 210 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
libsrc/interface/nginterface.cpp
View File

@ -651,14 +651,14 @@ int Ng_FindElementOfPoint (double * p, double * lami, int build_searchtree,
{ {
Point3d p3d(p[0], p[1], p[2]); Point3d p3d(p[0], p[1], p[2]);
ind = ind =
mesh->GetElementOfPoint(p3d, lami, dummy, build_searchtree != 0); mesh->GetElementOfPoint(p3d, lami, dummy, build_searchtree != 0) + 1;
} }
else else
{ {
double lam3[3]; double lam3[3];
Point3d p2d(p[0], p[1], 0); Point3d p2d(p[0], p[1], 0);
ind = ind =
mesh->GetElementOfPoint(p2d, lam3, dummy, build_searchtree != 0); mesh->GetSurfaceElementOfPoint(p2d, lam3, dummy, build_searchtree != 0) + 1;
if (ind > 0) 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]); Point3d p3d(p[0], p[1], p[2]);
ind = ind =
mesh->GetSurfaceElementOfPoint(p3d, lami, dummy, build_searchtree != 0); mesh->GetSurfaceElementOfPoint(p3d, lami, dummy, build_searchtree != 0) + 1;
} }
else else
{ {

27
libsrc/interface/nginterface_v2.cpp
View File

@ -1043,13 +1043,26 @@ namespace netgen
int * const indices, int numind) const int * const indices, int numind) const
{ {
if(build_searchtree)
mesh->BuildElementSearchTree(2);
Point<3> pp(p[0], p[1], 0.); Point<3> pp(p[0], p[1], 0.);
if(mesh->GetDimension() == 3) if(mesh->GetDimension() == 3)
pp[2] = p[2]; pp[2] = p[2];
NgArray<int> ind(numind, indices); FlatArray<int> ind(numind, indices);
return Find2dElement(*mesh, pp, lami, &ind, mesh->GetElementSearchTree(2)); 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 int * const indices, int numind) const
{ {
if(build_searchtree)
mesh->BuildElementSearchTree(3);
Point<3> pp(p[0], p[1], p[2]); Point<3> pp(p[0], p[1], p[2]);
NgArray<int> ind(numind, indices); FlatArray<int> ind(numind, indices);
return Find3dElement(*mesh, pp, lami, &ind, mesh->GetElementSearchTree(3)); return mesh->GetElementOfPoint(pp, lami, ind, build_searchtree);
} }
void Ngx_Mesh :: Curve (int order) void Ngx_Mesh :: Curve (int order)

11
libsrc/interface/writefluent.cpp
View File

@ -66,7 +66,7 @@ void WriteFluentFormat (const Mesh & mesh,
int i2, j2; int i2, j2;
NgArray<INDEX_3> surfaceelp; NgArray<INDEX_3> surfaceelp;
NgArray<int> surfaceeli; NgArray<int> surfaceeli;
NgArray<int> locels; Array<ElementIndex> locels;
//no cells=no tets //no cells=no tets
//no faces=2*tets //no faces=2*tets
@ -117,12 +117,9 @@ void WriteFluentFormat (const Mesh & mesh,
int eli2 = 0; int eli2 = 0;
int stopsig = 0; int stopsig = 0;
for (i2 = 1; i2 <= nel; i2++) for (auto locind : locels)
{ {
locind = locels.Get(i2); Element el2 = mesh[locind];
//cout << " locind=" << locind << endl;
Element el2 = mesh.VolumeElement(locind);
//if (inverttets) //if (inverttets)
// el2.Invert(); // el2.Invert();
@ -130,7 +127,7 @@ void WriteFluentFormat (const Mesh & mesh,
{ {
el2.GetFace(j2, face2); 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; if (stopsig) break;
} }

30
libsrc/meshing/fieldlines.cpp
View File

@ -9,6 +9,25 @@
namespace netgen 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() RKStepper :: ~RKStepper()
{ {
delete a; delete a;
@ -190,9 +209,9 @@ namespace netgen
for(int i=0; i<potential_startpoints.Size(); i++) for(int i=0; i<potential_startpoints.Size(); i++)
{ {
int elnr = mesh.GetElementOfPoint(potential_startpoints[i],lami,true) - 1; int elnr = GetVolElement(mesh, potential_startpoints[i], lami);
if(elnr == -1) if (elnr == -1)
continue; continue;
mesh.SetPointSearchStartElement(elnr); mesh.SetPointSearchStartElement(elnr);
@ -289,8 +308,7 @@ namespace netgen
dirstart.SetSize(0); dirstart.SetSize(0);
dirstart.Append(0); dirstart.Append(0);
int startelnr = GetVolElement(mesh, startpoint,startlami);
int startelnr = mesh.GetElementOfPoint(startpoint,startlami,true) - 1;
(*testout) << "p = " << startpoint << "; elnr = " << startelnr << endl; (*testout) << "p = " << startpoint << "; elnr = " << startelnr << endl;
if (startelnr == -1) if (startelnr == -1)
return; return;
@ -342,7 +360,7 @@ namespace netgen
Point<3> newp; Point<3> newp;
while(!stepper.GetNextData(newp,dummyt,h) && elnr != -1) while(!stepper.GetNextData(newp,dummyt,h) && elnr != -1)
{ {
elnr = mesh.GetElementOfPoint(newp,lami,true) - 1; elnr = GetVolElement(mesh, newp, lami);
if(elnr != -1) if(elnr != -1)
{ {
mesh.SetPointSearchStartElement(elnr); mesh.SetPointSearchStartElement(elnr);

264
libsrc/meshing/meshclass.cpp
View File

@ -12,15 +12,15 @@
namespace netgen namespace netgen
{ {
int Find3dElement (const Mesh& mesh, ElementIndex Find3dElement (const Mesh& mesh,
const netgen::Point<3> & p, const netgen::Point<3> & p,
double * lami, double * lami,
const NgArray<int> * const indices, optional<FlatArray<int>> indices,
BoxTree<3> * searchtree, BoxTree<3, ElementIndex> * searchtree,
const bool allowindex) const bool allowindex)
{ {
int ne = 0; int ne = 0;
NgArray<int> locels; Array<ElementIndex> locels;
if (searchtree) if (searchtree)
{ {
searchtree->GetIntersecting (p, p, locels); searchtree->GetIntersecting (p, p, locels);
@ -29,100 +29,101 @@ namespace netgen
else else
ne = mesh.GetNE(); ne = mesh.GetNE();
for (int i = 1; i <= ne; i++) for (auto i : Range(ne))
{ {
int ii; ElementIndex ei;
if (searchtree) if (searchtree)
ii = locels.Get(i); ei = locels[i];
else else
ii = i; ei = i;
if(indices != NULL && indices->Size() > 0) if(indices && indices->Size() > 0)
{ {
bool contained = indices->Contains(mesh.VolumeElement(ii).GetIndex()); bool contained = indices->Contains(mesh[ei].GetIndex());
if((allowindex && !contained) || (!allowindex && contained)) continue; if((allowindex && !contained) || (!allowindex && contained)) continue;
} }
if(mesh.PointContainedIn3DElement(p,lami,ii)) if(mesh.PointContainedIn3DElement(p,lami,ei))
return ii; return ei;
} }
// Not found, try uncurved variant: // Not found, try uncurved variant:
for (int i = 1; i <= ne; i++) for (auto i : Range(ne))
{ {
int ii; ElementIndex ei;
if (searchtree) if (searchtree)
ii = locels.Get(i); ei = locels[i];
else else
ii = i; ei = i;
if(indices != NULL && indices->Size() > 0) if(indices && indices->Size() > 0)
{ {
bool contained = indices->Contains(mesh.VolumeElement(ii).GetIndex()); bool contained = indices->Contains(mesh[ei].GetIndex());
if((allowindex && !contained) || (!allowindex && contained)) continue; if((allowindex && !contained) || (!allowindex && contained)) continue;
} }
if(mesh.PointContainedIn3DElementOld(p,lami,ii)) if(mesh.PointContainedIn3DElementOld(p,lami,ei))
{ {
(*testout) << "WARNING: found element of point " << p <<" only for uncurved mesh" << endl; (*testout) << "WARNING: found element of point " << p <<" only for uncurved mesh" << endl;
return ii; return ei;
} }
} }
return 0; return ElementIndex::INVALID;
} }
int Find2dElement (const Mesh& mesh, SurfaceElementIndex
const netgen::Point<3> & p, Find2dElement (const Mesh& mesh,
double * lami, const netgen::Point<3> & p,
const NgArray<int> * const indices, double * lami,
BoxTree<3> * searchtree, std::optional<FlatArray<int>> indices,
const bool allowindex) BoxTree<3, SurfaceElementIndex> * searchtree,
bool allowindex)
{ {
double vlam[3]; double vlam[3];
int velement = 0; cout << "find2delement " << p << endl;
ElementIndex velement = ElementIndex::INVALID;
if(mesh.GetNE()) if(mesh.GetNE())
{ {
if(searchtree) if(searchtree)
const_cast<Mesh&>(mesh).BuildElementSearchTree(3); const_cast<Mesh&>(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) << "p " << p << endl;
//(*testout) << "velement " << velement << endl; //(*testout) << "velement " << velement << endl;
// first try to find a volume element containing p and project to face // first try to find a volume element containing p and project to face
if(velement!=0) if(velement.IsValid())
{ {
auto & topology = mesh.GetTopology(); auto & topology = mesh.GetTopology();
// NgArray<int> faces; const auto & fnrs = topology.GetFaces(velement);
// topology.GetElementFaces(velement,faces); auto faces = ArrayMem<SurfaceElementIndex,4>();
auto faces = Array<int> (topology.GetFaces(ElementIndex(velement-1))); for(auto face : fnrs)
faces.Append(topology.GetFace2SurfaceElement(face));
//(*testout) << "faces " << faces << endl; cout << "faces = " << faces << endl;
for(int i=0; i<faces.Size(); i++)
faces[i] = topology.GetFace2SurfaceElement(faces[i])+1;
//(*testout) << "surfel " << faces << endl;
for(int i=0; i<faces.Size(); i++) for(int i=0; i<faces.Size(); i++)
{ {
if(faces[i] == 0) if(!faces[i].IsValid())
continue; continue;
auto sel = mesh.SurfaceElement(faces[i]); auto sel = mesh.SurfaceElement(faces[i]);
if(indices && indices->Size() != 0 && !indices->Contains(sel.GetIndex())) if(indices && indices->Size() > 0 && !indices->Contains(sel.GetIndex()))
continue; continue;
auto & el = mesh.VolumeElement(velement); auto & el = mesh[velement];
cout << "eltype = " << el.GetType() << endl;
if (el.GetType() == TET) if (el.GetType() == TET)
{ {
double lam4[4] = { vlam[0], vlam[1], vlam[2], 1.0-vlam[0]-vlam[1]-vlam[2] }; double lam4[4] = { vlam[0], vlam[1], vlam[2], 1.0-vlam[0]-vlam[1]-vlam[2] };
cout << "lam4 = " << lam4[0] << ", " << lam4[1] << ", " << lam4[2] << ", " << lam4[3] << endl;
double face_lam = lam4[i]; double face_lam = lam4[i];
cout << "face lam = " << face_lam << endl;
if(face_lam < 1e-5) if(face_lam < 1e-5)
{ {
// found volume point very close to a face -> use barycentric coordinates directly // found volume point very close to a face -> use barycentric coordinates directly
@ -131,20 +132,23 @@ namespace netgen
for(auto k : Range(4)) for(auto k : Range(4))
if(sel[j] == el[k]) if(sel[j] == el[k])
lami[j-1] = lam4[k]/(1.0-face_lam); lami[j-1] = lam4[k]/(1.0-face_lam);
return faces[i]; cout << "found close tet face = " << faces[i] << ", sel = " << sel << endl;
return SurfaceElementIndex(faces[i]);
} }
} }
if(mesh.PointContainedIn2DElement(p,lami,faces[i],true)) if(mesh.PointContainedIn2DElement(p,lami,faces[i],true))
return faces[i]; {
cout << "apoint contained in 2d el = " << faces[i] << ", sel = " << sel << endl;
return faces[i];
}
} }
} }
// Did't find any matching face of a volume element, search 2d elements directly // Did't find any matching face of a volume element, search 2d elements directly
int ne; int ne;
NgArray<int> locels; Array<SurfaceElementIndex> locels;
// TODO: build search tree for surface elements
if (searchtree) if (searchtree)
{ {
searchtree->GetIntersecting (p, p, locels); searchtree->GetIntersecting (p, p, locels);
@ -153,37 +157,42 @@ namespace netgen
else else
ne = mesh.GetNSE(); ne = mesh.GetNSE();
for (int i = 1; i <= ne; i++) for (auto i : Range(ne))
{ {
int ii; SurfaceElementIndex ii;
if (locels.Size()) if (locels.Size())
ii = locels.Get(i); ii = locels[i];
else else
ii = i; 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((allowindex && !contained) || (!allowindex && contained)) continue;
} }
if(mesh.PointContainedIn2DElement(p,lami,ii))
if(mesh.PointContainedIn2DElement(p,lami,ii)) return ii; {
cout << "point contained in 2d el = " << ii << ", sel = " << mesh[ii] << endl;
return ii;
}
} }
return 0; return 0;
} }
int Find1dElement (const Mesh& mesh, SegmentIndex Find1dElement (const Mesh& mesh,
const netgen::Point<3> & p, const netgen::Point<3> & p,
double * lami, double * lami,
const NgArray<int> * const indices, std::optional<FlatArray<int>> indices,
BoxTree<3> * searchtree, BoxTree<3> * searchtree,
const bool allowindex = true) const bool allowindex = true)
{ {
double vlam[3]; double vlam[3];
int velement = Find2dElement(mesh, p, vlam, NULL, searchtree, allowindex); if(searchtree)
if(velement == 0) const_cast<Mesh&>(mesh).BuildElementSearchTree(2);
auto velement = Find2dElement(mesh, p, vlam, nullopt, searchtree ? mesh.GetSurfaceElementSearchTree() : nullptr, allowindex);
if(!velement.IsValid())
return 0; return 0;
vlam[2] = 1.-vlam[0] - vlam[1]; vlam[2] = 1.-vlam[0] - vlam[1];
@ -5257,6 +5266,8 @@ namespace netgen
void Mesh :: BuildElementSearchTree (int dim) void Mesh :: BuildElementSearchTree (int dim)
{ {
if(dim < 2)
return;
if (elementsearchtreets[dim] == GetTimeStamp()) if (elementsearchtreets[dim] == GetTimeStamp())
return; return;
@ -5273,10 +5284,10 @@ namespace netgen
Point3d pmin, pmax; Point3d pmin, pmax;
GetBox(pmin, pmax); GetBox(pmin, pmax);
Box<3> box(pmin, pmax); Box<3> box(pmin, pmax);
if (dim > 1) if (dim == 3)
elementsearchtree[dim] = make_unique<BoxTree<3>>(box); elementsearchtree_vol = make_unique<BoxTree<3, ElementIndex>>(box);
else else
elementsearchtree[dim] = nullptr; // not yet implemented elementsearchtree_surf = make_unique<BoxTree<3, SurfaceElementIndex>>(box);
if (dim == 3) if (dim == 3)
{ {
@ -5316,7 +5327,7 @@ namespace netgen
} }
} }
box.Scale(1.2); box.Scale(1.2);
elementsearchtree[dim] -> Insert (box, ei+1); elementsearchtree_vol -> Insert (box, ei);
} }
} }
else if (dim == 2) else if (dim == 2)
@ -5341,7 +5352,7 @@ namespace netgen
} }
box.Scale(1.2); box.Scale(1.2);
} }
elementsearchtree[dim] -> Insert (box, ei+1); elementsearchtree_surf -> Insert (box, ei);
} }
} }
} }
@ -5381,7 +5392,7 @@ namespace netgen
bool Mesh :: PointContainedIn2DElement(const Point3d & p, bool Mesh :: PointContainedIn2DElement(const Point3d & p,
double lami[3], double lami[3],
const int element, SurfaceElementIndex ei,
bool consider3D) const bool consider3D) const
{ {
Vec3d col1, col2, col3; Vec3d col1, col2, col3;
@ -5392,9 +5403,9 @@ namespace netgen
//SZ //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 & p1 = Point(el.PNum(1));
const Point3d & p2 = Point(el.PNum(2)); const Point3d & p2 = Point(el.PNum(2));
@ -5413,7 +5424,7 @@ namespace netgen
int i = 0; int i = 0;
while(delta > 1e-16 && i < maxits) while(delta > 1e-16 && i < maxits)
{ {
curvedelems->CalcSurfaceTransformation(lam,element-1,x,Jac); curvedelems->CalcSurfaceTransformation(lam,ei,x,Jac);
rhs = p - x; rhs = p - x;
Jac.Solve(rhs,deltalam); Jac.Solve(rhs,deltalam);
lam += deltalam; lam += deltalam;
@ -5742,7 +5753,7 @@ namespace netgen
{ {
// SurfaceElement(element).GetTets (loctets); // SurfaceElement(element).GetTets (loctets);
loctrigs.SetSize(1); loctrigs.SetSize(1);
loctrigs.Elem(1) = SurfaceElement(element); loctrigs.Elem(1) = surfelements[ei];
@ -5777,11 +5788,11 @@ namespace netgen
//(*testout) << "col1 " << col1 << " col2 " << col2 << " col3 " << col3 << " rhs " << rhs << endl; //(*testout) << "col1 " << col1 << " col2 " << col2 << " col3 " << col3 << " rhs " << rhs << endl;
//(*testout) << "sol " << sol << 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(1./3,1./3);
netgen::Point<2> lam(sol.X(), sol.Y()); 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[0] = 1-sol.X()-sol.Y();
lam[1] = sol.X(); lam[1] = sol.X();
@ -5800,7 +5811,7 @@ namespace netgen
const int maxits = 30; const int maxits = 30;
while(delta > 1e-16 && i<maxits) while(delta > 1e-16 && i<maxits)
{ {
curvedelems->CalcSurfaceTransformation(lam,element-1,x,Jac); curvedelems->CalcSurfaceTransformation(lam,ei,x,Jac);
rhs = p-x; rhs = p-x;
Jac.Solve(rhs,deltalam); Jac.Solve(rhs,deltalam);
@ -5819,7 +5830,7 @@ namespace netgen
sol.X() = lam(0); sol.X() = lam(0);
sol.Y() = lam(1); sol.Y() = lam(1);
if (SurfaceElement(element).GetType() !=TRIG6 ) if (surfelements[ei].GetType() !=TRIG6 )
{ {
sol.Z() = sol.X(); sol.Z() = sol.X();
sol.X() = sol.Y(); sol.X() = sol.Y();
@ -5851,7 +5862,7 @@ namespace netgen
bool Mesh :: PointContainedIn3DElement(const Point3d & p, bool Mesh :: PointContainedIn3DElement(const Point3d & p,
double lami[3], double lami[3],
const int element) const ElementIndex ei) const
{ {
//bool oldresult = PointContainedIn3DElementOld(p,lami,element); //bool oldresult = PointContainedIn3DElementOld(p,lami,element);
//(*testout) << "old result: " << oldresult //(*testout) << "old result: " << oldresult
@ -5862,7 +5873,7 @@ namespace netgen
const double eps = 1.e-4; const double eps = 1.e-4;
const Element & el = VolumeElement(element); const Element & el = volelements[ei];
netgen::Point<3> lam = 0.0; netgen::Point<3> lam = 0.0;
@ -5897,7 +5908,7 @@ namespace netgen
const int maxits = 30; const int maxits = 30;
while(delta > 1e-16 && i<maxits) while(delta > 1e-16 && i<maxits)
{ {
curvedelems->CalcElementTransformation(lam,element-1,x,Jac); curvedelems->CalcElementTransformation(lam,ei,x,Jac);
rhs = p-x; rhs = p-x;
Jac.Solve(rhs,deltalam); Jac.Solve(rhs,deltalam);
@ -6019,91 +6030,72 @@ namespace netgen
} }
int Mesh :: GetElementOfPoint (const netgen::Point<3> & p, ElementIndex Mesh :: GetElementOfPoint (const netgen::Point<3> & p,
double lami[3], double* lami,
bool build_searchtree, bool build_searchtree,
const int index, int index,
const bool allowindex) const bool allowindex) const
{ {
if(index != -1) if(index != -1)
{ {
NgArray<int> dummy(1); Array<int> dummy(1);
dummy[0] = index; dummy[0] = index;
return GetElementOfPoint(p,lami,&dummy,build_searchtree,allowindex); return GetElementOfPoint(p,lami,dummy,build_searchtree,allowindex);
} }
else 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, ElementIndex Mesh :: GetElementOfPoint (const netgen::Point<3> & p,
double lami[3], double* lami,
const NgArray<int> * const indices, std::optional<FlatArray<int>> indices,
bool build_searchtree, bool build_searchtree,
const bool allowindex) const 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<Mesh&>(*this).BuildElementSearchTree (2);
return Find2dElement(*this, p, lami, indices, elementsearchtree[2].get(), allowindex);
}
if (build_searchtree) if (build_searchtree)
const_cast<Mesh&>(*this).BuildElementSearchTree (3); const_cast<Mesh&>(*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, SurfaceElementIndex Mesh ::
double lami[3], GetSurfaceElementOfPoint (const netgen::Point<3> & p,
bool build_searchtree, double* lami,
const int index, bool build_searchtree,
const bool allowindex) const int index,
bool allowindex) const
{ {
if(index != -1) if(index != -1)
{ {
NgArray<int> dummy(1); Array<int> dummy(1);
dummy[0] = index; dummy[0] = index;
return GetSurfaceElementOfPoint(p,lami,&dummy,build_searchtree,allowindex); return GetSurfaceElementOfPoint(p,lami,dummy,build_searchtree,allowindex);
} }
else else
return GetSurfaceElementOfPoint(p,lami,NULL,build_searchtree,allowindex); return GetSurfaceElementOfPoint(p,lami,nullopt,build_searchtree,allowindex);
} }
SurfaceElementIndex Mesh ::
GetSurfaceElementOfPoint (const netgen::Point<3> & p,
double* lami,
int Mesh :: GetSurfaceElementOfPoint (const netgen::Point<3> & p, std::optional<FlatArray<int>> indices,
double lami[3], bool build_searchtree,
const NgArray<int> * const indices, bool allowindex) const
bool build_searchtree,
const bool allowindex) const
{ {
if (dimension == 2) if (build_searchtree)
return Find1dElement(*this, p, lami, indices, nullptr, allowindex); const_cast<Mesh&>(*this).BuildElementSearchTree(2);
else return Find2dElement(*this, p, lami, indices, elementsearchtree_surf.get(), allowindex);
{
if (build_searchtree)
const_cast<Mesh&>(*this).BuildElementSearchTree(2);
return Find2dElement(*this, p, lami, indices, elementsearchtree[2].get(), allowindex);
}
return 0;
} }
void Mesh::GetIntersectingVolEls(const Point3d& p1, const Point3d& p2, void Mesh::GetIntersectingVolEls(const Point3d& p1, const Point3d& p2,
NgArray<int> & locels) const Array<ElementIndex> & locels) const
{ {
elementsearchtree[3]->GetIntersecting (p1, p2, locels); elementsearchtree_vol->GetIntersecting (p1, p2, locels);
} }
void Mesh :: SplitIntoParts() void Mesh :: SplitIntoParts()

77
libsrc/meshing/meshclass.hpp
View File

@ -148,7 +148,8 @@ namespace netgen
int numvertices; int numvertices;
/// geometric search tree for interval intersection search /// geometric search tree for interval intersection search
unique_ptr<BoxTree<3>> elementsearchtree[4] = {nullptr,nullptr,nullptr,nullptr}; unique_ptr<BoxTree<3, ElementIndex>> elementsearchtree_vol;
unique_ptr<BoxTree<3, SurfaceElementIndex>> elementsearchtree_surf;
/// time stamp for tree /// time stamp for tree
mutable size_t elementsearchtreets[4]; mutable size_t elementsearchtreets[4];
@ -200,11 +201,11 @@ namespace netgen
DLL_HEADER bool PointContainedIn2DElement(const Point3d & p, DLL_HEADER bool PointContainedIn2DElement(const Point3d & p,
double lami[3], double lami[3],
const int element, SurfaceElementIndex element,
bool consider3D = false) const; bool consider3D = false) const;
DLL_HEADER bool PointContainedIn3DElement(const Point3d & p, DLL_HEADER bool PointContainedIn3DElement(const Point3d & p,
double lami[3], double lami[3],
const int element) const; ElementIndex element) const;
DLL_HEADER bool PointContainedIn3DElementOld(const Point3d & p, DLL_HEADER bool PointContainedIn3DElementOld(const Point3d & p,
double lami[3], double lami[3],
const int element) const; const int element) const;
@ -664,38 +665,47 @@ namespace netgen
/// build box-search tree /// build box-search tree
DLL_HEADER void BuildElementSearchTree (int dim); 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;} void SetPointSearchStartElement(const int el) const {ps_startelement = el;}
/// gives element of point, barycentric coordinates /// gives element of point, barycentric coordinates
DLL_HEADER int GetElementOfPoint (const netgen::Point<3> & p, DLL_HEADER ElementIndex
double * lami, GetElementOfPoint (const netgen::Point<3> & p,
bool build_searchtree = 0, double * lami,
const int index = -1, bool build_searchtree = false,
const bool allowindex = true) const; int index = -1,
DLL_HEADER int GetElementOfPoint (const netgen::Point<3> & p, bool allowindex = true) const;
double * lami, DLL_HEADER ElementIndex
const NgArray<int> * const indices, GetElementOfPoint (const netgen::Point<3> & p,
bool build_searchtree = 0, double * lami,
const bool allowindex = true) const; std::optional<FlatArray<int>> indices,
DLL_HEADER int GetSurfaceElementOfPoint (const netgen::Point<3> & p, bool build_searchtree = 0,
double * lami, bool allowindex = true) const;
bool build_searchtree = 0, DLL_HEADER SurfaceElementIndex
const int index = -1, GetSurfaceElementOfPoint (const netgen::Point<3> & p,
const bool allowindex = true) const; double * lami,
DLL_HEADER int GetSurfaceElementOfPoint (const netgen::Point<3> & p, bool build_searchtree = false,
double * lami, int index = -1,
const NgArray<int> * const indices, bool allowindex = true) const;
bool build_searchtree = 0, DLL_HEADER SurfaceElementIndex
const bool allowindex = true) const; GetSurfaceElementOfPoint (const netgen::Point<3> & p,
double * lami,
std::optional<FlatArray<int>> indices,
bool build_searchtree = false,
bool allowindex = true) const;
/// give list of vol elements which are int the box(p1,p2) /// give list of vol elements which are int the box(p1,p2)
void GetIntersectingVolEls(const Point3d& p1, const Point3d& p2, void GetIntersectingVolEls(const Point3d& p1, const Point3d& p2,
NgArray<int> & locels) const; Array<ElementIndex> & locels) const;
/// ///
int AddFaceDescriptor(const FaceDescriptor& fd) int AddFaceDescriptor(const FaceDescriptor& fd)
@ -1040,21 +1050,6 @@ namespace netgen
} }
DLL_HEADER void AddFacesBetweenDomains(Mesh & mesh); DLL_HEADER void AddFacesBetweenDomains(Mesh & mesh);
int Find2dElement (const Mesh& mesh,
const netgen::Point<3> & p,
double * lami,
const NgArray<int> * const indices,
BoxTree<3> * searchtree,
const bool allowindex = true);
int Find3dElement (const Mesh& mesh,
const netgen::Point<3> & p,
double * lami,
const NgArray<int> * const indices,
BoxTree<3> * searchtree,
const bool allowindex = true);
} }
#endif // NETGEN_MESHCLASS_HPP #endif // NETGEN_MESHCLASS_HPP

7
libsrc/meshing/meshfunc.cpp
View File

@ -808,17 +808,16 @@ namespace netgen
double lam[3]; double lam[3];
ElementIndex ei_start = mesh.GetElementOfPoint(p_new, lam, false, domain); 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"); PrintMessage(1, "Could not find volume element with new point");
return; return;
} }
ei_start -= 1;
if(mesh[ei_start].IsDeleted()) if(mesh[ei_start].IsDeleted())
return; return;
double max_inside = -1.; 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", // search for adjacent volume element, where the new point is "most inside",
// i.e. the minimal barycentric coordinate is maximal // i.e. the minimal barycentric coordinate is maximal
@ -828,7 +827,7 @@ namespace netgen
if(mesh[ei1].IsDeleted()) if(mesh[ei1].IsDeleted())
return; return;
if(!mesh.PointContainedIn3DElement(p_new, lam, ei1+1)) if(!mesh.PointContainedIn3DElement(p_new, lam, ei1))
continue; continue;
double inside = min(min(lam[0], lam[1]), min(lam[2], 1.0-lam[0]-lam[1])); double inside = min(min(lam[0], lam[1]), min(lam[2], 1.0-lam[0]-lam[1]));

10
libsrc/visualization/vssolution.cpp
View File

@ -4337,7 +4337,7 @@ namespace netgen
} }
double lami[3]; 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) if (elnr != -1)
{ {
@ -4858,18 +4858,18 @@ namespace netgen
if(sol.iscomplex && rcomponent != 0) if(sol.iscomplex && rcomponent != 0)
{ {
rcomponent = 2 * ((rcomponent-1)/2) + 1; 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); imag);
comp = (scalcomp-1)/2 + 1; 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); printScalValue(sol, comp, val, imag, sol.iscomplex && comp > 0);
} }
if(vecfunction!=-1 && soldata[vecfunction]->draw_volume) if(vecfunction!=-1 && soldata[vecfunction]->draw_volume)
{ {
auto & sol = *soldata[vecfunction]; auto & sol = *soldata[vecfunction];
ArrayMem<double, 10> values(sol.components); ArrayMem<double, 10> 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); printVecValue(sol, values);
} }
return; return;
@ -4880,7 +4880,7 @@ namespace netgen
double lami[3] = {0.0, 0.0, 0.0}; 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) // Check if unprojected Point is close to surface element (eps of 1e-3 due to z-Buffer accuracy)
bool found_2del = false; 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 // Found it, use coordinates of point projected to surface element
mesh->GetCurvedElements().CalcSurfaceTransformation({1.0-lami[0]-lami[1], lami[0]}, selelement-1, p); mesh->GetCurvedElements().CalcSurfaceTransformation({1.0-lami[0]-lami[1], lami[0]}, selelement-1, p);