diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index afdf8fa1..3e42ca32 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -42,11 +42,11 @@ namespace netgen return 3; } - INDEX_3 GetFace (int i) const + PointIndices<3> GetFace (int i) const { - return INDEX_3 (pnums[deltetfaces[i][0]], - pnums[deltetfaces[i][1]], - pnums[deltetfaces[i][2]]); + return { pnums[deltetfaces[i][0]], + pnums[deltetfaces[i][1]], + pnums[deltetfaces[i][2]] }; } void GetFace (int i, Element2d & face) const @@ -372,10 +372,10 @@ namespace netgen } else { - INDEX_3 i3 = tempels.Get(helind).GetFace (k); - const Point<3> & p1 = mesh[PointIndex (i3.I1())]; - const Point<3> & p2 = mesh[PointIndex (i3.I2())]; - const Point<3> & p3 = mesh[PointIndex (i3.I3())]; + PointIndices<3> i3 = tempels.Get(helind).GetFace (k); + const Point<3> & p1 = mesh[i3[0]]; + const Point<3> & p2 = mesh[i3[1]]; + const Point<3> & p3 = mesh[i3[2]]; Vec<3> n = Cross (p2-p1, p3-p1); n /= n.Length(); @@ -613,9 +613,12 @@ namespace netgen for (auto & face : adfront->Faces()) for (PointIndex pi : face.Face().PNums()) usep[pi] = true; - + + /* for (size_t i = oldnp + PointIndex::BASE; i < np + PointIndex::BASE; i++) + */ + for (auto i : mesh.Points().Range().Modify(oldnp, -4)) usep[i] = true; for (PointIndex pi : mesh.LockedPoints()) @@ -682,24 +685,25 @@ namespace netgen while (np % prims[i] == 0) i++; prim = prims[i]; } - + // for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End()-4; pi++) for (PointIndex pi : mesh.Points().Range().Modify(0, -4)) - mixed[pi] = PointIndex ( (prim * pi) % np + PointIndex::BASE ); + // mixed[pi] = PointIndex ( (prim * pi) % np + PointIndex::BASE ); + mixed[pi] = (prim * (pi-IndexBASE()+1)) % np + IndexBASE() ; Array newels; // for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End()-4; pi++) for (PointIndex pi : mesh.Points().Range().Modify(0, -4)) { - if (pi % 1000 == 0) + if ((pi-IndexBASE()) % 1000 == 0) { - if (pi % 10000 == 0) + if ((pi-IndexBASE()) % 10000 == 0) PrintDot ('+'); else PrintDot ('.'); } - multithread.percent = 100.0 * pi / np; + multithread.percent = 100.0 * (pi-IndexBASE()) / np; if (multithread.terminate) break; @@ -719,7 +723,7 @@ namespace netgen } for (int i = tempels.Size(); i >= 1; i--) - if (tempels.Get(i)[0] <= 0) + if (!tempels.Get(i)[0].IsValid()) tempels.DeleteElement (i); PrintDot ('\n'); @@ -745,7 +749,7 @@ namespace netgen { static Timer tdegenerated("Delaunay - remove degenerated"); RegionTimer rt(tdegenerated); - NgBitArray badnode(points.Size()); + TBitArray badnode(points.Size()); badnode.Clear(); int ndeg = 0; @@ -767,13 +771,15 @@ namespace netgen double h = v1.Length() + v2.Length() + v3.Length(); if (fabs (vol) < 1e-8 * (h * h * h) && - (el[0] <= np && el[1] <= np && - el[2] <= np && el[3] <= np) ) // old: 1e-12 + (el[0] < IndexBASE()+np && + el[1] < IndexBASE()+np && + el[2] < IndexBASE()+np && + el[3] < IndexBASE()+np) ) // old: 1e-12 { - badnode.Set(el[0]); - badnode.Set(el[1]); - badnode.Set(el[2]); - badnode.Set(el[3]); + badnode.SetBitAtomic(el[0]); + badnode.SetBitAtomic(el[1]); + badnode.SetBitAtomic(el[2]); + badnode.SetBitAtomic(el[3]); ndeg++; (*testout) << "vol = " << vol << " h = " << h << endl; } @@ -803,7 +809,7 @@ namespace netgen static Timer topenel("Delaunay - find openel"); RegionTimer rt(topenel); // find surface triangles which are no face of any tet - BitArray bnd_points( mesh.GetNP() + PointIndex::BASE ); + TBitArray bnd_points( mesh.GetNP() + PointIndex::BASE ); bnd_points.Clear(); for (int i = 1; i <= mesh.GetNOpenElements(); i++) @@ -849,7 +855,7 @@ namespace netgen for(auto ei : tets_with_3_bnd_points) for(auto j : Range(4)) { - auto i3_ = tempels[ei].GetFace (j); + PointIndices<3> i3_ = tempels[ei].GetFace (j); ngcore::IVec<3> i3 = {i3_[0], i3_[1], i3_[2]}; if(bnd_points[i3[0]] && bnd_points[i3[1]] && bnd_points[i3[2]]) { @@ -865,7 +871,7 @@ namespace netgen for (int i = 1; i <= mesh.GetNOpenElements(); i++) { const Element2d & tri = mesh.OpenElement(i); - ngcore::IVec<3> i3(tri[0], tri[1], tri[2]); + ngcore::IVec<3,PointIndex> i3(tri[0], tri[1], tri[2]); i3.Sort(); if(!face_table.Used(i3)) openels.Append(i); @@ -883,7 +889,7 @@ namespace netgen table.Add(tri[2], openel_i); }, mesh.GetNP()); - ngcore::BitArray badnode(mesh.GetNP()+PointIndex::BASE); + TBitArray badnode(mesh.GetNP()+PointIndex::BASE); badnode.Clear(); ngcore::ParallelForRange(openels.Size(), [&] (auto myrange) { @@ -941,8 +947,8 @@ namespace netgen double h = v1.Length() + v2.Length() + v3.Length(); if (fabs (vol) < 1e-4 * (h * h * h)) // old: 1e-12 { - badnode.SetBitAtomic(pi2); - badnode.SetBitAtomic(pi3); + badnode.SetBitAtomic(pi2); + badnode.SetBitAtomic(pi3); } break; } @@ -1106,7 +1112,7 @@ namespace netgen */ for (const Element2d & tri : mesh.OpenElements()) { - INDEX_3 i3 (tri[0], tri[1], tri[2]); + PointIndices<3> i3 (tri[0], tri[1], tri[2]); i3.Sort(); boundaryfaces.PrepareSet (i3); } @@ -1114,7 +1120,7 @@ namespace netgen for (int i = 1; i <= mesh.GetNOpenElements(); i++) { const Element2d & tri = mesh.OpenElement(i); - INDEX_3 i3 (tri[0], tri[1], tri[2]); + PointIndices<3> i3 (tri[0], tri[1], tri[2]); i3.Sort(); boundaryfaces.Set (i3, 1); } @@ -1136,7 +1142,7 @@ namespace netgen */ for (const DelaunayTet & el : tempels) { - INDEX_4 i4(el[0], el[1], el[2], el[3]); + PointIndices<4> i4(el[0], el[1], el[2], el[3]); i4.Sort(); elsonpoint.IncSizePrepare (i4.I1()); elsonpoint.IncSizePrepare (i4.I2()); @@ -1147,7 +1153,7 @@ namespace netgen for (int i = 0; i < tempels.Size(); i++) { const DelaunayTet & el = tempels[i]; - INDEX_4 i4(el[0], el[1], el[2], el[3]); + PointIndices<4> i4(el[0], el[1], el[2], el[3]); i4.Sort(); elsonpoint.Add (i4.I1(), i+1); elsonpoint.Add (i4.I2(), i+1); @@ -1176,7 +1182,7 @@ namespace netgen if (hel[0] == pi) { - INDEX_3 i3(hel[0], hel[1], hel[2]); + PointIndices<3> i3(hel[0], hel[1], hel[2]); if (!boundaryfaces.Used (i3)) { @@ -1191,7 +1197,7 @@ namespace netgen { hel.Invert(); hel.NormalizeNumbering(); - INDEX_3 i3i(hel[0], hel[1], hel[2]); + PointIndices<3> i3i(hel[0], hel[1], hel[2]); INDEX_2 i2(i, j); faceht.Set (i3i, i2); } diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 939b1e7b..017432a5 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -234,7 +234,6 @@ namespace netgen friend bool operator== (PointIndex a, PointIndex b); friend bool operator!= (PointIndex a, PointIndex b); */ - public: constexpr PointIndex (t_invalid inv) : i(PointIndex::BASE-1) { ; } // PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; } @@ -258,7 +257,6 @@ namespace netgen void DoArchive (Archive & ar) { ar & 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); } @@ -321,6 +319,21 @@ namespace netgen template PointIndex get() const { return PointIndex(INDEX_3::operator[](J)); } }; + + template <> class PointIndices<4> : public INDEX_4 + { + public: + PointIndices () = default; + PointIndices (INDEX_4 i4) : INDEX_4(i4) { ; } + PointIndices (PointIndex i1, PointIndex i2, PointIndex i3, PointIndex i4) : INDEX_4(i1,i2,i3,i4) { ; } + PointIndex operator[] (int i) const { return PointIndex(INDEX_4::operator[](i)); } + PointIndex & operator[] (int i) { return reinterpret_cast(INDEX_4::operator[](i)); } + using INDEX_4::Sort; + // static PointIndices Sort(PointIndex i1, PointIndex i2, PointIndex i3, PointIndex i4) { return INDEX_4::Sort(i1, i2, i3, i4); } + template + PointIndex get() const { return PointIndex(INDEX_4::operator[](J)); } + }; + } namespace std diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp index ac4f7128..053b8ac9 100644 --- a/libsrc/occ/occgenmesh.cpp +++ b/libsrc/occ/occgenmesh.cpp @@ -300,7 +300,7 @@ namespace netgen } for(const auto& vert : geom.GetFaceVertices(geom.GetFace(k-1))) { - PointIndex pi = vert->nr + 1; + PointIndex pi = vert->nr + IndexBASE(); if(glob2loc[pi] == 0) { auto gi = occface.Project(mesh[pi]); @@ -398,7 +398,7 @@ namespace netgen } for(const auto& vert : geom.GetFaceVertices(geom.GetFace(k-1))) { - PointIndex pi = vert->nr + 1; + PointIndex pi = vert->nr + IndexBASE(); if(glob2loc[pi] == 0) { auto gi = occface.Project(mesh[pi]);