code cleanup

This commit is contained in:
Joachim Schoeberl 2024-12-15 18:00:50 +01:00
parent b560719a47
commit 8f73a00d2d
17 changed files with 164 additions and 141 deletions

View File

@ -182,9 +182,13 @@ void AdFront3 :: DeleteFace (INDEX fi)
{ {
nff--; nff--;
/*
for (int i = 1; i <= faces.Get(fi).Face().GetNP(); i++) for (int i = 1; i <= faces.Get(fi).Face().GetNP(); i++)
{ {
PointIndex pi = faces.Get(fi).Face().PNum(i); PointIndex pi = faces.Get(fi).Face().PNum(i);
*/
for (PointIndex pi : faces.Get(fi).Face().PNums())
{
points[pi].RemoveFace(); points[pi].RemoveFace();
if (!points[pi].Valid()) if (!points[pi].Valid())
delpointl.Append (pi); delpointl.Append (pi);
@ -397,16 +401,10 @@ void AdFront3 :: RebuildInternalTables ()
int negvol = 0; bool negvol = false;
/*
for (int i = PointIndex::BASE;
i < clvol.Size()+PointIndex::BASE; i++)
*/
for (auto i : clvol.Range()) for (auto i : clvol.Range())
{ if (clvol[i] < 0)
if (clvol[i] < 0) negvol = true;
negvol = 1;
}
if (negvol) if (negvol)
{ {
@ -427,8 +425,6 @@ void AdFront3 :: RebuildInternalTables ()
int AdFront3 :: SelectBaseElement () int AdFront3 :: SelectBaseElement ()
{ {
int i, hi, fstind;
/* /*
static int minval = -1; static int minval = -1;
static int lasti = 0; static int lasti = 0;
@ -453,12 +449,12 @@ int AdFront3 :: SelectBaseElement ()
} }
*/ */
fstind = 0; int fstind = 0;
for (i = lasti+1; i <= faces.Size() && !fstind; i++) for (int i = lasti+1; i <= faces.Size() && !fstind; i++)
if (faces.Elem(i).Valid()) if (faces.Elem(i).Valid())
{ {
hi = faces.Get(i).QualClass() + int hi = faces.Get(i).QualClass() +
points[faces.Get(i).Face().PNum(1)].FrontNr() + points[faces.Get(i).Face().PNum(1)].FrontNr() +
points[faces.Get(i).Face().PNum(2)].FrontNr() + points[faces.Get(i).Face().PNum(2)].FrontNr() +
points[faces.Get(i).Face().PNum(3)].FrontNr(); points[faces.Get(i).Face().PNum(3)].FrontNr();
@ -474,10 +470,10 @@ int AdFront3 :: SelectBaseElement ()
if (!fstind) if (!fstind)
{ {
minval = INT_MAX; minval = INT_MAX;
for (i = 1; i <= faces.Size(); i++) for (int i = 1; i <= faces.Size(); i++)
if (faces.Elem(i).Valid()) if (faces.Elem(i).Valid())
{ {
hi = faces.Get(i).QualClass() + int hi = faces.Get(i).QualClass() +
points[faces.Get(i).Face().PNum(1)].FrontNr() + points[faces.Get(i).Face().PNum(1)].FrontNr() +
points[faces.Get(i).Face().PNum(2)].FrontNr() + points[faces.Get(i).Face().PNum(2)].FrontNr() +
points[faces.Get(i).Face().PNum(3)].FrontNr(); points[faces.Get(i).Face().PNum(3)].FrontNr();
@ -499,9 +495,9 @@ int AdFront3 :: SelectBaseElement ()
int AdFront3 :: GetLocals (int fstind, int AdFront3 :: GetLocals (int fstind,
Array<Point3d, PointIndex> & locpoints, Array<Point3d, PointIndex> & locpoints,
NgArray<MiniElement2d> & locfaces, // local index Array<MiniElement2d> & locfaces, // local index
Array<PointIndex, PointIndex> & pindex, Array<PointIndex, PointIndex> & pindex,
NgArray<INDEX> & findex, Array<INDEX> & findex,
INDEX_2_HASHTABLE<int> & getconnectedpairs, INDEX_2_HASHTABLE<int> & getconnectedpairs,
float xh, float xh,
float relh, float relh,
@ -518,7 +514,7 @@ int AdFront3 :: GetLocals (int fstind,
hashcreated=1; hashcreated=1;
} }
INDEX i, j; INDEX i;
PointIndex pstind; PointIndex pstind;
Point3d midp, p0; Point3d midp, p0;
@ -599,28 +595,37 @@ int AdFront3 :: GetLocals (int fstind,
invpindex.SetSize (points.Size()); invpindex.SetSize (points.Size());
/*
for (i = 1; i <= locfaces.Size(); i++) for (i = 1; i <= locfaces.Size(); i++)
for (j = 1; j <= locfaces.Get(i).GetNP(); j++) for (j = 1; j <= locfaces.Get(i).GetNP(); j++)
{ {
PointIndex pi = locfaces.Get(i).PNum(j); PointIndex pi = locfaces.Get(i).PNum(j);
invpindex[pi] = PointIndex::INVALID; invpindex[pi] = PointIndex::INVALID;
} }
*/
for (auto & f : locfaces)
for (int j = 1; j <= f.GetNP(); j++)
{
PointIndex pi = f.PNum(j);
invpindex[pi] = PointIndex::INVALID;
}
for (i = 1; i <= locfaces.Size(); i++) // for (i = 1; i <= locfaces.Size(); i++)
for (auto & f : locfaces)
{ {
for (j = 1; j <= locfaces.Get(i).GetNP(); j++) // for (j = 1; j <= locfaces.Get(i).GetNP(); j++)
for (int j = 1; j <= f.GetNP(); j++)
{ {
PointIndex pi = locfaces.Get(i).PNum(j); // PointIndex pi = locfaces.Get(i).PNum(j);
PointIndex pi = f.PNum(j);
if (!invpindex[pi].IsValid()) if (!invpindex[pi].IsValid())
{ {
pindex.Append (pi); pindex.Append (pi);
locpoints.Append (points[pi].P()); locpoints.Append (points[pi].P());
invpindex[pi] = pindex.Size()-1+IndexBASE<PointIndex>(); invpindex[pi] = pindex.Size()-1+IndexBASE<PointIndex>();
} }
// locfaces.Elem(i).PNum(j) = locpoints.Append (points[pi].P()); // locfaces.Elem(i).PNum(j) = invpindex[pi];
// } f.PNum(j) = invpindex[pi];
// else
locfaces.Elem(i).PNum(j) = invpindex[pi];
} }
} }
@ -675,10 +680,10 @@ int AdFront3 :: GetLocals (int fstind,
// returns all points connected with fi // returns all points connected with fi
void AdFront3 :: GetGroup (int fi, void AdFront3 :: GetGroup (int fi,
NgArray<MeshPoint, PointIndex::BASE> & grouppoints, Array<MeshPoint, PointIndex> & grouppoints,
NgArray<MiniElement2d> & groupelements, Array<MiniElement2d> & groupelements,
Array<PointIndex, PointIndex> & pindex, Array<PointIndex, PointIndex> & pindex,
NgArray<INDEX> & findex) Array<INDEX> & findex)
{ {
// static NgArray<char> pingroup; // static NgArray<char> pingroup;
int changed; int changed;
@ -796,8 +801,8 @@ void AdFront3 :: SetStartFront (int /* baseelnp */)
*/ */
} }
bool AdFront3 :: PointInsideGroup(const Array<PointIndex, PointIndex> &grouppindex, bool AdFront3 :: PointInsideGroup(const Array<PointIndex, PointIndex> &grouppindex,
const NgArray<MiniElement2d> &groupfaces) const const Array<MiniElement2d> &groupfaces) const
{ {
for(auto pi : Range(points)) for(auto pi : Range(points))
{ {

View File

@ -95,7 +95,8 @@ public:
const PointIndex PNum (int i) const { return pnum[i-1]; } const PointIndex PNum (int i) const { return pnum[i-1]; }
PointIndex & PNum (int i) { return pnum[i-1]; } PointIndex & PNum (int i) { return pnum[i-1]; }
const PointIndex PNumMod (int i) const { return pnum[(i-1)%np]; } const PointIndex PNumMod (int i) const { return pnum[(i-1)%np]; }
auto PNums() const { return NgFlatArray<const PointIndex> (np, &pnum[0]); } auto PNums() { return FlatArray<PointIndex> (np, &pnum[0]); }
auto PNums() const { return FlatArray<const PointIndex> (np, &pnum[0]); }
void Delete () { deleted = true; for (PointIndex & p : pnum) p.Invalidate(); } void Delete () { deleted = true; for (PointIndex & p : pnum) p.Invalidate(); }
bool IsDeleted () const { return deleted; } bool IsDeleted () const { return deleted; }
}; };
@ -238,9 +239,9 @@ public:
/// ///
int GetNF() const int GetNF() const
{ return nff; } { return nff; }
/// /// 1-based
const MiniElement2d & GetFace (int i) const const MiniElement2d & GetFace (int i) const
{ return faces.Get(i).Face(); } { return faces[i-1].Face(); }
const auto & Faces() const { return faces; } const auto & Faces() const { return faces; }
/// ///
void Print () const; void Print () const;
@ -265,7 +266,7 @@ public:
NgArray<int> & ifaces) const; NgArray<int> & ifaces) const;
bool PointInsideGroup(const Array<PointIndex, PointIndex> &grouppindex, bool PointInsideGroup(const Array<PointIndex, PointIndex> &grouppindex,
const NgArray<MiniElement2d>& groupfaces) const; const Array<MiniElement2d>& groupfaces) const;
/// ///
void GetFaceBoundingBox (int i, Box3d & box) const; void GetFaceBoundingBox (int i, Box3d & box) const;
@ -273,9 +274,9 @@ public:
/// ///
int GetLocals (int baseelement, int GetLocals (int baseelement,
Array<Point3d, PointIndex> & locpoints, Array<Point3d, PointIndex> & locpoints,
NgArray<MiniElement2d> & locfaces, // local index Array<MiniElement2d> & locfaces, // local index
Array<PointIndex, PointIndex> & pindex, Array<PointIndex, PointIndex> & pindex,
NgArray<INDEX> & findex, Array<INDEX> & findex,
INDEX_2_HASHTABLE<int> & connectedpairs, INDEX_2_HASHTABLE<int> & connectedpairs,
float xh, float xh,
float relh, float relh,
@ -283,10 +284,10 @@ public:
/// ///
void GetGroup (int fi, void GetGroup (int fi,
NgArray<MeshPoint, PointIndex::BASE> & grouppoints, Array<MeshPoint, PointIndex> & grouppoints,
NgArray<MiniElement2d> & groupelements, Array<MiniElement2d> & groupelements,
Array<PointIndex, PointIndex> & pindex, Array<PointIndex, PointIndex> & pindex,
NgArray<INDEX> & findex); Array<INDEX> & findex);
/// ///
void DeleteFace (INDEX fi); void DeleteFace (INDEX fi);
@ -298,11 +299,11 @@ public:
INDEX AddConnectedPair (PointIndices<2> pair); INDEX AddConnectedPair (PointIndices<2> pair);
/// ///
void IncrementClass (INDEX fi) void IncrementClass (INDEX fi)
{ faces.Elem(fi).IncrementQualClass(); } { faces[fi-1].IncrementQualClass(); }
/// ///
void ResetClass (INDEX fi) void ResetClass (INDEX fi)
{ faces.Elem(fi).ResetQualClass(); } { faces[fi-1].ResetQualClass(); }
/// ///
void SetStartFront (int baseelnp = 0); void SetStartFront (int baseelnp = 0);

View File

@ -527,12 +527,12 @@ namespace netgen
{ {
PointIndices<2> e1(el2d[i], el2d[(i+1) % el2d.GetNP()]); PointIndices<2> e1(el2d[i], el2d[(i+1) % el2d.GetNP()]);
e1.Sort(); e1.Sort();
INDEX_2 e2; PointIndices<2> e2;
for(k = 0; k < idmaps.Size(); k++) for(k = 0; k < idmaps.Size(); k++)
{ {
e2.I1() = (*idmaps[k])[e1.I1()]; e2[0] = (*idmaps[k])[e1[0]];
e2.I2() = (*idmaps[k])[e1.I2()]; e2[1] = (*idmaps[k])[e1[1]];
if(e2.I1() == 0 || e2.I2() == 0 || if(e2.I1() == 0 || e2.I2() == 0 ||
e1.I1() == e2.I1() || e1.I2() == e2.I2()) e1.I1() == e2.I1() || e1.I2() == e2.I2())
@ -985,7 +985,7 @@ namespace netgen
if(j == 0 || mi.pnums[j+mi.np] < min2) if(j == 0 || mi.pnums[j+mi.np] < min2)
min2 = mi.pnums[j+mi.np]; min2 = mi.pnums[j+mi.np];
identified = (mi.pnums[j+mi.np] != 0 && mi.pnums[j+mi.np] != mi.pnums[j]); identified = (mi.pnums[j+mi.np].IsValid() && mi.pnums[j+mi.np] != mi.pnums[j]);
} }
identified = identified && (min1 < min2); identified = identified && (min1 < min2);
@ -1929,13 +1929,14 @@ namespace netgen
ist >> size; ist >> size;
mtets.SetSize(size); mtets.SetSize(size);
constexpr auto PI0 = IndexBASE<PointIndex>();
for(int i=0; i<size; i++) for(int i=0; i<size; i++)
{ {
ist >> mtets[i]; ist >> mtets[i];
if(mtets[i].pnums[0] > mesh.GetNV() || if(mtets[i].pnums[0] >= PI0+mesh.GetNV() ||
mtets[i].pnums[1] > mesh.GetNV() || mtets[i].pnums[1] >= PI0+mesh.GetNV() ||
mtets[i].pnums[2] > mesh.GetNV() || mtets[i].pnums[2] >= PI0+mesh.GetNV() ||
mtets[i].pnums[3] > mesh.GetNV()) mtets[i].pnums[3] >= PI0+mesh.GetNV())
return false; return false;
} }

View File

@ -146,7 +146,7 @@ namespace netgen
auto & seg = mesh[segi]; auto & seg = mesh[segi];
if(seg.edgenr != edgenr+1) if(seg.edgenr != edgenr+1)
continue; continue;
PointIndex other = seg[0]+seg[1]-pi; PointIndex other = seg[0]-pi+seg[1];
if(!pts.Contains(other)) if(!pts.Contains(other))
pts.Append(other); pts.Append(other);
} }
@ -1353,7 +1353,8 @@ namespace netgen
nel[1] = p2; nel[1] = p2;
nel[2] = p3; nel[2] = p3;
nel[3] = p4; nel[3] = p4;
nel[4] = el[0] + el[1] + el[2] + el[3] - fixed[0] - moved[0] - moved[1]; // nel[4] = el[0] + el[1] + el[2] + el[3] - fixed[0] - moved[0] - moved[1];
nel[4] = el[0]-fixed[0] + el[1]-moved[0] + el[2]-moved[1] + el[3];
if(Cross(mesh[p2]-mesh[p1], mesh[p4]-mesh[p1]) * (mesh[nel[4]]-mesh[nel[1]]) > 0) if(Cross(mesh[p2]-mesh[p1], mesh[p4]-mesh[p1]) * (mesh[nel[4]]-mesh[nel[1]]) > 0)
Swap(nel[1], nel[3]); Swap(nel[1], nel[3]);
nel.SetIndex(el.GetIndex()); nel.SetIndex(el.GetIndex());
@ -1399,7 +1400,8 @@ namespace netgen
nel[1] = p2; nel[1] = p2;
nel[2] = p3; nel[2] = p3;
nel[3] = p4; nel[3] = p4;
nel[4] = el[0] + el[1] + el[2] + el[3] + el[4] - fixed[0] - fixed[1] - moved[0] - moved[1]; // nel[4] = el[0] + el[1] + el[2] + el[3] + el[4] - fixed[0] - fixed[1] - moved[0] - moved[1];
nel[4] = el[0]-fixed[0] + el[1]-fixed[1] + el[2]-moved[0] + el[3]-moved[1] + el[4];
if(Cross(mesh[p2] - mesh[p1], mesh[p4]-mesh[p1]) * (mesh[nel[4]]-mesh[nel[1]]) > 0) if(Cross(mesh[p2] - mesh[p1], mesh[p4]-mesh[p1]) * (mesh[nel[4]]-mesh[nel[1]]) > 0)
Swap(nel[1], nel[3]); Swap(nel[1], nel[3]);
nel.SetIndex(el.GetIndex()); nel.SetIndex(el.GetIndex());

View File

@ -10,6 +10,7 @@ HPREF_ELEMENT_TYPE ClassifyTet(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges
HPREF_ELEMENT_TYPE type = HP_NONE; HPREF_ELEMENT_TYPE type = HP_NONE;
int debug = 0; int debug = 0;
/*
for (int j = 0;j < 4; j++) for (int j = 0;j < 4; j++)
{ {
if (el.pnums[j] == 444) debug++; if (el.pnums[j] == 444) debug++;
@ -18,7 +19,7 @@ HPREF_ELEMENT_TYPE ClassifyTet(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges
if (el.pnums[j] == 281) debug++; if (el.pnums[j] == 281) debug++;
} }
if (debug < 4) debug = 0; if (debug < 4) debug = 0;
*/
// *testout << "new el" << endl; // *testout << "new el" << endl;

View File

@ -23,6 +23,7 @@ namespace netgen
// static int timer2 = NgProfiler::CreateTimer ("clusters2"); // static int timer2 = NgProfiler::CreateTimer ("clusters2");
// static int timer3 = NgProfiler::CreateTimer ("clusters3"); // static int timer3 = NgProfiler::CreateTimer ("clusters3");
RegionTimer reg (timer); RegionTimer reg (timer);
constexpr auto PI0 = IndexBASE<PointIndex>();
const MeshTopology & top = mesh.GetTopology(); const MeshTopology & top = mesh.GetTopology();
@ -103,7 +104,7 @@ namespace netgen
nnums.SetSize(elnv+elned+elnfa+1); nnums.SetSize(elnv+elned+elnfa+1);
for (int j = 0; j < elnv; j++) for (int j = 0; j < elnv; j++)
nnums[j] = el[j]+1-PointIndex::BASE; nnums[j] = el[j]+1-PI0;
for (int j = 0; j < elned; j++) for (int j = 0; j < elned; j++)
nnums[elnv+j] = nv+ednums[j]+1; nnums[elnv+j] = nv+ednums[j]+1;
for (int j = 0; j < elnfa; j++) for (int j = 0; j < elnfa; j++)
@ -131,7 +132,7 @@ namespace netgen
nnums.SetSize(elnv+elned+1); nnums.SetSize(elnv+elned+1);
for (int j = 1; j <= elnv; j++) for (int j = 1; j <= elnv; j++)
nnums.Elem(j) = el.PNum(j)+1-PointIndex::BASE; nnums.Elem(j) = el.PNum(j)+1-PI0;
for (int j = 1; j <= elned; j++) for (int j = 1; j <= elned; j++)
nnums.Elem(elnv+j) = nv+ednums.Elem(j); nnums.Elem(elnv+j) = nv+ednums.Elem(j);
nnums.Elem(elnv+elned+1) = fanum; nnums.Elem(elnv+elned+1) = fanum;
@ -162,7 +163,7 @@ namespace netgen
nnums.SetSize(elnv+elned+1); nnums.SetSize(elnv+elned+1);
for (int j = 0; j < elnv; j++) for (int j = 0; j < elnv; j++)
nnums[j] = el[j]+1-PointIndex::BASE; nnums[j] = el[j]+1-PI0;
for (int j = 0; j < elned; j++) for (int j = 0; j < elned; j++)
nnums[elnv+j] = nv+ednums[j]+1; nnums[elnv+j] = nv+ednums[j]+1;
nnums[elnv+elned] = fanum; nnums[elnv+elned] = fanum;
@ -219,7 +220,6 @@ namespace netgen
{ 2, 3, 1, 1, 4, 5, 1, 6, 4, 5, 5, 4, 7, 7, 7 }; { 2, 3, 1, 1, 4, 5, 1, 6, 4, 5, 5, 4, 7, 7, 7 };
// int cnt = 0; // int cnt = 0;
do do
{ {
static Timer t("update cluster, identify"); static Timer t("update cluster, identify");
@ -247,23 +247,23 @@ namespace netgen
break; break;
case TET: case TET:
case TET10: case TET10:
if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) == if (cluster_reps.Get(el.PNum(1)+1-PI0) ==
cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE)) cluster_reps.Get(el.PNum(2)+1-PI0))
clustertab = tet_cluster12; clustertab = tet_cluster12;
else if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) == else if (cluster_reps.Get(el.PNum(1)+1-PI0) ==
cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE)) cluster_reps.Get(el.PNum(3)+1-PI0))
clustertab = tet_cluster13; clustertab = tet_cluster13;
else if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) == else if (cluster_reps.Get(el.PNum(1)+1-PI0) ==
cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE)) cluster_reps.Get(el.PNum(4)+1-PI0))
clustertab = tet_cluster14; clustertab = tet_cluster14;
else if (cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE) == else if (cluster_reps.Get(el.PNum(2)+1-PI0) ==
cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE)) cluster_reps.Get(el.PNum(3)+1-PI0))
clustertab = tet_cluster23; clustertab = tet_cluster23;
else if (cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE) == else if (cluster_reps.Get(el.PNum(2)+1-PI0) ==
cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE)) cluster_reps.Get(el.PNum(4)+1-PI0))
clustertab = tet_cluster24; clustertab = tet_cluster24;
else if (cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE) == else if (cluster_reps.Get(el.PNum(3)+1-PI0) ==
cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE)) cluster_reps.Get(el.PNum(4)+1-PI0))
clustertab = tet_cluster34; clustertab = tet_cluster34;
else else
@ -286,7 +286,7 @@ namespace netgen
nnums.SetSize(elnv+elned+elnfa+1); nnums.SetSize(elnv+elned+elnfa+1);
for (int j = 0; j < elnv; j++) for (int j = 0; j < elnv; j++)
nnums[j] = el[j]+1-PointIndex::BASE; nnums[j] = el[j]+1-IndexBASE<PointIndex>();
for (int j = 0; j < elned; j++) for (int j = 0; j < elned; j++)
nnums[elnv+j] = nv+ednums[j]+1; nnums[elnv+j] = nv+ednums[j]+1;
for (int j = 0; j < elnfa; j++) for (int j = 0; j < elnfa; j++)

View File

@ -1194,8 +1194,9 @@ namespace netgen
// if (el.GetType() == TRIG && order >= 3) // if (el.GetType() == TRIG && order >= 3)
if (top.GetFaceType(facenr+1) == TRIG && order >= 3) if (top.GetFaceType(facenr+1) == TRIG && order >= 3)
{ {
NgArrayMem<int, 3> verts(3); // NgArrayMem<PointIndex, 3> verts(3);
top.GetFaceVertices (facenr+1, verts); // top.GetFaceVertices (facenr+1, verts);
auto verts = top.GetFaceVertices(facenr);
int fnums[] = { 0, 1, 2 }; int fnums[] = { 0, 1, 2 };
/* /*

View File

@ -34,12 +34,10 @@ namespace netgen
int NB(int i) const { return nb[i]; } int NB(int i) const { return nb[i]; }
int FaceNr (INDEX_3 & face) const // which face nr is it ? int FaceNr (const PointIndices<3> & face) const // which face nr is it ?
{ {
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
if (pnums[i] != face.I1() && if (pnums[i] != face[0] && pnums[i] != face[1] && pnums[i] != face[2])
pnums[i] != face.I2() &&
pnums[i] != face.I3())
return i; return i;
return 3; return 3;
} }
@ -673,7 +671,7 @@ namespace netgen
IndexSet closesphere(mesh.GetNP()); IndexSet closesphere(mesh.GetNP());
// "random" reordering of points (speeds a factor 3 - 5 !!!) // "random" reordering of points (speeds a factor 3 - 5 !!!)
NgArray<PointIndex, PointIndex::BASE, PointIndex> mixed(np); Array<PointIndex, PointIndex> mixed(np);
// int prims[] = { 11, 13, 17, 19, 23, 29, 31, 37 }; // int prims[] = { 11, 13, 17, 19, 23, 29, 31, 37 };
// int prims[] = { 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 }; // int prims[] = { 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 };
int prims[] = { 211, 223, 227, 229, 233, 239, 241, 251, 257, 263 }; int prims[] = { 211, 223, 227, 229, 233, 239, 241, 251, 257, 263 };
@ -919,7 +917,7 @@ namespace netgen
{ {
auto & tri_other = mesh.OpenElement(i_other); auto & tri_other = mesh.OpenElement(i_other);
PointIndex pi2 = tri[(edge+2)%3]; PointIndex pi2 = tri[(edge+2)%3];
PointIndex pi3 = tri_other[0]+tri_other[1]+tri_other[2] - pi0 - pi1; PointIndex pi3 = tri_other[0]-pi0+tri_other[1]-pi1+tri_other[2];
if(pi2>pi3) if(pi2>pi3)
Swap(pi2, pi3); Swap(pi2, pi3);

View File

@ -819,7 +819,7 @@ namespace netgen
for(auto sei : els) for(auto sei : els)
{ {
auto & el = tempmesh[sei]; auto & el = tempmesh[sei];
PointIndex pi2 = el[0]+el[1]+el[2] - seg[0] - seg[1]; PointIndex pi2 = el[0]-seg[0]+el[1]-seg[1]+el[2];
bool is_left = ::netgen::Area(P2(tempmesh[seg[0]]), P2(tempmesh[seg[1]]), P2(tempmesh[pi2]))>0.0; bool is_left = ::netgen::Area(P2(tempmesh[seg[0]]), P2(tempmesh[seg[1]]), P2(tempmesh[pi2]))>0.0;
POSITION pos; POSITION pos;

View File

@ -190,7 +190,7 @@ namespace netgen
MinCoords(maxextreal,maxp); MinCoords(maxextreal,maxp);
int cluster = faces->Get(fstind).Cluster(); PointIndex cluster = faces->Get(fstind).Cluster();
int sx = int((minp.X()-minext.X())/elemsize.X()+1.); int sx = int((minp.X()-minext.X())/elemsize.X()+1.);
int ex = int((maxp.X()-minext.X())/elemsize.X()+1.); int ex = int((maxp.X()-minext.X())/elemsize.X()+1.);

View File

@ -226,7 +226,7 @@ namespace netgen
*/ */
for (const Element2d & el : rule.oldels) for (const Element2d & el : rule.oldels)
for (PointIndex pi : el.PNums()) for (PointIndex pi : el.PNums())
rule.incelsonnode[pi-PointIndex::BASE]--; rule.incelsonnode[pi-IndexBASE<PointIndex>()]--;
for (int j = 1; j <= rule.newels.Size(); j++) for (int j = 1; j <= rule.newels.Size(); j++)
{ {

View File

@ -1585,11 +1585,11 @@ void MeshOptimize3d :: SwapImproveSurface (
bool found = false; bool found = false;
for(int k=0; !found && k<used_idmaps->Size(); k++) for(int k=0; !found && k<used_idmaps->Size(); k++)
{ {
if(pi2 < (*used_idmaps)[k]->Size() + PointIndex::BASE) if(pi2 < (*used_idmaps)[k]->Size() + IndexBASE<PointIndex>())
{ {
pi1other = (*(*used_idmaps)[k])[pi1]; pi1other = (*(*used_idmaps)[k])[pi1];
pi2other = (*(*used_idmaps)[k])[pi2]; pi2other = (*(*used_idmaps)[k])[pi2];
found = (pi1other != 0 && pi2other != 0 && pi1other != pi1 && pi2other != pi2); found = (pi1other.IsValid() && pi2other.IsValid() && pi1other != pi1 && pi2other != pi2);
if(found) if(found)
idnum = k; idnum = k;
} }

View File

@ -182,10 +182,10 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
Array<Point3d, PointIndex> locpoints; // local points Array<Point3d, PointIndex> locpoints; // local points
NgArray<MiniElement2d> locfaces; // local faces Array<MiniElement2d> locfaces; // local faces
Array<PointIndex, PointIndex> pindex; // mapping from local to front point numbering Array<PointIndex, PointIndex> pindex; // mapping from local to front point numbering
Array<int, PointIndex> allowpoint; // point is allowed ? Array<int, PointIndex> allowpoint; // point is allowed ?
NgArray<INDEX> findex; // mapping from local to front face numbering Array<INDEX> findex; // mapping from local to front face numbering
//INDEX_2_HASHTABLE<int> connectedpairs(100); // connecgted pairs for prism meshing //INDEX_2_HASHTABLE<int> connectedpairs(100); // connecgted pairs for prism meshing
Array<Point3d, PointIndex> plainpoints; // points in reference coordinates Array<Point3d, PointIndex> plainpoints; // points in reference coordinates
@ -195,7 +195,6 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
int j, oldnp, oldnf; int j, oldnp, oldnf;
int found; int found;
referencetransform trans; referencetransform trans;
int rotind;
Point3d inp; Point3d inp;
float err; float err;
@ -211,10 +210,10 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
// for star-shaped domain meshing // for star-shaped domain meshing
NgArray<MeshPoint, PointIndex::BASE> grouppoints; Array<MeshPoint, PointIndex> grouppoints;
NgArray<MiniElement2d> groupfaces; Array<MiniElement2d> groupfaces;
Array<PointIndex, PointIndex> grouppindex; Array<PointIndex, PointIndex> grouppindex;
NgArray<INDEX> groupfindex; Array<INDEX> groupfindex;
float minerr; float minerr;
@ -313,7 +312,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
if (loktestmode) if (loktestmode)
{ {
(*testout) << "baseel = " << baseelem << ", ind = " << findex.Get(1) << endl; (*testout) << "baseel = " << baseelem << ", ind = " << findex[0] << endl;
int pi1 = pindex[locfaces[0].PNum(1)]; int pi1 = pindex[locfaces[0].PNum(1)];
int pi2 = pindex[locfaces[0].PNum(2)]; int pi2 = pindex[locfaces[0].PNum(2)];
int pi3 = pindex[locfaces[0].PNum(3)]; int pi3 = pindex[locfaces[0].PNum(3)];
@ -368,9 +367,10 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
grouppindex, groupfindex); grouppindex, groupfindex);
bool onlytri = 1; bool onlytri = 1;
for (auto i : groupfaces.Range()) for (auto & f : groupfaces)
if (groupfaces[i].GetNP() != 3) if (f.GetNP() != 3)
onlytri = 0; onlytri = 0;
if (onlytri && groupfaces.Size() <= 20 + 2*stat.qualclass && if (onlytri && groupfaces.Size() <= 20 + 2*stat.qualclass &&
FindInnerPoint (grouppoints, groupfaces, inp) && FindInnerPoint (grouppoints, groupfaces, inp) &&
@ -379,11 +379,11 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
(*testout) << "inner point found" << endl; (*testout) << "inner point found" << endl;
for(int i = 1; i <= groupfaces.Size(); i++) for(int i = 1; i <= groupfaces.Size(); i++)
adfront -> DeleteFace (groupfindex.Get(i)); adfront -> DeleteFace (groupfindex[i-1]);
for(int i = 1; i <= groupfaces.Size(); i++) for(int i = 1; i <= groupfaces.Size(); i++)
for (j = 1; j <= locfaces.Size(); j++) for (j = 1; j <= locfaces.Size(); j++)
if (findex.Get(j) == groupfindex.Get(i)) if (findex[j-1] == groupfindex[i-1])
delfaces.Append (j); delfaces.Append (j);
@ -396,13 +396,13 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
newel.SetNP(4); newel.SetNP(4);
newel.PNum(4) = npi; newel.PNum(4) = npi;
for(int i = 1; i <= groupfaces.Size(); i++) for(int i = 0; i < groupfaces.Size(); i++)
{ {
for (j = 1; j <= 3; j++) for (j = 1; j <= 3; j++)
{ {
newel.PNum(j) = newel.PNum(j) =
adfront->GetGlobalIndex adfront->GetGlobalIndex
(grouppindex[groupfaces.Get(i).PNum(j)]); (grouppindex[groupfaces[i].PNum(j)]);
} }
mesh.AddVolumeElement (newel); mesh.AddVolumeElement (newel);
} }
@ -437,7 +437,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
bool impossible = 1; bool impossible = 1;
for (rotind = 1; rotind <= locfaces[0].GetNP(); rotind++) for (int rotind = 1; rotind <= locfaces[0].GetNP(); rotind++)
{ {
// set transformatino to reference coordinates // set transformatino to reference coordinates
@ -687,26 +687,26 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
stat.cntelem++; stat.cntelem++;
} }
for(int i = oldnf+1; i <= locfaces.Size(); i++) for(int i = oldnf; i < locfaces.Size(); i++)
{ {
for (j = 1; j <= locfaces.Get(i).GetNP(); j++) for (j = 1; j <= locfaces[i].GetNP(); j++)
locfaces.Elem(i).PNum(j) = locfaces[i].PNum(j) =
pindex[locfaces.Get(i).PNum(j)]; pindex[locfaces[i].PNum(j)];
// (*testout) << "add face " << locfaces.Get(i) << endl; // (*testout) << "add face " << locfaces.Get(i) << endl;
adfront->AddFace (locfaces.Get(i)); adfront->AddFace (locfaces[i]);
} }
for(int i = 1; i <= delfaces.Size(); i++) for(int i = 1; i <= delfaces.Size(); i++)
adfront->DeleteFace (findex.Get(delfaces.Get(i))); adfront->DeleteFace (findex[delfaces.Get(i)-1]);
} }
else else
{ {
adfront->IncrementClass (findex.Get(1)); adfront->IncrementClass (findex[0]);
if (impossible && mp.check_impossible) if (impossible && mp.check_impossible)
{ {
(*testout) << "skip face since it is impossible" << endl; (*testout) << "skip face since it is impossible" << endl;
for (j = 0; j < 100; j++) for (j = 0; j < 100; j++)
adfront->IncrementClass (findex.Get(1)); adfront->IncrementClass (findex[0]);
} }
} }
@ -938,7 +938,7 @@ void Meshing3 :: BlockFill (Mesh & mesh, double gh)
for(int i = 1; i <= n; i++) for(int i = 1; i <= n; i++)
{ {
pointnr.Elem(i) = 0; pointnr.Elem(i) = PointIndex::INVALID;
frontpointnr.Elem(i) = 0; frontpointnr.Elem(i) = 0;
} }
@ -954,7 +954,7 @@ void Meshing3 :: BlockFill (Mesh & mesh, double gh)
for (j3 = i3; j3 <= i3+1; j3++) for (j3 = i3; j3 <= i3+1; j3++)
{ {
j = j3 + (j2-1) * n3 + (j1-1) * n2 * n3; j = j3 + (j2-1) * n3 + (j1-1) * n2 * n3;
if (pointnr.Get(j) == 0) if (!pointnr.Get(j).IsValid())
{ {
Point3d hp(xmin + (j1-1) * gh, Point3d hp(xmin + (j1-1) * gh,
ymin + (j2-1) * gh, ymin + (j2-1) * gh,

View File

@ -47,7 +47,7 @@ public:
/// ///
int ApplyRules (Array<Point3d, PointIndex> & lpoints, int ApplyRules (Array<Point3d, PointIndex> & lpoints,
Array<int, PointIndex> & allowpoint, Array<int, PointIndex> & allowpoint,
NgArray<MiniElement2d> & lfaces, INDEX lfacesplit, Array<MiniElement2d> & lfaces, INDEX lfacesplit,
INDEX_2_HASHTABLE<int> & connectedpairs, INDEX_2_HASHTABLE<int> & connectedpairs,
NgArray<Element> & elements, NgArray<Element> & elements,
NgArray<INDEX> & delfaces, int tolerance, NgArray<INDEX> & delfaces, int tolerance,

View File

@ -221,7 +221,6 @@ namespace netgen
friend istream & operator>> (istream &, PointIndex &); friend istream & operator>> (istream &, PointIndex &);
friend ostream & operator<< (ostream &, const 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);
@ -235,7 +234,6 @@ 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(PointIndex::BASE-1) { ; } constexpr PointIndex (t_invalid inv) : i(PointIndex::BASE-1) { ; }
// PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; } // PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; }
@ -250,6 +248,7 @@ namespace netgen
PointIndex operator+= (int add) { i += add; return *this; } PointIndex operator+= (int add) { i += add; return *this; }
void Invalidate() { i = PointIndex::BASE-1; } void Invalidate() { i = PointIndex::BASE-1; }
bool IsValid() const { return i != PointIndex::BASE-1; } bool IsValid() const { return i != PointIndex::BASE-1; }
// operator bool() const { return IsValid(); }
#ifdef BASE0 #ifdef BASE0
static constexpr size_t BASE = 0; static constexpr size_t BASE = 0;
#else #else
@ -308,6 +307,19 @@ namespace netgen
template <size_t J> template <size_t J>
PointIndex get() const { return PointIndex(INDEX_2::operator[](J)); } PointIndex get() const { return PointIndex(INDEX_2::operator[](J)); }
}; };
template <> class PointIndices<3> : public INDEX_3
{
public:
PointIndices () = default;
PointIndices (INDEX_3 i3) : INDEX_3(i3) { ; }
PointIndices (PointIndex i1, PointIndex i2, PointIndex i3) : INDEX_3(i1,i2,i3) { ; }
PointIndex operator[] (int i) const { return PointIndex(INDEX_3::operator[](i)); }
PointIndex & operator[] (int i) { return reinterpret_cast<PointIndex&>(INDEX_3::operator[](i)); }
using INDEX_3::Sort;
static PointIndices Sort(PointIndex i1, PointIndex i2, PointIndex i3) { return INDEX_3::Sort(i1, i2, i3); }
template <size_t J>
PointIndex get() const { return PointIndex(INDEX_3::operator[](J)); }
};
} }
namespace std namespace std

View File

@ -51,7 +51,7 @@ int Meshing3 :: ApplyRules
( (
Array<Point3d, PointIndex> & lpoints, // in: local points, out: old+new local points Array<Point3d, PointIndex> & lpoints, // in: local points, out: old+new local points
Array<int, PointIndex> & allowpoint, // in: 2 .. it is allowed to use pointi, 1..will be allowed later, 0..no means Array<int, PointIndex> & allowpoint, // in: 2 .. it is allowed to use pointi, 1..will be allowed later, 0..no means
NgArray<MiniElement2d> & lfaces, // in: local faces, out: old+new local faces Array<MiniElement2d> & lfaces, // in: local faces, out: old+new local faces
INDEX lfacesplit, // for local faces in outer radius INDEX lfacesplit, // for local faces in outer radius
INDEX_2_HASHTABLE<int> & connectedpairs, // connected pairs for prism-meshing INDEX_2_HASHTABLE<int> & connectedpairs, // connected pairs for prism-meshing
NgArray<Element> & elements, // out: new elements NgArray<Element> & elements, // out: new elements
@ -314,7 +314,7 @@ int Meshing3 :: ApplyRules
if (fnearness.Get(locfi) > rule->GetFNearness (nfok) || if (fnearness.Get(locfi) > rule->GetFNearness (nfok) ||
fused.Get(locfi) || fused.Get(locfi) ||
actfnp != lfaces.Get(locfi).GetNP() ) actfnp != lfaces[locfi-1].GetNP() )
{ {
// face not feasible in any rotation // face not feasible in any rotation
@ -325,7 +325,7 @@ int Meshing3 :: ApplyRules
ok = 1; ok = 1;
locface = &lfaces.Get(locfi); locface = &lfaces[locfi-1];
// reference point already mapped differently ? // reference point already mapped differently ?
@ -528,7 +528,7 @@ int Meshing3 :: ApplyRules
(*testout) << pi << " "; (*testout) << pi << " ";
(*testout) << endl; (*testout) << endl;
snprintf (problems.Elem(ri), 255, "mapping found"); snprintf (problems.Elem(ri), 255, "mapping found");
(*testout) << rule->GetNP(1) << " = " << lfaces.Get(1).GetNP() << endl; (*testout) << rule->GetNP(1) << " = " << lfaces[0].GetNP() << endl;
} }
ok = 1; ok = 1;
@ -691,7 +691,7 @@ int Meshing3 :: ApplyRules
if (!fused.Get(i)) if (!fused.Get(i))
{ {
int triin; int triin;
const MiniElement2d & lfacei = lfaces.Get(i); const MiniElement2d & lfacei = lfaces[i-1];
if (!triboxes.Elem(i).Intersect (rule->fzbox)) if (!triboxes.Elem(i).Intersect (rule->fzbox))
triin = 0; triin = 0;
@ -744,19 +744,19 @@ int Meshing3 :: ApplyRules
if (loktestmode) if (loktestmode)
{ {
(*testout) << "El with " << lfaces.Get(i).GetNP() << " points in freezone: " (*testout) << "El with " << lfaces[i-1].GetNP() << " points in freezone: "
<< lfaces.Get(i).PNum(1) << " - " << lfaces[i-1].PNum(1) << " - "
<< lfaces.Get(i).PNum(2) << " - " << lfaces[i-1].PNum(2) << " - "
<< lfaces.Get(i).PNum(3) << " - " << lfaces[i-1].PNum(3) << " - "
<< lfaces.Get(i).PNum(4) << endl; << lfaces[i-1].PNum(4) << endl;
for (int lj = 1; lj <= lfaces.Get(i).GetNP(); lj++) for (int lj = 1; lj <= lfaces[i-1].GetNP(); lj++)
(*testout) << lpoints[lfaces.Get(i).PNum(lj)] << " "; (*testout) << lpoints[lfaces[i-1].PNum(lj)] << " ";
(*testout) << endl; (*testout) << endl;
sprintf (problems.Elem(ri), "triangle (%d, %d, %d) in Freezone", sprintf (problems.Elem(ri), "triangle (%d, %d, %d) in Freezone",
lfaces.Get(i).PNum(1), lfaces.Get(i).PNum(2), lfaces[i-1].PNum(1), lfaces[i-1].PNum(2),
lfaces.Get(i).PNum(3)); lfaces[i-1].PNum(3));
} }
#else #else
if (loktestmode) if (loktestmode)
@ -789,9 +789,9 @@ int Meshing3 :: ApplyRules
<< endl; << endl;
snprintf (problems.Elem(ri), 255, "triangle (%d, %d, %d) in Freezone", snprintf (problems.Elem(ri), 255, "triangle (%d, %d, %d) in Freezone",
int(lfaces.Get(i).PNum(1)), int(lfaces[i-1].PNum(1)),
int(lfaces.Get(i).PNum(2)), int(lfaces[i-1].PNum(2)),
int(lfaces.Get(i).PNum(3))); int(lfaces[i-1].PNum(3)));
} }
hc = 0; hc = 0;
@ -802,9 +802,9 @@ int Meshing3 :: ApplyRules
rule->GetPointNr(k, 3) <= rule->GetNOldP()) rule->GetPointNr(k, 3) <= rule->GetNOldP())
{ {
for (int j = 1; j <= 3; j++) for (int j = 1; j <= 3; j++)
if (lfaces.Get(i).PNumMod(j ) == pmap.Get(rule->GetPointNr(k, 1)) && if (lfaces[i-1].PNumMod(j ) == pmap.Get(rule->GetPointNr(k, 1)) &&
lfaces.Get(i).PNumMod(j+1) == pmap.Get(rule->GetPointNr(k, 3)) && lfaces[i-1].PNumMod(j+1) == pmap.Get(rule->GetPointNr(k, 3)) &&
lfaces.Get(i).PNumMod(j+2) == pmap.Get(rule->GetPointNr(k, 2))) lfaces[i-1].PNumMod(j+2) == pmap.Get(rule->GetPointNr(k, 2)))
{ {
fmapi.Elem(k) = i; fmapi.Elem(k) = i;
hc = 1; hc = 1;
@ -827,14 +827,14 @@ int Meshing3 :: ApplyRules
if (loktestmode) if (loktestmode)
{ {
(*testout) << "Triangle in freezone: " (*testout) << "Triangle in freezone: "
<< lfaces.Get(i).PNum(1) << " - " << lfaces[i-1].PNum(1) << " - "
<< lfaces.Get(i).PNum(2) << " - " << lfaces[i-1].PNum(2) << " - "
<< lfaces.Get(i).PNum(3) << endl; << lfaces[i-1].PNum(3) << endl;
snprintf (problems.Elem(ri), 255, "triangle (%d, %d, %d) in Freezone", snprintf (problems.Elem(ri), 255, "triangle (%d, %d, %d) in Freezone",
int (lfaces.Get(i).PNum(1)), int (lfaces[i-1].PNum(1)),
int (lfaces.Get(i).PNum(2)), int (lfaces[i-1].PNum(2)),
int (lfaces.Get(i).PNum(3))); int (lfaces[i-1].PNum(3)));
} }
ok = 0; ok = 0;
} }

View File

@ -137,6 +137,8 @@ public:
DLL_HEADER void GetFaceVertices (int fnr, NgArray<int> & vertices) const; DLL_HEADER void GetFaceVertices (int fnr, NgArray<int> & vertices) const;
DLL_HEADER void GetFaceVertices (int fnr, int * vertices) const; DLL_HEADER void GetFaceVertices (int fnr, int * vertices) const;
auto GetFaceVertices (int fnr) const
{ return FlatArray (face2vert[fnr][3].IsValid() ? 4 : 3, &face2vert[fnr][0]); }
[[deprecated("use GetEdgeVertices -> tupe(v0,v1) instead")]] [[deprecated("use GetEdgeVertices -> tupe(v0,v1) instead")]]
DLL_HEADER void GetEdgeVertices (int enr, int & v1, int & v2) const; DLL_HEADER void GetEdgeVertices (int enr, int & v1, int & v2) const;
[[deprecated("use GetEdgeVertices -> tupe(v0,v1) instead")]] [[deprecated("use GetEdgeVertices -> tupe(v0,v1) instead")]]