PointIndex in bisect

This commit is contained in:
Joachim Schoeberl 2024-12-23 19:24:48 +01:00
parent 31ed810144
commit 73bcb1bd29
10 changed files with 99 additions and 77 deletions

View File

@ -210,6 +210,24 @@ private:
NGCORE_API std::ostream & operator<<(std::ostream & s, const BitArray & ba); NGCORE_API std::ostream & operator<<(std::ostream & s, const BitArray & ba);
template <typename IndexType>
class TBitArray : public BitArray
{
public:
using BitArray::BitArray;
void SetBit (IndexType i) { BitArray::SetBit(i-IndexBASE<IndexType>()); }
void Clear () { BitArray::Clear(); }
void Clear (IndexType i) { BitArray::Clear(i-IndexBASE<IndexType>()); }
void SetBitAtomic (IndexType i) { BitArray::SetBitAtomic(i-IndexBASE<IndexType>()); }
bool Test (IndexType i) const { return BitArray::Test(i-IndexBASE<IndexType>()); }
bool operator[] (IndexType i) const { return Test(i); }
};
} // namespace ngcore } // namespace ngcore
#endif // NETGEN_CORE_BITARRAY #endif // NETGEN_CORE_BITARRAY

View File

@ -349,12 +349,9 @@ NGX_INLINE DLL_HEADER Ng_Buffer<int[2]> Ngx_Mesh :: GetPeriodicVertices(int idnr
NGX_INLINE void Ngx_Mesh :: GetParentNodes (int ni, int * parents) const NGX_INLINE void Ngx_Mesh :: GetParentNodes (int ni, int * parents) const
{ {
ni++; if (ni < mesh->mlbetweennodes.Size())
if (ni <= mesh->mlbetweennodes.Size()) for (int j = 0; j < 2; j++)
{ parents[j] = mesh->mlbetweennodes[IndexBASE<PointIndex>()+ni][j] - IndexBASE<PointIndex>();
parents[0] = mesh->mlbetweennodes.Get(ni).I1()-1;
parents[1] = mesh->mlbetweennodes.Get(ni).I2()-1;
}
else else
parents[0] = parents[1] = -1; parents[0] = parents[1] = -1;
} }

View File

@ -1736,8 +1736,8 @@ void Ng_GetParentNodes (int ni, int * parents)
{ {
if (ni <= mesh->mlbetweennodes.Size()) if (ni <= mesh->mlbetweennodes.Size())
{ {
parents[0] = mesh->mlbetweennodes.Get(ni).I1(); parents[0] = mesh->mlbetweennodes[ni].I1();
parents[1] = mesh->mlbetweennodes.Get(ni).I2(); parents[1] = mesh->mlbetweennodes[ni].I2();
} }
else else
parents[0] = parents[1] = 0; parents[0] = parents[1] = 0;

View File

@ -3230,7 +3230,7 @@ namespace netgen
if (opt.refine_hp) if (opt.refine_hp)
{ {
PrintMessage(3,"refine hp"); PrintMessage(3,"refine hp");
NgBitArray singv(np); TBitArray<PointIndex> singv(np);
singv.Clear(); singv.Clear();
if (mesh.GetDimension() == 3) if (mesh.GetDimension() == 3)
@ -3238,8 +3238,8 @@ namespace netgen
for (int i = 1; i <= mesh.GetNSeg(); i++) for (int i = 1; i <= mesh.GetNSeg(); i++)
{ {
const Segment & seg = mesh.LineSegment(i); const Segment & seg = mesh.LineSegment(i);
singv.Set (seg[0]); singv.SetBit (seg[0]);
singv.Set (seg[1]); singv.SetBit (seg[1]);
} }
/* /*
for ( i=1; i<= mesh.GetNSE(); i++) for ( i=1; i<= mesh.GetNSE(); i++)
@ -3253,18 +3253,18 @@ namespace netgen
else else
{ {
// vertices with 2 different bnds // vertices with 2 different bnds
NgArray<int> bndind(np); Array<int,PointIndex> bndind(np);
bndind = 0; bndind = 0;
for (int i = 1; i <= mesh.GetNSeg(); i++) for (int i = 1; i <= mesh.GetNSeg(); i++)
{ {
const Segment & seg = mesh.LineSegment(i); const Segment & seg = mesh.LineSegment(i);
for (int j = 0; j < 2; j++) for (int j = 0; j < 2; j++)
{ {
int pi = (j == 0) ? seg[0] : seg[1]; PointIndex pi = (j == 0) ? seg[0] : seg[1];
if (bndind.Elem(pi) == 0) if (bndind[pi] == 0)
bndind.Elem(pi) = seg.edgenr; bndind[pi] = seg.edgenr;
else if (bndind.Elem(pi) != seg.edgenr) else if (bndind[pi] != seg.edgenr)
singv.Set (pi); singv.SetBit (pi);
} }
} }
} }
@ -3339,7 +3339,7 @@ namespace netgen
{ {
MarkedTet oldtet = mtets[i]; MarkedTet oldtet = mtets[i];
INDEX_2 edge(oldtet.pnums[oldtet.tetedge1], PointIndices<2> edge(oldtet.pnums[oldtet.tetedge1],
oldtet.pnums[oldtet.tetedge2]); oldtet.pnums[oldtet.tetedge2]);
edge.Sort(); edge.Sort();
@ -3350,8 +3350,7 @@ namespace netgen
} }
else else
{ {
Point<3> npt = Center (mesh.Point (edge.I1()), Point<3> npt = Center (mesh[edge[0]], mesh[edge[1]]);
mesh.Point (edge.I2()));
newp = mesh.AddPoint (npt); newp = mesh.AddPoint (npt);
cutedges.Set (edge, newp); cutedges.Set (edge, newp);
} }
@ -3379,9 +3378,9 @@ namespace netgen
pi1++; pi1++;
int pi2 = 3-pi1-oldprism.markededge; int pi2 = 3-pi1-oldprism.markededge;
INDEX_2 edge1(oldprism.pnums[pi1], PointIndices<2> edge1(oldprism.pnums[pi1],
oldprism.pnums[pi2]); oldprism.pnums[pi2]);
INDEX_2 edge2(oldprism.pnums[pi1+3], PointIndices<2> edge2(oldprism.pnums[pi1+3],
oldprism.pnums[pi2+3]); oldprism.pnums[pi2+3]);
edge1.Sort(); edge1.Sort();
edge2.Sort(); edge2.Sort();
@ -3390,8 +3389,7 @@ namespace netgen
newp1 = cutedges.Get(edge1); newp1 = cutedges.Get(edge1);
else else
{ {
Point<3> npt = Center (mesh.Point (edge1.I1()), Point<3> npt = Center (mesh[edge1[0]], mesh[edge1[1]]);
mesh.Point (edge1.I2()));
newp1 = mesh.AddPoint (npt); newp1 = mesh.AddPoint (npt);
cutedges.Set (edge1, newp1); cutedges.Set (edge1, newp1);
} }
@ -3399,8 +3397,7 @@ namespace netgen
newp2 = cutedges.Get(edge2); newp2 = cutedges.Get(edge2);
else else
{ {
Point<3> npt = Center (mesh.Point (edge2.I1()), Point<3> npt = Center (mesh[edge2[0]], mesh[edge2[1]]);
mesh.Point (edge2.I2()));
newp2 = mesh.AddPoint (npt); newp2 = mesh.AddPoint (npt);
cutedges.Set (edge2, newp2); cutedges.Set (edge2, newp2);
} }
@ -3422,18 +3419,22 @@ namespace netgen
oldid = mids.Get(i); oldid = mids.Get(i);
NgArray<INDEX_2> edges; NgArray<PointIndices<2>> edges;
edges.Append(INDEX_2(oldid.pnums[oldid.markededge], edges.Append( {
oldid.pnums[(oldid.markededge+1)%oldid.np])); oldid.pnums[oldid.markededge],
edges.Append(INDEX_2(oldid.pnums[oldid.markededge + oldid.np], oldid.pnums[(oldid.markededge+1)%oldid.np] } );
oldid.pnums[(oldid.markededge+1)%oldid.np + oldid.np])); edges.Append( {
oldid.pnums[oldid.markededge + oldid.np],
oldid.pnums[(oldid.markededge+1)%oldid.np + oldid.np] } );
if(oldid.np == 4) if(oldid.np == 4)
{ {
edges.Append(INDEX_2(oldid.pnums[(oldid.markededge+2)%oldid.np], edges.Append( {
oldid.pnums[(oldid.markededge+3)%oldid.np])); oldid.pnums[(oldid.markededge+2)%oldid.np],
edges.Append(INDEX_2(oldid.pnums[(oldid.markededge+2)%oldid.np + oldid.np], oldid.pnums[(oldid.markededge+3)%oldid.np]} );
oldid.pnums[(oldid.markededge+3)%oldid.np + oldid.np])); edges.Append( {
oldid.pnums[(oldid.markededge+2)%oldid.np + oldid.np],
oldid.pnums[(oldid.markededge+3)%oldid.np + oldid.np] } );
} }
for (int j = 0; j < edges.Size(); j++) for (int j = 0; j < edges.Size(); j++)
{ {
@ -3479,7 +3480,7 @@ namespace netgen
MarkedTri oldtri = mtris[i]; MarkedTri oldtri = mtris[i];
PointIndex oldpi1 = oldtri.pnums[(oldtri.markededge+1)%3]; PointIndex oldpi1 = oldtri.pnums[(oldtri.markededge+1)%3];
PointIndex oldpi2 = oldtri.pnums[(oldtri.markededge+2)%3]; PointIndex oldpi2 = oldtri.pnums[(oldtri.markededge+2)%3];
INDEX_2 edge(oldpi1, oldpi2); PointIndices<2> edge(oldpi1, oldpi2);
edge.Sort(); edge.Sort();
int si = mesh.GetFaceDescriptor (oldtri.surfid).SurfNr(); int si = mesh.GetFaceDescriptor (oldtri.surfid).SurfNr();
@ -3529,21 +3530,21 @@ namespace netgen
INDEX_2 edge2(oldquad.pnums[2], INDEX_2 edge2(oldquad.pnums[2],
oldquad.pnums[3]); oldquad.pnums[3]);
*/ */
INDEX_2 edge1, edge2; PointIndices<2> edge1, edge2;
PointGeomInfo pgi11, pgi12, pgi21, pgi22; PointGeomInfo pgi11, pgi12, pgi21, pgi22;
if (oldquad.markededge==0 || oldquad.markededge==2) if (oldquad.markededge==0 || oldquad.markededge==2)
{ {
edge1.I1()=oldquad.pnums[0]; pgi11=oldquad.pgeominfo[0]; edge1[0] = oldquad.pnums[0]; pgi11=oldquad.pgeominfo[0];
edge1.I2()=oldquad.pnums[1]; pgi12=oldquad.pgeominfo[1]; edge1[1] = oldquad.pnums[1]; pgi12=oldquad.pgeominfo[1];
edge2.I1()=oldquad.pnums[2]; pgi21=oldquad.pgeominfo[2]; edge2[0] = oldquad.pnums[2]; pgi21=oldquad.pgeominfo[2];
edge2.I2()=oldquad.pnums[3]; pgi22=oldquad.pgeominfo[3]; edge2[1] = oldquad.pnums[3]; pgi22=oldquad.pgeominfo[3];
} }
else // 3 || 1 else // 3 || 1
{ {
edge1.I1()=oldquad.pnums[0]; pgi11=oldquad.pgeominfo[0]; edge1[0] = oldquad.pnums[0]; pgi11=oldquad.pgeominfo[0];
edge1.I2()=oldquad.pnums[2]; pgi12=oldquad.pgeominfo[2]; edge1[1] = oldquad.pnums[2]; pgi12=oldquad.pgeominfo[2];
edge2.I1()=oldquad.pnums[1]; pgi21=oldquad.pgeominfo[1]; edge2[0] = oldquad.pnums[1]; pgi21=oldquad.pgeominfo[1];
edge2.I2()=oldquad.pnums[3]; pgi22=oldquad.pgeominfo[3]; edge2[1] = oldquad.pnums[3]; pgi22=oldquad.pgeominfo[3];
} }
edge1.Sort(); edge1.Sort();
@ -3555,8 +3556,8 @@ namespace netgen
} }
else else
{ {
Point<3> np1 = Center (mesh.Point (edge1.I1()), Point<3> np1 = Center (mesh.Point (edge1[0]),
mesh.Point (edge1.I2())); mesh.Point (edge1[1]));
newp1 = mesh.AddPoint (np1); newp1 = mesh.AddPoint (np1);
cutedges.Set (edge1, newp1); cutedges.Set (edge1, newp1);
} }
@ -3567,8 +3568,8 @@ namespace netgen
} }
else else
{ {
Point<3> np2 = Center (mesh.Point (edge2.I1()), Point<3> np2 = Center (mesh.Point (edge2[0]),
mesh.Point (edge2.I2())); mesh.Point (edge2[1]));
newp2 = mesh.AddPoint (np2); newp2 = mesh.AddPoint (np2);
cutedges.Set (edge2, newp2); cutedges.Set (edge2, newp2);
} }
@ -3606,7 +3607,7 @@ namespace netgen
for (int i = 1; i <= nseg; i++) for (int i = 1; i <= nseg; i++)
{ {
Segment & seg = mesh.LineSegment (i); Segment & seg = mesh.LineSegment (i);
INDEX_2 edge(seg[0], seg[1]); PointIndices<2> edge(seg[0], seg[1]);
edge.Sort(); edge.Sort();
if (cutedges.Used (edge)) if (cutedges.Used (edge))
{ {
@ -3614,7 +3615,7 @@ namespace netgen
Segment nseg1 = seg; Segment nseg1 = seg;
Segment nseg2 = seg; Segment nseg2 = seg;
int newpi = cutedges.Get(edge); PointIndex newpi = cutedges.Get(edge);
nseg1[1] = newpi; nseg1[1] = newpi;
nseg2[0] = newpi; nseg2[0] = newpi;
@ -3670,7 +3671,7 @@ namespace netgen
if (opt.refine_hp) if (opt.refine_hp)
{ {
// //
NgArray<int> v_order (mesh.GetNP()); Array<int,PointIndex> v_order (mesh.GetNP());
v_order = 0; v_order = 0;
if (mesh.GetDimension() == 3) if (mesh.GetDimension() == 3)
{ {
@ -3680,12 +3681,12 @@ namespace netgen
for (int i = 0; i < mtets.Size(); i++) for (int i = 0; i < mtets.Size(); i++)
for (int j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
if (int(mtets[i].order) > v_order.Elem(mtets[i].pnums[j])) if (int(mtets[i].order) > v_order[mtets[i].pnums[j]])
v_order.Elem(mtets[i].pnums[j]) = mtets[i].order; v_order[mtets[i].pnums[j]] = mtets[i].order;
for (int i = 0; i < mtets.Size(); i++) for (int i = 0; i < mtets.Size(); i++)
for (int j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
if (int(mtets[i].order) < v_order.Elem(mtets[i].pnums[j])-1) if (int(mtets[i].order) < v_order[mtets[i].pnums[j]]-1)
mtets[i].order = v_order.Elem(mtets[i].pnums[j])-1; mtets[i].order = v_order[mtets[i].pnums[j]]-1;
} }
else else
{ {
@ -3697,13 +3698,13 @@ namespace netgen
for (int i = 0; i < mtris.Size(); i++) for (int i = 0; i < mtris.Size(); i++)
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
if (int(mtris[i].order) > v_order.Elem(mtris[i].pnums[j])) if (int(mtris[i].order) > v_order[mtris[i].pnums[j]])
v_order.Elem(mtris[i].pnums[j]) = mtris[i].order; v_order[mtris[i].pnums[j]] = mtris[i].order;
for (int i = 0; i < mtris.Size(); i++) for (int i = 0; i < mtris.Size(); i++)
{ {
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
if (int(mtris[i].order) < v_order.Elem(mtris[i].pnums[j])-1) if (int(mtris[i].order) < v_order[mtris[i].pnums[j]]-1)
mtris[i].order = v_order.Elem(mtris[i].pnums[j])-1; mtris[i].order = v_order[mtris[i].pnums[j]]-1;
} }
} }
} }
@ -3866,11 +3867,15 @@ namespace netgen
if (mesh.level_nv.Size() <= 1) if (mesh.level_nv.Size() <= 1)
{ {
PrintMessage(4,"RESETTING mlbetweennodes"); PrintMessage(4,"RESETTING mlbetweennodes");
/*
for (int i = 1; i <= np; i++) for (int i = 1; i <= np; i++)
{ {
mesh.mlbetweennodes.Elem(i).I1() = 0; mesh.mlbetweennodes.Elem(i).I1() = 0;
mesh.mlbetweennodes.Elem(i).I2() = 0; mesh.mlbetweennodes.Elem(i).I2() = 0;
} }
*/
for (auto i : mesh.mlbetweennodes.Range())
mesh.mlbetweennodes[i] = { PointIndex::INVALID, PointIndex::INVALID };
} }
mesh.level_nv.Append (np); mesh.level_nv.Append (np);
@ -3888,7 +3893,7 @@ namespace netgen
*/ */
NgBitArray isnewpoint(np); TBitArray<PointIndex> isnewpoint(np);
isnewpoint.Clear(); isnewpoint.Clear();
{ {
@ -3899,8 +3904,8 @@ namespace netgen
INDEX_2 edge; INDEX_2 edge;
PointIndex newpi; PointIndex newpi;
cutedges.GetData0 (i, edge, newpi); cutedges.GetData0 (i, edge, newpi);
isnewpoint.Set(newpi); isnewpoint.SetBit(newpi);
mesh.mlbetweennodes.Elem(newpi) = edge; mesh.mlbetweennodes[newpi] = edge;
} }
} }
/* /*
@ -4034,7 +4039,8 @@ namespace netgen
if(noprojection) if(noprojection)
{ {
do_repair = false; do_repair = false;
for(int ii=1; ii<=mesh.GetNP(); ii++) // for(int ii=1; ii<=mesh.GetNP(); ii++)
for (auto ii : mesh.Points().Range())
{ {
if(isnewpoint.Test(ii) && mesh.mlbetweennodes[ii][0].IsValid()) if(isnewpoint.Test(ii) && mesh.mlbetweennodes[ii][0].IsValid())
{ {

View File

@ -222,7 +222,7 @@ namespace netgen
// number of vertices on each refinement level: // number of vertices on each refinement level:
NgArray<size_t> level_nv; NgArray<size_t> level_nv;
/// refinement hierarchy /// refinement hierarchy
NgArray<PointIndices<2>,PointIndex::BASE> mlbetweennodes; Array<PointIndices<2>,PointIndex> mlbetweennodes;
/// parent element of volume element /// parent element of volume element
NgArray<int> mlparentelement; NgArray<int> mlparentelement;
/// parent element of surface element /// parent element of surface element

View File

@ -217,11 +217,11 @@ namespace netgen
// throw Exception("illegal PointIndex, use PointIndex::INVALID instead"); // throw Exception("illegal PointIndex, use PointIndex::INVALID instead");
#endif #endif
} }
/*
friend constexpr netgen::PointIndex ngcore::IndexBASE<netgen::PointIndex> (); friend constexpr netgen::PointIndex ngcore::IndexBASE<netgen::PointIndex> ();
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);
friend PointIndex operator+ (int, PointIndex); friend PointIndex operator+ (int, PointIndex);
@ -234,6 +234,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(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; }

View File

@ -955,7 +955,7 @@ namespace netgen
self.openelements = NgArray<Element2d>(0); self.openelements = NgArray<Element2d>(0);
self.opensegments = NgArray<Segment>(0); self.opensegments = NgArray<Segment>(0);
self.numvertices = 0; self.numvertices = 0;
self.mlbetweennodes = NgArray<PointIndices<2>,PointIndex::BASE> (0); self.mlbetweennodes = Array<PointIndices<2>,PointIndex> (0);
self.mlparentelement = NgArray<int>(0); self.mlparentelement = NgArray<int>(0);
self.mlparentsurfaceelement = NgArray<int>(0); self.mlparentsurfaceelement = NgArray<int>(0);
self.curvedelems = make_unique<CurvedElements> (self); self.curvedelems = make_unique<CurvedElements> (self);

View File

@ -184,7 +184,7 @@ namespace netgen
if (mesh.mlbetweennodes.Size() == mesh.Points().Size()) if (mesh.mlbetweennodes.Size() == mesh.Points().Size())
{ {
NgArray<PointIndices<2>,PointIndex::BASE> hml { mesh.mlbetweennodes }; Array<PointIndices<2>,PointIndex> hml { mesh.mlbetweennodes };
for (PointIndex pi : Range(mesh.Points())) for (PointIndex pi : Range(mesh.Points()))
mesh.mlbetweennodes[inv_index[pi]] = hml[pi]; mesh.mlbetweennodes[inv_index[pi]] = hml[pi];
} }

View File

@ -6,7 +6,7 @@
namespace netgen namespace netgen
{ {
void GetPureBadness(Mesh & mesh, NgArray<double> & pure_badness, void GetPureBadness(Mesh & mesh, NgArray<double> & pure_badness,
const NgBitArray & isnewpoint) const TBitArray<PointIndex> & isnewpoint)
{ {
//const int ne = mesh.GetNE(); //const int ne = mesh.GetNE();
const int np = mesh.GetNP(); const int np = mesh.GetNP();
@ -152,7 +152,7 @@ namespace netgen
void RepairBisection(Mesh & mesh, NgArray<ElementIndex> & bad_elements, void RepairBisection(Mesh & mesh, NgArray<ElementIndex> & bad_elements,
const NgBitArray & isnewpoint, const Refinement & refinement, const TBitArray<PointIndex> & isnewpoint, const Refinement & refinement,
const NgArray<double> & pure_badness, const NgArray<double> & pure_badness,
double max_worsening, const bool uselocalworsening, double max_worsening, const bool uselocalworsening,
const NgArray< idmap_type* > & idmaps) const NgArray< idmap_type* > & idmaps)

View File

@ -5,13 +5,13 @@ namespace netgen
{ {
void GetPureBadness(Mesh & mesh, NgArray<double> & pure_badness, void GetPureBadness(Mesh & mesh, NgArray<double> & pure_badness,
const NgBitArray & isnewpoint); const TBitArray<PointIndex> & isnewpoint);
double Validate(const Mesh & mesh, NgArray<ElementIndex> & bad_elements, double Validate(const Mesh & mesh, NgArray<ElementIndex> & bad_elements,
const NgArray<double> & pure_badness, const NgArray<double> & pure_badness,
double max_worsening, const bool uselocalworsening, double max_worsening, const bool uselocalworsening,
NgArray<double> * quality_loss = NULL); NgArray<double> * quality_loss = NULL);
void RepairBisection(Mesh & mesh, NgArray<ElementIndex> & bad_elements, void RepairBisection(Mesh & mesh, NgArray<ElementIndex> & bad_elements,
const NgBitArray & isnewpoint, const Refinement & refinement, const TBitArray<PointIndex> & isnewpoint, const Refinement & refinement,
const NgArray<double> & pure_badness, const NgArray<double> & pure_badness,
double max_worsening, const bool uselocalworsening, double max_worsening, const bool uselocalworsening,
const NgArray< idmap_type* > & idmaps); const NgArray< idmap_type* > & idmaps);