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

Merge branch 'cleanup_searchtree_elementindex' into 'master'

Browse Source 
move all searchtrees to use elementindex

See merge request ngsolve/netgen!702
This commit is contained in:
Lackner, Christopher 2025-03-13 19:55:29 +01:00
commit 8944322e60
8 changed files with 208 additions and 211 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]);
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
{

27
libsrc/interface/nginterface_v2.cpp
View File

@ -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<int> ind(numind, indices);
return Find2dElement(*mesh, pp, lami, &ind, mesh->GetElementSearchTree(2));
FlatArray<int> 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<int> ind(numind, indices);
return Find3dElement(*mesh, pp, lami, &ind, mesh->GetElementSearchTree(3));
FlatArray<int> ind(numind, indices);
return mesh->GetElementOfPoint(pp, lami, ind, build_searchtree);
}
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;
NgArray<INDEX_3> surfaceelp;
NgArray<int> surfaceeli;
NgArray<int> locels;
Array<ElementIndex> 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;
}

30
libsrc/meshing/fieldlines.cpp
View File

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

251
libsrc/meshing/meshclass.cpp
View File

@ -12,15 +12,15 @@
namespace netgen
{
int Find3dElement (const Mesh& mesh,
const netgen::Point<3> & p,
double * lami,
const NgArray<int> * const indices,
BoxTree<3> * searchtree,
const bool allowindex)
ElementIndex Find3dElement (const Mesh& mesh,
const netgen::Point<3> & p,
double * lami,
optional<FlatArray<int>> indices,
BoxTree<3, ElementIndex> * searchtree,
const bool allowindex)
{
int ne = 0;
NgArray<int> locels;
Array<ElementIndex> locels;
if (searchtree)
{
searchtree->GetIntersecting (p, p, locels);
@ -29,96 +29,90 @@ namespace netgen
else
ne = mesh.GetNE();
for (int i = 1; i <= ne; i++)
for (auto i : Range(ne))
{
int ii;
ElementIndex ei;
if (searchtree)
ii = locels.Get(i);
ei = locels[i];
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(mesh.PointContainedIn3DElement(p,lami,ii))
return ii;
if(mesh.PointContainedIn3DElement(p,lami,ei))
return ei;
}
// Not found, try uncurved variant:
for (int i = 1; i <= ne; i++)
for (auto i : Range(ne))
{
int ii;
ElementIndex ei;
if (searchtree)
ii = locels.Get(i);
ei = locels[i];
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(mesh.PointContainedIn3DElementOld(p,lami,ii))
if(mesh.PointContainedIn3DElementOld(p,lami,ei))
{
(*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,
const netgen::Point<3> & p,
double * lami,
const NgArray<int> * const indices,
BoxTree<3> * searchtree,
const bool allowindex)
SurfaceElementIndex
Find2dElement (const Mesh& mesh,
const netgen::Point<3> & p,
double * lami,
std::optional<FlatArray<int>> 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&>(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<int> faces;
// topology.GetElementFaces(velement,faces);
auto faces = Array<int> (topology.GetFaces(ElementIndex(velement-1)));
//(*testout) << "faces " << faces << endl;
for(int i=0; i<faces.Size(); i++)
faces[i] = topology.GetFace2SurfaceElement(faces[i])+1;
//(*testout) << "surfel " << faces << endl;
const auto & fnrs = topology.GetFaces(velement);
auto faces = ArrayMem<SurfaceElementIndex,4>();
for(auto face : fnrs)
faces.Append(topology.GetFace2SurfaceElement(face));
for(int i=0; i<faces.Size(); i++)
{
if(faces[i] == 0)
if(!faces[i].IsValid())
continue;
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;
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<int> locels;
// TODO: build search tree for surface elements
Array<SurfaceElementIndex> 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<int> * const indices,
BoxTree<3> * searchtree,
const bool allowindex = true)
SegmentIndex Find1dElement (const Mesh& mesh,
const netgen::Point<3> & p,
double * lami,
std::optional<FlatArray<int>> 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&>(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<BoxTree<3>>(box);
if (dim == 3)
elementsearchtree_vol = make_unique<BoxTree<3, ElementIndex>>(box);
else
elementsearchtree[dim] = nullptr; // not yet implemented
elementsearchtree_surf = make_unique<BoxTree<3, SurfaceElementIndex>>(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 && i<maxits)
{
curvedelems->CalcSurfaceTransformation(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 && i<maxits)
{
curvedelems->CalcElementTransformation(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<int> dummy(1);
Array<int> 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<int> * const indices,
bool build_searchtree,
const bool allowindex) const
ElementIndex Mesh :: GetElementOfPoint (const netgen::Point<3> & p,
double* lami,
std::optional<FlatArray<int>> 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<Mesh&>(*this).BuildElementSearchTree (2);
return Find2dElement(*this, p, lami, indices, elementsearchtree[2].get(), allowindex);
}
if (build_searchtree)
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,
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<int> dummy(1);
Array<int> 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<int> * const indices,
bool build_searchtree,
const bool allowindex) const
SurfaceElementIndex Mesh ::
GetSurfaceElementOfPoint (const netgen::Point<3> & p,
double* lami,
std::optional<FlatArray<int>> indices,
bool build_searchtree,
bool allowindex) const
{
if (dimension == 2)
return Find1dElement(*this, p, lami, indices, nullptr, allowindex);
else
{
if (build_searchtree)
const_cast<Mesh&>(*this).BuildElementSearchTree(2);
return Find2dElement(*this, p, lami, indices, elementsearchtree[2].get(), allowindex);
}
return 0;
if (build_searchtree)
const_cast<Mesh&>(*this).BuildElementSearchTree(2);
return Find2dElement(*this, p, lami, indices, elementsearchtree_surf.get(), allowindex);
}
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()

77
libsrc/meshing/meshclass.hpp
View File

@ -148,7 +148,8 @@ namespace netgen
int numvertices;
/// 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
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<int> * 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<int> * 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<FlatArray<int>> 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<FlatArray<int>> 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<int> & locels) const;
Array<ElementIndex> & 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<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

7
libsrc/meshing/meshfunc.cpp
View File

@ -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]));

10
libsrc/visualization/vssolution.cpp
View File

@ -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<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);
}
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);