From 75032f99059c8741d748a1c7b7ac5e6dbaef7954 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Sat, 28 Dec 2024 21:26:05 +0100 Subject: [PATCH] operators +/- for PointIndex --- libsrc/geom2d/genmesh2d.cpp | 16 ++-- libsrc/meshing/bisect.cpp | 2 +- libsrc/meshing/boundarylayer.cpp | 10 +-- libsrc/meshing/boundarylayer2d.cpp | 6 +- libsrc/meshing/boundarylayer_interpolate.cpp | 2 +- libsrc/meshing/boundarylayer_limiter.hpp | 12 +-- libsrc/meshing/meshfunc.cpp | 23 ++++-- libsrc/meshing/meshtype.hpp | 60 ++++++++++---- libsrc/meshing/parallelmesh.cpp | 4 +- libsrc/meshing/parser3.cpp | 16 ++-- libsrc/meshing/smoothing3.cpp | 27 +++--- libsrc/meshing/topology.cpp | 87 ++++++++++---------- libsrc/meshing/topology.hpp | 2 +- libsrc/visualization/vssolution.cpp | 4 +- 14 files changed, 155 insertions(+), 116 deletions(-) diff --git a/libsrc/geom2d/genmesh2d.cpp b/libsrc/geom2d/genmesh2d.cpp index 54f624d0..cd3def6c 100644 --- a/libsrc/geom2d/genmesh2d.cpp +++ b/libsrc/geom2d/genmesh2d.cpp @@ -139,7 +139,9 @@ namespace netgen mark = spline.GetPoint (edgelength); { - PointIndex pi1 = -1, pi2 = -1; + PointIndex pi1{PointIndex::INVALID}; + PointIndex pi2{PointIndex::INVALID}; + Point3d mark3(mark(0), mark(1), 0); Point3d oldmark3(oldmark(0), oldmark(1), 0); @@ -157,12 +159,12 @@ namespace netgen if ( mesh[PointIndex(locsearch[k])].GetLayer() == spline.layer) pi2 = locsearch[k]; - if (pi1 == -1) + if (!pi1.IsValid()) { pi1 = mesh.AddPoint(oldmark3, spline.layer); searchtree.Insert (oldmark3, pi1); } - if (pi2 == -1) + if (!pi2.IsValid()) { pi2 = mesh.AddPoint(mark3, spline.layer); searchtree.Insert (mark3, pi2); @@ -292,13 +294,13 @@ namespace netgen Point<2> hnewp = (j == 1) ? splines[i]->StartPI() : splines[i]->EndPI(); Point<3> newp(hnewp(0), hnewp(1), 0); int layer = GetSpline(i).layer; - int npi = -1; - for (PointIndex pi = PointIndex::BASE; - pi < mesh2d.GetNP()+PointIndex::BASE; pi++) + PointIndex npi(PointIndex::INVALID); + for (PointIndex pi = IndexBASE(); + pi < mesh2d.GetNP()+IndexBASE(); pi++) if (Dist2 (mesh2d.Point(pi), newp) < 1e-12 * diam2 && mesh2d.Point(pi).GetLayer() == layer) npi = pi; - if (npi == -1) + if (!npi.IsValid()) { npi = mesh2d.AddPoint (newp, layer); searchtree.Insert (newp, npi); diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp index 038e30af..9eeb9419 100644 --- a/libsrc/meshing/bisect.cpp +++ b/libsrc/meshing/bisect.cpp @@ -547,7 +547,7 @@ namespace netgen e2[0] = (*idmaps[k])[e1[0]]; e2[1] = (*idmaps[k])[e1[1]]; - if(e2.I1() == 0 || e2.I2() == 0 || + if(!e2.I1().IsValid() || !e2.I2().IsValid() || e1.I1() == e2.I1() || e1.I2() == e2.I2()) continue; diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index 671601f3..ec92bc9b 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -87,7 +87,7 @@ Vec<3> CalcGrowthVector (FlatArray> ns) auto [maxpos1, maxpos2] = FindCloseVectors(ns); Array> new_normals; new_normals = ns; - const auto dot = ns[maxpos1] * ns[maxpos2]; + // const auto dot = ns[maxpos1] * ns[maxpos2]; auto average = 0.5 * (ns[maxpos1] + ns[maxpos2]); average.Normalize(); new_normals[maxpos1] = average; @@ -154,7 +154,7 @@ Vec<3> BoundaryLayerTool ::getEdgeTangent(PointIndex pi, int edgenr, FlatArray= np_old) || mapto[pi].Size() > 0 || special_boundary_points.count(pi); + return (pi - IndexBASE() >= np_old) || mapto[pi].Size() > 0 || special_boundary_points.count(pi); }; std::set points_set; diff --git a/libsrc/meshing/boundarylayer_limiter.hpp b/libsrc/meshing/boundarylayer_limiter.hpp index dacbdc09..d30c5445 100644 --- a/libsrc/meshing/boundarylayer_limiter.hpp +++ b/libsrc/meshing/boundarylayer_limiter.hpp @@ -288,7 +288,7 @@ struct GrowthVectorLimiter return; for (PointIndex pi : IntRange(tool.np, mesh.GetNP())) { - auto pi_from = map_from[pi]; + // auto pi_from = map_from[pi]; std::set pis; for (auto sei : p2sel[pi]) for (auto pi_ : tool.new_sels[sei].PNums()) @@ -360,7 +360,7 @@ struct GrowthVectorLimiter if (sel.GetNP() == 4) continue; - const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex()); + // const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex()); auto np = sel.GetNP(); double shift = 1.0; @@ -447,7 +447,7 @@ struct GrowthVectorLimiter for (auto sei : SurfaceElementsRange()) { const auto& sel = Get(sei); - auto sel_index = sel.GetIndex(); + // auto sel_index = sel.GetIndex(); Box<3> box(Box<3>::EMPTY_BOX); for (auto pi : sel.PNums()) @@ -466,7 +466,7 @@ struct GrowthVectorLimiter RegionTimer rt(t); BuildSearchTree(trig_shift); auto np_new = mesh.Points().Size(); - int counter = 0; + // int counter = 0; for (auto i : IntRange(tool.np, np_new)) { PointIndex pi_to = i + PointIndex::BASE; @@ -478,7 +478,7 @@ struct GrowthVectorLimiter continue; Box<3> box(Box<3>::EMPTY_BOX); - auto seg = GetSeg(pi_to, seg_shift); + // auto seg = GetSeg(pi_to, seg_shift); box.Add(GetPoint(pi_to, 0)); box.Add(GetPoint(pi_to, GetLimit(pi_from))); @@ -488,7 +488,7 @@ struct GrowthVectorLimiter return false; if (sel.PNums().Contains(pi_to)) return false; - counter++; + // counter++; f(pi_to, sei); return false; }); diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 1170a7aa..8fae3dc8 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -59,6 +59,11 @@ namespace netgen auto num_points = mesh.GetNP(); auto num_facedescriptors = mesh.GetNFD(); + + constexpr PointIndex state0 = IndexBASE()-1; + constexpr PointIndex state1 = state0+1; + constexpr PointIndex state2 = state0+2; + for(auto i : Range(ret)) { auto & md = ret[i]; @@ -73,7 +78,7 @@ namespace netgen m.SetLocalH(mesh.GetLocalH()); ipmap[i].SetSize(num_points); - ipmap[i] = 0; // PointIndex::INVALID; + ipmap[i] = state0; // 0; // PointIndex::INVALID; m.SetDimension( mesh.GetDimension() ); m.SetGeometry( mesh.GetGeometry() ); @@ -86,8 +91,8 @@ namespace netgen { if(seg.domin > 0 && seg.domin == seg.domout) { - ipmap[seg.domin-1][seg[0]] = 1; - ipmap[seg.domin-1][seg[1]] = 1; + ipmap[seg.domin-1][seg[0]] = state1; // 1; + ipmap[seg.domin-1][seg[1]] = state1; // 1; } } @@ -105,7 +110,7 @@ namespace netgen auto & sels = ret[dom-1].mesh->SurfaceElements(); for(auto pi : sel.PNums()) - ipmap[dom-1][pi] = 1; + ipmap[dom-1][pi] = state1; // 1; sels.Append(sel); } } @@ -114,17 +119,17 @@ namespace netgen for(const auto & el : mesh.VolumeElements()) { auto dom = el.GetIndex(); - + auto & els = ret[dom-1].mesh->VolumeElements(); for(auto pi : el.PNums()) - ipmap[dom-1][pi] = 1; + ipmap[dom-1][pi] = state1; // 1; els.Append(el); } // mark locked/fixed points for each domain TODO: domain bounding box to add only relevant points? for(auto pi : mesh.LockedPoints()) for(auto i : Range(ret)) - ipmap[i][pi] = 2; + ipmap[i][pi] = state2; // 2; // add used points to domain mesh, build point mapping for(auto i : Range(ret)) @@ -133,11 +138,11 @@ namespace netgen auto & pmap = ret[i].pmap; for(auto pi : Range(ipmap[i])) - if(ipmap[i][pi]) + if(ipmap[i][pi] != state0) { const auto& mp = mesh[pi]; auto pi_new = m.AddPoint( mp, mp.GetLayer(), mp.Type() ); - if(ipmap[i][pi] == 2) + if(ipmap[i][pi] == state2) // 2) mesh.AddLockedPoint(pi_new); ipmap[i][pi] = pi_new; pmap.Append( pi ); diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 7034f6ce..de93b229 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -223,18 +223,18 @@ namespace netgen friend constexpr netgen::PointIndex ngcore::IndexBASE (); template friend class PointIndices; + friend constexpr PointIndex operator+ (PointIndex, int); + friend constexpr PointIndex operator+ (PointIndex, size_t); + friend constexpr PointIndex operator+ (int, PointIndex); + friend constexpr PointIndex operator+ (size_t, PointIndex); + friend constexpr PointIndex operator- (PointIndex, int); + friend constexpr int operator- (PointIndex, PointIndex); /* - friend PointIndex operator+ (PointIndex, int); - friend PointIndex operator+ (PointIndex, size_t); - friend PointIndex operator+ (int, PointIndex); - friend PointIndex operator+ (size_t, PointIndex); - friend PointIndex operator- (PointIndex, int); - friend int operator- (PointIndex, PointIndex); 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); - friend bool operator!= (PointIndex a, PointIndex b); + friend constexpr bool operator== (PointIndex a, PointIndex b); + friend constexpr bool operator!= (PointIndex a, PointIndex b); */ public: @@ -261,18 +261,19 @@ namespace netgen void DoArchive (Archive & ar) { ar & i; } }; + constexpr inline PointIndex operator+ (PointIndex pi, int i) { return PointIndex(pi.i+i); } + constexpr inline PointIndex operator+ (PointIndex pi, size_t i) { return PointIndex(pi.i+i); } + constexpr inline PointIndex operator+ (int i, PointIndex pi) { return PointIndex(pi.i+i); } + constexpr inline PointIndex operator+ (size_t i, PointIndex pi) { return PointIndex(pi.i+i); } + constexpr inline PointIndex operator- (PointIndex pi, int i) { return PointIndex(pi.i-i); } + constexpr inline int operator- (PointIndex pa, PointIndex pb) { return PointIndex(pa.i-pb.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+ (int i, PointIndex pi) { return PointIndex(pi.i+i); } - inline PointIndex operator+ (size_t i, PointIndex pi) { return PointIndex(pi.i+i); } - inline PointIndex operator- (PointIndex pi, int i) { return PointIndex(pi.i-i); } - inline int operator- (PointIndex pa, PointIndex pb) { return PointIndex(pa.i-pb.i); } inline bool operator< (PointIndex a, PointIndex b) { return a.i < b.i; } inline bool operator> (PointIndex a, PointIndex b) { return a.i > b.i; } inline bool operator>= (PointIndex a, PointIndex b) { return a.i >= b.i; } - inline bool operator== (PointIndex a, PointIndex b) { return a.i == b.i; } - inline bool operator!= (PointIndex a, PointIndex b) { return a.i != b.i; } + inline constexpr bool operator== (PointIndex a, PointIndex b) { return a.i == b.i; } + inline constexpr bool operator!= (PointIndex a, PointIndex b) { return a.i != b.i; } */ } @@ -316,6 +317,12 @@ namespace netgen constexpr PointIndices (PointIndex i1, PointIndex i2) : INDEX_2(i1,i2) { ; } PointIndex operator[] (int i) const { return PointIndex(INDEX_2::operator[](i)); } PointIndex & operator[] (int i) { return reinterpret_cast(INDEX_2::operator[](i)); } + + PointIndex & I1 () { return (*this)[0]; } + PointIndex & I2 () { return (*this)[1]; } + PointIndex I1 () const { return (*this)[0]; } + PointIndex I2 () const { return (*this)[1]; } + using INDEX_2::Sort; static PointIndices Sort(PointIndex i1, PointIndex i2) { return INDEX_2::Sort(i1, i2); } template @@ -333,6 +340,13 @@ namespace netgen constexpr 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(INDEX_3::operator[](i)); } + PointIndex & I1 () { return (*this)[0]; } + PointIndex & I2 () { return (*this)[1]; } + PointIndex & I3 () { return (*this)[2]; } + PointIndex I1 () const { return (*this)[0]; } + PointIndex I2 () const { return (*this)[1]; } + PointIndex I3 () const { return (*this)[2]; } + using INDEX_3::Sort; static PointIndices Sort(PointIndex i1, PointIndex i2, PointIndex i3) { return INDEX_3::Sort(i1, i2, i3); } template @@ -347,6 +361,16 @@ namespace netgen 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)); } + + PointIndex & I1 () { return (*this)[0]; } + PointIndex & I2 () { return (*this)[1]; } + PointIndex & I3 () { return (*this)[2]; } + PointIndex & I4 () { return (*this)[3]; } + PointIndex I1 () const { return (*this)[0]; } + PointIndex I2 () const { return (*this)[1]; } + PointIndex I3 () const { return (*this)[2]; } + PointIndex I4 () const { return (*this)[3]; } + using INDEX_4::Sort; // static PointIndices Sort(PointIndex i1, PointIndex i2, PointIndex i3, PointIndex i4) { return INDEX_4::Sort(i1, i2, i3, i4); } template @@ -1646,9 +1670,9 @@ namespace netgen /// int haltnode; /// - int haltsegmentp1; + PointIndex haltsegmentp1; /// - int haltsegmentp2; + PointIndex haltsegmentp2; /// int haltexistingline; /// diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp index 767e8ca7..826e2f7d 100644 --- a/libsrc/meshing/parallelmesh.cpp +++ b/libsrc/meshing/parallelmesh.cpp @@ -329,7 +329,9 @@ namespace netgen /** The same table as per_verts, but TRANSITIVE!! **/ auto iterate_per_verts_trans = [&](auto f){ NgArray allvs; - for (int k = PointIndex::BASE; k < GetNV()+PointIndex::BASE; k++) + // for (int k = PointIndex::BASE; k < GetNV()+PointIndex::BASE; k++) + for (PointIndex k = IndexBASE(); + k < GetNV()+IndexBASE(); k++) { allvs.SetSize(0); allvs.Append(per_verts[k]); diff --git a/libsrc/meshing/parser3.cpp b/libsrc/meshing/parser3.cpp index 6c2ce74a..0eeee374 100644 --- a/libsrc/meshing/parser3.cpp +++ b/libsrc/meshing/parser3.cpp @@ -879,16 +879,16 @@ void vnetrule :: LoadRule (istream & ist) // NgArray & freeset = *freesets.Get(fs); NgArray & freesetedges = *freeedges.Last(); NgArray & freesetfaces = *freefaces.Get(fs); - int k,l; - INDEX ind; + // int k,l; + // INDEX ind; - for (k = 1; k <= freesetfaces.Size(); k++) + for (int k = 1; k <= freesetfaces.Size(); k++) { // threeint tr = freesetfaces.Get(k); - for (l = k+1; l <= freesetfaces.Size(); l++) + for (int l = k+1; l <= freesetfaces.Size(); l++) { - ind = NeighbourTrianglePoint(freesetfaces.Get(k), freesetfaces.Get(l)); + INDEX ind = NeighbourTrianglePoint(freesetfaces.Get(k), freesetfaces.Get(l)); if (!ind) continue; INDEX_3 f1(freesetfaces.Get(k).i1, @@ -897,7 +897,7 @@ void vnetrule :: LoadRule (istream & ist) INDEX_3 f2(freesetfaces.Get(l).i1, freesetfaces.Get(l).i2, freesetfaces.Get(l).i3); - INDEX_2 ed(0, 0); + PointIndices<2> ed(PointIndex::INVALID, PointIndex::INVALID); for (int f11 = 1; f11 <= 3; f11++) for (int f12 = 1; f12 <= 3; f12++) if (f11 != f12) @@ -916,8 +916,8 @@ void vnetrule :: LoadRule (istream & ist) { for (int elr = 1; elr <= 4; elr++) { - if (GetPointNrMod (eli, elr) == ed.I(1) && - GetPointNrMod (eli, elr+2) == ed.I(2)) + if (GetPointNrMod (eli, elr) == ed[0] && + GetPointNrMod (eli, elr+2) == ed[1]) { /* (*testout) << "ed is diagonal of rectangle" << endl; diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp index 8d3ef551..4e5f5446 100644 --- a/libsrc/meshing/smoothing3.cpp +++ b/libsrc/meshing/smoothing3.cpp @@ -1731,26 +1731,31 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp, pf.SetPointIndex (pi); - PointIndex brother (-1); + constexpr PointIndex state0(PointIndex::INVALID); + constexpr PointIndex statem1 = state0-1; + + PointIndex brother = statem1; if(usesum) { - for(int j=0; brother == -1 && jSize(); j++) + for(int j=0; brother == statem1 && jSize(); j++) { if(pi < (*used_idmaps)[j]->Size() + IndexBASE()) { brother = (*(*used_idmaps)[j])[pi]; - if(brother == pi || brother == 0) - brother = -1; + if(brother == pi || brother == state0) + brother = statem1; } } - if(brother >= pi) + // if(brother >= pi) + if(brother-pi >= 0) { pf2ptr->SetPointIndex(brother); pf2ptr->SetNV(*nv[brother-1]); } } - if(usesum && brother < pi) + // if(usesum && brother < pi) + if(usesum && (brother-pi < 0)) continue; //pf.UnSetNV(); x = 0; @@ -1759,12 +1764,12 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp, pf.SetNV(*nv[pi-1]); x = 0; - int pok = (brother == -1) ? (pf.Func (x) < 1e10) : (pf_sum.Func (x) < 1e10); + int pok = (brother == statem1) ? (pf.Func (x) < 1e10) : (pf_sum.Func (x) < 1e10); if (pok) { - if(brother == -1) + if(brother == statem1) BFGS (x, pf, par); else BFGS (x, pf_sum, par); @@ -1773,7 +1778,7 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp, for(int j=0; j<3; j++) points[pi](j) += x(j);// - scal*nv[i-1].X(j); - if(brother != -1) + if(brother != statem1) for(int j=0; j<3; j++) points[brother](j) += x(j);// - scal*nv[brother-1].X(j); @@ -1783,8 +1788,8 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp, { cout << "el not ok" << endl; (*testout) << "el not ok" << endl - << " func " << ((brother == -1) ? pf.Func(x) : pf_sum.Func (x)) << endl; - if(brother != -1) + << " func " << ((brother == statem1) ? pf.Func(x) : pf_sum.Func (x)) << endl; + if(brother != statem1) (*testout) << " func1 " << pf.Func(x) << endl << " func2 " << pf2ptr->Func(x) << endl; } diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp index bf953140..fc178b9f 100644 --- a/libsrc/meshing/topology.cpp +++ b/libsrc/meshing/topology.cpp @@ -102,11 +102,11 @@ namespace netgen auto eledges = MeshTopology::GetEdges (el.GetType()); for (int k = 0; k < eledges.Size(); k++) { - INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); + PointIndices<2> edge(el[eledges[k][0]], el[eledges[k][1]]); - int edgedir = (edge.I1() > edge.I2()); - if (edgedir) swap (edge.I1(), edge.I2()); - if (edge.I1() != v) continue; + bool edgedir = edge[0] > edge[1]; + if (edgedir) swap (edge[0], edge[1]); + if (edge[0] != v) continue; func (edge, elnr, k, 3); } @@ -119,7 +119,7 @@ namespace netgen auto eledges = MeshTopology::GetEdges (el.GetType()); for (int k = 0; k < eledges.Size(); k++) { - INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); + PointIndices<2> edge(el[eledges[k][0]], el[eledges[k][1]]); int edgedir = (edge.I1() > edge.I2()); if (edgedir) swap (edge.I1(), edge.I2()); @@ -133,7 +133,7 @@ namespace netgen for (SegmentIndex elnr : top.GetVertexSegments(v)) { const Segment & el = mesh[elnr]; - INDEX_2 edge(el[0], el[1]); + PointIndices<2> edge(el[0], el[1]); int edgedir = (edge.I1() > edge.I2()); if (edgedir) swap (edge.I1(), edge.I2()); @@ -159,8 +159,8 @@ namespace netgen if (elfaces[j][3] < 0) { // triangle - INDEX_4 face(el[elfaces[j][0]], el[elfaces[j][1]], - el[elfaces[j][2]], 0); + PointIndices<4> face(el[elfaces[j][0]], el[elfaces[j][1]], + el[elfaces[j][2]], 0); [[maybe_unused]] int facedir = 0; @@ -197,8 +197,8 @@ namespace netgen { // quad // int facenum; - INDEX_4 face4(el[elfaces[j][0]], el[elfaces[j][1]], - el[elfaces[j][2]], el[elfaces[j][3]]); + PointIndices<4> face4(el[elfaces[j][0]], el[elfaces[j][1]], + el[elfaces[j][2]], el[elfaces[j][3]]); // int facedir = 0; if (min2 (face4.I1(), face4.I2()) > @@ -264,24 +264,24 @@ namespace netgen // int facenum; // int facedir; - INDEX_4 face(el.PNum(elfaces[0][0]), - el.PNum(elfaces[0][1]), - el.PNum(elfaces[0][2]),0); + PointIndices<4> face(el.PNum(elfaces[0][0]), + el.PNum(elfaces[0][1]), + el.PNum(elfaces[0][2]),0); // facedir = 0; - if (face.I1() > face.I2()) + if (face[0] > face[1]) { - swap (face.I1(), face.I2()); + swap (face[0], face[1]); // facedir += 1; } - if (face.I2() > face.I3()) + if (face[1] > face[2]) { - swap (face.I2(), face.I3()); + swap (face[1], face[2]); // facedir += 2; } - if (face.I1() > face.I2()) + if (face.I1() > face[1]) { - swap (face.I1(), face.I2()); + swap (face.I1(), face[1]); // facedir += 4; } @@ -313,10 +313,10 @@ namespace netgen // int facenum; // int facedir; - INDEX_4 face4(el.PNum(elfaces[0][0]), - el.PNum(elfaces[0][1]), - el.PNum(elfaces[0][2]), - el.PNum(elfaces[0][3])); + PointIndices<4> face4(el.PNum(elfaces[0][0]), + el.PNum(elfaces[0][1]), + el.PNum(elfaces[0][2]), + el.PNum(elfaces[0][3])); // facedir = 0; if (min2 (face4.I1(), face4.I2()) > @@ -698,7 +698,7 @@ namespace netgen // edge is splitting edge in middle of triangle: for (int j = 1; j <= 2; j++) { - IVec<2> paedge1, paedge2, paedge3; + IVec<2,PointIndex> paedge1, paedge2, paedge3; int orient_inner = 0; if (j == 1) { @@ -722,7 +722,7 @@ namespace netgen Swap (paedge3[0], paedge3[1]); // if first vertex number is -1, then don't try to find entry in node2edge hash table - if ( paedge1[0] == PointIndex::BASE-1 || paedge2[0] == PointIndex::BASE-1 ) + if ( !paedge1[0].IsValid() || !paedge2[0].IsValid() ) continue; int paedgenr1=-1, paedgenr2=-1, paedgenr3=-1, orient1 = 0, orient2 = 0; @@ -751,21 +751,21 @@ namespace netgen if (!bisect_edge) // not a bisect edge (then a red edge) { - IVec<2> paedge1, paedge2, paedge3; + IVec<2,PointIndex> paedge1, paedge2, paedge3; int orient1 = 0, orient2 = 0, orient3=0; // int orient_inner = 0; - paedge1 = IVec<2> (pa0[0], pa0[1]); - paedge2 = IVec<2> (pa1[0], pa1[1]); + paedge1 = IVec<2,PointIndex> (pa0[0], pa0[1]); + paedge2 = IVec<2,PointIndex> (pa1[0], pa1[1]); // find common vertex and the third pa edge if (pa0[0]==pa1[0]){// 00 //orient1 = 0; orient2 = 1; if (pa0[1] (pa0[1], pa1[1]); + paedge3 = IVec<2,PointIndex> (pa0[1], pa1[1]); }else{ //orient3 = 0; - paedge3 = IVec<2> (pa1[1], pa0[1]); + paedge3 = IVec<2,PointIndex> (pa1[1], pa0[1]); } } else if (pa0[0]==pa1[1]){//01 @@ -773,10 +773,10 @@ namespace netgen //orient2 = 0; if (pa0[1] (pa0[1], pa1[0]); + paedge3 = IVec<2,PointIndex> (pa0[1], pa1[0]); }else{ //orient3 = 0; - paedge3 = IVec<2> (pa1[0], pa0[1]); + paedge3 = IVec<2,PointIndex> (pa1[0], pa0[1]); } } else if (pa0[1]==pa1[0]){//10 @@ -784,10 +784,10 @@ namespace netgen orient2 = 1; if (pa0[0] (pa0[0], pa1[1]); + paedge3 = IVec<2,PointIndex> (pa0[0], pa1[1]); }else{ //orient3 = 0; - paedge3 = IVec<2> (pa1[1], pa0[0]); + paedge3 = IVec<2,PointIndex> (pa1[1], pa0[0]); } } else if (pa0[1]==pa1[1]){//11 @@ -795,10 +795,10 @@ namespace netgen //orient2 = 0; if (pa0[0] (pa0[0], pa1[0]); + paedge3 = IVec<2,PointIndex> (pa0[0], pa1[0]); }else{ //orient3 = 0; - paedge3 = IVec<2> (pa1[0], pa0[0]); + paedge3 = IVec<2,PointIndex> (pa1[0], pa0[0]); } } @@ -1836,7 +1836,7 @@ namespace netgen for(auto & face : elfaces) { auto v = face2vert[face-1]; - if(v[3]!=0) + if(v[3].IsValid()) cerr << "GetElementFaces with orientation currently not supported for quads" << endl; int classnr = 0; @@ -2270,12 +2270,13 @@ namespace netgen void MeshTopology :: GetFaceEdges (int fnr, NgArray & fedges, bool withorientation) const { - NgArrayMem pi(4); + // NgArrayMem pi(4); // NgArrayMem eledges; fedges.SetSize (0); - GetFaceVertices(fnr, pi); - + // GetFaceVertices(fnr, pi); + auto pi = GetFaceVertices(fnr-1); + // Sort Edges according to global vertex numbers // e1 = fmax, f2 // e2 = fmax, f1 @@ -2345,8 +2346,8 @@ namespace netgen // fedges.Append (eledges[j]); for(int k=0;k 0) diff --git a/libsrc/meshing/topology.hpp b/libsrc/meshing/topology.hpp index 761158cf..b8a00268 100644 --- a/libsrc/meshing/topology.hpp +++ b/libsrc/meshing/topology.hpp @@ -207,7 +207,7 @@ public: FlatArray GetVertexPointElements (PointIndex vnr) const { return vert2pointelement[vnr]; } - DLL_HEADER int GetVerticesEdge ( int v1, int v2) const; + DLL_HEADER int GetVerticesEdge ( PointIndex v1, PointIndex v2) const; void GetSegmentVolumeElements ( int segnr, NgArray & els ) const; void GetSegmentSurfaceElements ( int segnr, NgArray & els ) const; diff --git a/libsrc/visualization/vssolution.cpp b/libsrc/visualization/vssolution.cpp index 74501c78..4634e46e 100644 --- a/libsrc/visualization/vssolution.cpp +++ b/libsrc/visualization/vssolution.cpp @@ -248,7 +248,7 @@ namespace netgen const Element2d & el = (*mesh)[sei]; surf_ost << el.GetNP(); for (int j = 0; j < el.GetNP(); j++) - surf_ost << " " << el[j] - PointIndex::BASE; + surf_ost << " " << el[j] - IndexBASE(); surf_ost << "\n"; } surf_ost << "\nCELL_TYPES " << mesh->GetNSE() << "\n"; @@ -291,7 +291,7 @@ namespace netgen const Element & el = (*mesh)[ei]; ost << el.GetNP(); for (int j = 0; j < el.GetNP(); j++) - ost << " " << el[j] - PointIndex::BASE; + ost << " " << el[j] - IndexBASE(); ost << "\n"; } ost << "\nCELL_TYPES " << mesh->GetNE() << "\n";