ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-02-11 12:53:08 +05:00
Code Issues Packages Projects Releases Wiki Activity

improve3 with consistent PointIndex

Browse Source
This commit is contained in:
Joachim Schoeberl 2024-12-28 13:01:20 +01:00
parent 10a56a9e86
commit 00edc92c00
5 changed files with 61 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
libsrc/meshing/improve3.cpp
View File

@ -100,7 +100,7 @@ static double SplitElementBadness (const Mesh::T_POINTS & points, const MeshingP
{ {
double badness = 0; double badness = 0;
auto np = old.GetNP(); auto np = old.GetNP();
PointIndex dummy{-1}; PointIndex dummy{PointIndex::INVALID};
if(np == 4) if(np == 4)
{ {
// Split tet into two tets // Split tet into two tets
@ -451,7 +451,7 @@ void MeshOptimize3d :: CombineImprove ()
double MeshOptimize3d :: SplitImproveEdge (Table<ElementIndex,PointIndex> & elementsonnode, NgArray<INDEX_3> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only) double MeshOptimize3d :: SplitImproveEdge (Table<ElementIndex,PointIndex> & elementsonnode, NgArray<PointIndices<3>> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only)
{ {
double d_badness = 0.0; double d_badness = 0.0;
// int cnt = 0; // int cnt = 0;
@ -512,11 +512,11 @@ double MeshOptimize3d :: SplitImproveEdge (Table<ElementIndex,PointIndex> & elem
for (int l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
if (el[l] == pi1 || el[l] == pi2) if (el[l] == pi1 || el[l] == pi2)
{ {
INDEX_3 i3; PointIndices<3> i3;
Element2d face(TRIG); Element2d face(TRIG);
el.GetFace (l+1, face); el.GetFace (l+1, face);
for (int kk = 1; kk <= 3; kk++) for (int kk = 0; kk < 3; kk++)
i3.I(kk) = face.PNum(kk); i3[kk] = face[kk];
locfaces.Append (i3); locfaces.Append (i3);
} }
} }
@ -536,7 +536,7 @@ double MeshOptimize3d :: SplitImproveEdge (Table<ElementIndex,PointIndex> & elem
{ {
int pok = pf.Func (px) < 1e10; int pok = pf.Func (px) < 1e10;
if (!pok) if (!pok)
pok = FindInnerPoint (mesh.Points(), locfaces, pnew); pok = FindInnerPoint (mesh.Points(), locfaces, pnew);
if(pok) if(pok)
{ {
@ -647,7 +647,7 @@ void MeshOptimize3d :: SplitImprove ()
tsearch.Start(); tsearch.Start();
ParallelForRange(Range(edges), [&] (auto myrange) ParallelForRange(Range(edges), [&] (auto myrange)
{ {
NgArray<INDEX_3> locfaces; NgArray<PointIndices<3>> locfaces;
for(auto i : myrange) for(auto i : myrange)
{ {
@ -671,7 +671,7 @@ void MeshOptimize3d :: SplitImprove ()
// Apply actual optimizations // Apply actual optimizations
topt.Start(); topt.Start();
int cnt = 0; int cnt = 0;
NgArray<INDEX_3> locfaces; NgArray<PointIndices<3>> locfaces;
for(auto [d_badness, ei] : edges_with_improvement) for(auto [d_badness, ei] : edges_with_improvement)
{ {
auto [p0,p1] = edges[ei]; auto [p0,p1] = edges[ei];
@ -871,7 +871,7 @@ double MeshOptimize3d :: SwapImproveEdge (
if ((goal == OPT_CONFORM) && NotTooBad(bad1, bad2)) if ((goal == OPT_CONFORM) && NotTooBad(bad1, bad2))
{ {
INDEX_3 face(pi3, pi4, pi5); PointIndices<3> face(pi3, pi4, pi5);
face.Sort(); face.Sort();
if (faces.Used(face)) if (faces.Used(face))
{ {
@ -1237,9 +1237,9 @@ double MeshOptimize3d :: SwapImproveEdge (
for (int k = l+1; k <= nsuround + l - 2; k++) for (int k = l+1; k <= nsuround + l - 2; k++)
{ {
INDEX_3 hi3(suroundpts[l], PointIndices<3> hi3(suroundpts[l],
suroundpts[k % nsuround], suroundpts[k % nsuround],
suroundpts[(k+1) % nsuround]); suroundpts[(k+1) % nsuround]);
hi3.Sort(); hi3.Sort();
if (faces.Used(hi3)) if (faces.Used(hi3))
{ {
@ -1336,7 +1336,7 @@ void MeshOptimize3d :: SwapImprove (const BitArray * working_elements)
// int ne = mesh.GetNE(); // int ne = mesh.GetNE();
mesh.BuildBoundaryEdges(false); mesh.BuildBoundaryEdges(false);
BitArray free_points(mesh.GetNP()+PointIndex::BASE); TBitArray<PointIndex> free_points(mesh.GetNP());
free_points.Clear(); free_points.Clear();
ParallelForRange(mesh.VolumeElements().Range(), [&] (auto myrange) ParallelForRange(mesh.VolumeElements().Range(), [&] (auto myrange)
@ -1372,7 +1372,7 @@ void MeshOptimize3d :: SwapImprove (const BitArray * working_elements)
for (int i = 1; i <= mesh.GetNOpenElements(); i++) for (int i = 1; i <= mesh.GetNOpenElements(); i++)
{ {
const Element2d & hel = mesh.OpenElement(i); const Element2d & hel = mesh.OpenElement(i);
INDEX_3 face(hel[0], hel[1], hel[2]); PointIndices<3> face(hel[0], hel[1], hel[2]);
face.Sort(); face.Sort();
faces.Set (face, i); faces.Set (face, i);
} }
@ -1439,7 +1439,7 @@ void MeshOptimize3d :: SwapImprove (const BitArray * working_elements)
{ {
Element2d sel; Element2d sel;
el.GetFace(i, sel); el.GetFace(i, sel);
INDEX_3 face(sel[0], sel[1], sel[2]); PointIndices<3> face(sel[0], sel[1], sel[2]);
face.Sort(); face.Sort();
if(faces.Used(face)) if(faces.Used(face))
open_els[faces.Get(face)-1].Delete(); open_els[faces.Get(face)-1].Delete();
@ -1502,9 +1502,9 @@ void MeshOptimize3d :: SwapImproveSurface (
// contains at least all elements at node // contains at least all elements at node
TABLE<ElementIndex,PointIndex::BASE> elementsonnode(np); DynamicTable<ElementIndex,PointIndex> elementsonnode(np);
TABLE<SurfaceElementIndex,PointIndex::BASE> surfaceelementsonnode(np); DynamicTable<SurfaceElementIndex,PointIndex> surfaceelementsonnode(np);
TABLE<int,PointIndex::BASE> surfaceindicesonnode(np); DynamicTable<int,PointIndex> surfaceindicesonnode(np);
NgArray<ElementIndex> hasbothpoints; NgArray<ElementIndex> hasbothpoints;
NgArray<ElementIndex> hasbothpointsother; NgArray<ElementIndex> hasbothpointsother;
@ -1607,14 +1607,14 @@ void MeshOptimize3d :: SwapImproveSurface (
othermattype = -1; othermattype = -1;
INDEX_2 i2 (pi1, pi2); PointIndices<2> i2 (pi1, pi2);
i2.Sort(); i2.Sort();
if (edgeused.Used(i2)) continue; if (edgeused.Used(i2)) continue;
edgeused.Set (i2, 1); edgeused.Set (i2, 1);
if(periodic) if(periodic)
{ {
i2.I1() = pi1other; i2[0] = pi1other;
i2.I2() = pi2other; i2[1] = pi2other;
i2.Sort(); i2.Sort();
edgeused.Set(i2,1); edgeused.Set(i2,1);
} }
@ -1770,7 +1770,7 @@ void MeshOptimize3d :: SwapImproveSurface (
//(*testout) << "sel1 " << sel1 << " sel2 " << sel2 << " el " << mesh[sel1] << " resp. " << mesh[sel2] << endl; //(*testout) << "sel1 " << sel1 << " sel2 " << sel2 << " el " << mesh[sel1] << " resp. " << mesh[sel2] << endl;
PointIndex sp1(0), sp2(0); PointIndex sp1(PointIndex::INVALID), sp2(PointIndex::INVALID);
PointIndex sp1other, sp2other; PointIndex sp1other, sp2other;
for(int l=0; l<mesh[sel1].GetNP(); l++) for(int l=0; l<mesh[sel1].GetNP(); l++)
if(mesh[sel1][l] != pi1 && mesh[sel1][l] != pi2) if(mesh[sel1][l] != pi1 && mesh[sel1][l] != pi2)
@ -2319,8 +2319,8 @@ void MeshOptimize3d :: SwapImproveSurface (
*/ */
double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face, double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
Table<ElementIndex, PointIndex> & elementsonnode, Table<ElementIndex, PointIndex> & elementsonnode,
TABLE<SurfaceElementIndex, PointIndex::BASE> & belementsonnode, bool check_only ) DynamicTable<SurfaceElementIndex, PointIndex> & belementsonnode, bool check_only )
{ {
PointIndex pi1, pi2, pi3, pi4, pi5; PointIndex pi1, pi2, pi3, pi4, pi5;
Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET); Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET);
@ -2509,7 +2509,7 @@ void MeshOptimize3d :: SwapImprove2 ()
int nse = mesh.GetNSE(); int nse = mesh.GetNSE();
// contains at least all elements at node // contains at least all elements at node
TABLE<SurfaceElementIndex, PointIndex::BASE> belementsonnode(np); DynamicTable<SurfaceElementIndex, PointIndex> belementsonnode(np);
PrintMessage (3, "SwapImprove2 "); PrintMessage (3, "SwapImprove2 ");
(*testout) << "\n" << "Start SwapImprove2" << "\n"; (*testout) << "\n" << "Start SwapImprove2" << "\n";

8
libsrc/meshing/improve3.hpp
View File

@ -35,7 +35,7 @@ public:
void CombineImprove (); void CombineImprove ();
void SplitImprove (); void SplitImprove ();
double SplitImproveEdge (Table<ElementIndex,PointIndex> & elementsonnode, NgArray<INDEX_3> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only=false); double SplitImproveEdge (Table<ElementIndex,PointIndex> & elementsonnode, NgArray<PointIndices<3>> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only=false);
void SplitImprove2 (); void SplitImprove2 ();
double SplitImprove2Element (ElementIndex ei, const Table<ElementIndex, PointIndex> & elements_of_point, bool check_only); double SplitImprove2Element (ElementIndex ei, const Table<ElementIndex, PointIndex> & elements_of_point, bool check_only);
@ -46,7 +46,7 @@ public:
void SwapImproveSurface (const BitArray * working_elements = NULL, void SwapImproveSurface (const BitArray * working_elements = NULL,
const NgArray< idmap_type* > * idmaps = NULL); const NgArray< idmap_type* > * idmaps = NULL);
void SwapImprove2 (); void SwapImprove2 ();
double SwapImprove2 (ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, TABLE<SurfaceElementIndex, PointIndex::BASE> & belementsonnode, bool check_only=false ); double SwapImprove2 (ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, DynamicTable<SurfaceElementIndex, PointIndex> & belementsonnode, bool check_only=false );
void ImproveMesh() { mesh.ImproveMesh(mp, goal); } void ImproveMesh() { mesh.ImproveMesh(mp, goal); }
@ -108,12 +108,12 @@ public:
class PointFunction1 : public MinFunction class PointFunction1 : public MinFunction
{ {
Mesh::T_POINTS & points; Mesh::T_POINTS & points;
const NgArray<INDEX_3> & faces; const NgArray<PointIndices<3>> & faces;
const MeshingParameters & mp; const MeshingParameters & mp;
double h; double h;
public: public:
PointFunction1 (Mesh::T_POINTS & apoints, PointFunction1 (Mesh::T_POINTS & apoints,
const NgArray<INDEX_3> & afaces, const NgArray<PointIndices<3>> & afaces,
const MeshingParameters & amp, const MeshingParameters & amp,
double ah); double ah);

40
libsrc/meshing/meshclass.cpp
View File

@ -1052,7 +1052,7 @@ namespace netgen
// for (PointIndex pi = points.Begin(); pi < points.End(); pi++) // for (PointIndex pi = points.Begin(); pi < points.End(); pi++)
for (PointIndex pi : points.Range()) for (PointIndex pi : points.Range())
if ((*this)[pi].Singularity()>=1.) if ((*this)[pi].Singularity()>=1.)
outfile << int(pi) << "\t" << (*this)[pi].Singularity() << endl; outfile << pi << "\t" << (*this)[pi].Singularity() << endl;
} }
cnt_sing = 0; cnt_sing = 0;
@ -1341,7 +1341,7 @@ namespace netgen
el.SetNP(nep); el.SetNP(nep);
el.SetCurved (nep != 4); el.SetCurved (nep != 4);
for (int j = 0; j < nep; j++) for (int j = 0; j < nep; j++)
infile >> (int&)(el[j]); infile >> el[j];
if (inverttets) if (inverttets)
el.Invert(); el.Invert();
@ -2113,7 +2113,7 @@ namespace netgen
for (int j = 0; j < nep; j++) for (int j = 0; j < nep; j++)
{ {
infile >> (int&)(el[j]); infile >> el[j];
el[j] = el[j]+oldnp; el[j] = el[j]+oldnp;
} }
@ -2182,7 +2182,7 @@ namespace netgen
for (ElementIndex ei = 0; ei < volelements.Size(); ei++) for (ElementIndex ei = 0; ei < volelements.Size(); ei++)
{ {
for (int j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
if ( (*this)[ei][j] <= PointIndex::BASE-1) if ( !(*this)[ei][j].IsValid())
{ {
(*testout) << "El " << ei << " has 0 nodes: "; (*testout) << "El " << ei << " has 0 nodes: ";
for (int k = 0; k < 4; k++) for (int k = 0; k < 4; k++)
@ -2223,9 +2223,9 @@ namespace netgen
if (sel.GetNP() <= 4) if (sel.GetNP() <= 4)
for (int j = 0; j < sel.GetNP(); j++) for (int j = 0; j < sel.GetNP(); j++)
{ {
INDEX_2 i2; PointIndices<2> i2;
i2.I1() = sel.PNumMod(j+1); i2[0] = sel.PNumMod(j+1);
i2.I2() = sel.PNumMod(j+2); i2[1] = sel.PNumMod(j+2);
i2.Sort(); i2.Sort();
boundaryedges->Set (i2, 1); boundaryedges->Set (i2, 1);
} }
@ -2233,9 +2233,9 @@ namespace netgen
{ {
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
{ {
INDEX_2 i2; PointIndices<2> i2;
i2.I1() = sel[j]; i2[0] = sel[j];
i2.I2() = sel[(j+1)%3]; i2[1] = sel[(j+1)%3];
i2.Sort(); i2.Sort();
boundaryedges->Set (i2, 1); boundaryedges->Set (i2, 1);
} }
@ -2298,7 +2298,7 @@ namespace netgen
static Timer tn2se("Mesh::CalcSurfacesOfNode - surf on node"); static Timer tn2se("Mesh::CalcSurfacesOfNode - surf on node");
static Timer tht("Mesh::CalcSurfacesOfNode - surfelementht"); static Timer tht("Mesh::CalcSurfacesOfNode - surfelementht");
// surfacesonnode.SetSize (GetNP()); // surfacesonnode.SetSize (GetNP());
TABLE<int,PointIndex::BASE> surfacesonnode(GetNP()); DynamicTable<int,PointIndex> surfacesonnode(GetNP());
// delete boundaryedges; // delete boundaryedges;
// boundaryedges = NULL; // boundaryedges = NULL;
@ -2379,10 +2379,10 @@ namespace netgen
const Element2d & sel = surfelements[sei]; const Element2d & sel = surfelements[sei];
if (sel.IsDeleted()) continue; if (sel.IsDeleted()) continue;
INDEX_3 i3; PointIndices<3> i3;
i3.I1() = sel.PNum(1); i3[0] = sel.PNum(1);
i3.I2() = sel.PNum(2); i3[1] = sel.PNum(2);
i3.I3() = sel.PNum(3); i3[2] = sel.PNum(3);
i3.Sort(); i3.Sort();
surfelementht -> Set (i3, sei); // war das wichtig ??? sel.GetIndex()); surfelementht -> Set (i3, sei); // war das wichtig ??? sel.GetIndex());
} }
@ -2479,7 +2479,7 @@ namespace netgen
for (int i = 0; i < GetNSeg(); i++) for (int i = 0; i < GetNSeg(); i++)
{ {
const Segment & seg = segments[i]; const Segment & seg = segments[i];
INDEX_2 i2(seg[0], seg[1]); PointIndices<2> i2(seg[0], seg[1]);
i2.Sort(); i2.Sort();
//boundaryedges -> Set (i2, 2); //boundaryedges -> Set (i2, 2);
@ -2817,7 +2817,7 @@ namespace netgen
hel.NormalizeNumbering(); hel.NormalizeNumbering();
if (hel.PNum(1) == pi) if (hel.PNum(1) == pi)
{ {
INDEX_3 i3(hel[0], hel[1], hel[2]); PointIndices<3> i3(hel[0], hel[1], hel[2]);
tval i2; tval i2;
i2.index = GetFaceDescriptor(ind).DomainIn(); i2.index = GetFaceDescriptor(ind).DomainIn();
i2.p4 = (hel.GetNP() == 3) i2.p4 = (hel.GetNP() == 3)
@ -2833,7 +2833,7 @@ namespace netgen
hel.NormalizeNumbering(); hel.NormalizeNumbering();
if (hel.PNum(1) == pi) if (hel.PNum(1) == pi)
{ {
INDEX_3 i3(hel[0], hel[1], hel[2]); PointIndices<3> i3(hel[0], hel[1], hel[2]);
tval i2; tval i2;
i2.index = GetFaceDescriptor(ind).DomainOut(); i2.index = GetFaceDescriptor(ind).DomainOut();
i2.p4 = (hel.GetNP() == 3) i2.p4 = (hel.GetNP() == 3)
@ -2860,7 +2860,7 @@ namespace netgen
if (hel[0] == pi) if (hel[0] == pi)
{ {
INDEX_3 i3(hel[0], hel[1], hel[2]); PointIndices<3> i3(hel[0], hel[1], hel[2]);
if (faceht.Used (i3)) if (faceht.Used (i3))
{ {
@ -2889,7 +2889,7 @@ namespace netgen
{ {
hel.Invert(); hel.Invert();
hel.NormalizeNumbering(); hel.NormalizeNumbering();
INDEX_3 i3(hel[0], hel[1], hel[2]); PointIndices<3> i3(hel[0], hel[1], hel[2]);
tval i2; tval i2;
i2.index = el.GetIndex(); i2.index = el.GetIndex();

14
libsrc/meshing/meshtype.hpp
View File

@ -220,11 +220,10 @@ namespace netgen
#endif #endif
} }
/*
friend constexpr netgen::PointIndex ngcore::IndexBASE<netgen::PointIndex> (); friend constexpr netgen::PointIndex ngcore::IndexBASE<netgen::PointIndex> ();
// friend istream & operator>> (istream &, PointIndex &);
// friend ostream & operator<< (ostream &, const PointIndex &);
template <int N> friend class PointIndices; template <int N> friend class PointIndices;
/*
friend PointIndex operator+ (PointIndex, int); friend PointIndex operator+ (PointIndex, int);
friend PointIndex operator+ (PointIndex, size_t); friend PointIndex operator+ (PointIndex, size_t);
friend PointIndex operator+ (int, PointIndex); friend PointIndex operator+ (int, PointIndex);
@ -237,6 +236,7 @@ namespace netgen
friend bool operator== (PointIndex a, PointIndex b); friend bool operator== (PointIndex a, PointIndex b);
friend bool operator!= (PointIndex a, PointIndex b); friend bool operator!= (PointIndex a, PointIndex b);
*/ */
public: public:
constexpr PointIndex (t_invalid inv) : i(long(PointIndex::BASE)-1) { ; } constexpr PointIndex (t_invalid inv) : i(long(PointIndex::BASE)-1) { ; }
// PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; } // PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; }
@ -260,6 +260,7 @@ namespace netgen
void DoArchive (Archive & ar) { ar & i; } void DoArchive (Archive & ar) { ar & i; }
}; };
/* /*
inline PointIndex operator+ (PointIndex pi, int i) { return PointIndex(pi.i+i); } inline PointIndex operator+ (PointIndex pi, int i) { return PointIndex(pi.i+i); }
inline PointIndex operator+ (PointIndex pi, size_t i) { return PointIndex(pi.i+i); } inline PointIndex operator+ (PointIndex pi, size_t i) { return PointIndex(pi.i+i); }
@ -367,12 +368,19 @@ namespace netgen
{ Sort(); } { Sort(); }
}; };
constexpr inline size_t HashValue2 (const PointIndex & ind, size_t mask)
{ return (ind-IndexBASE<PointIndex>()) & mask; }
} }
namespace ngcore namespace ngcore
{ {
template <typename T> constexpr inline T InvalidHash(); template <typename T> constexpr inline T InvalidHash();
template <>
constexpr inline netgen::PointIndex InvalidHash<netgen::PointIndex> ()
{ return netgen::PointIndex::INVALID; }
template <> template <>
constexpr inline netgen::PointIndices<2> InvalidHash<netgen::PointIndices<2>> () constexpr inline netgen::PointIndices<2> InvalidHash<netgen::PointIndices<2>> ()

2
libsrc/meshing/smoothing3.cpp
View File

@ -91,7 +91,7 @@ namespace netgen
} }
PointFunction1 :: PointFunction1 (Mesh::T_POINTS & apoints, PointFunction1 :: PointFunction1 (Mesh::T_POINTS & apoints,
const NgArray<INDEX_3> & afaces, const NgArray<PointIndices<3>> & afaces,
const MeshingParameters & amp, const MeshingParameters & amp,
double ah) double ah)
: points(apoints), faces(afaces), mp(amp) : points(apoints), faces(afaces), mp(amp)