ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-01-27 05:10:34 +05:00
Code Issues Packages Projects Releases Wiki Activity

operators +/- for PointIndex

Browse Source
This commit is contained in:
Joachim Schoeberl 2024-12-28 21:26:05 +01:00
parent 00e3a3490b
commit 75032f9905
14 changed files with 155 additions and 116 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
libsrc/geom2d/genmesh2d.cpp
View File

@ -139,7 +139,9 @@ namespace netgen
mark = spline.GetPoint (edgelength); 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 mark3(mark(0), mark(1), 0);
Point3d oldmark3(oldmark(0), oldmark(1), 0); Point3d oldmark3(oldmark(0), oldmark(1), 0);
@ -157,12 +159,12 @@ namespace netgen
if ( mesh[PointIndex(locsearch[k])].GetLayer() == spline.layer) if ( mesh[PointIndex(locsearch[k])].GetLayer() == spline.layer)
pi2 = locsearch[k]; pi2 = locsearch[k];
if (pi1 == -1) if (!pi1.IsValid())
{ {
pi1 = mesh.AddPoint(oldmark3, spline.layer); pi1 = mesh.AddPoint(oldmark3, spline.layer);
searchtree.Insert (oldmark3, pi1); searchtree.Insert (oldmark3, pi1);
} }
if (pi2 == -1) if (!pi2.IsValid())
{ {
pi2 = mesh.AddPoint(mark3, spline.layer); pi2 = mesh.AddPoint(mark3, spline.layer);
searchtree.Insert (mark3, pi2); searchtree.Insert (mark3, pi2);
@ -292,13 +294,13 @@ namespace netgen
Point<2> hnewp = (j == 1) ? splines[i]->StartPI() : splines[i]->EndPI(); Point<2> hnewp = (j == 1) ? splines[i]->StartPI() : splines[i]->EndPI();
Point<3> newp(hnewp(0), hnewp(1), 0); Point<3> newp(hnewp(0), hnewp(1), 0);
int layer = GetSpline(i).layer; int layer = GetSpline(i).layer;
int npi = -1; PointIndex npi(PointIndex::INVALID);
for (PointIndex pi = PointIndex::BASE; for (PointIndex pi = IndexBASE<PointIndex>();
pi < mesh2d.GetNP()+PointIndex::BASE; pi++) pi < mesh2d.GetNP()+IndexBASE<PointIndex>(); pi++)
if (Dist2 (mesh2d.Point(pi), newp) < 1e-12 * diam2 && mesh2d.Point(pi).GetLayer() == layer) if (Dist2 (mesh2d.Point(pi), newp) < 1e-12 * diam2 && mesh2d.Point(pi).GetLayer() == layer)
npi = pi; npi = pi;
if (npi == -1) if (!npi.IsValid())
{ {
npi = mesh2d.AddPoint (newp, layer); npi = mesh2d.AddPoint (newp, layer);
searchtree.Insert (newp, npi); searchtree.Insert (newp, npi);

2
libsrc/meshing/bisect.cpp
View File

@ -547,7 +547,7 @@ namespace netgen
e2[0] = (*idmaps[k])[e1[0]]; e2[0] = (*idmaps[k])[e1[0]];
e2[1] = (*idmaps[k])[e1[1]]; 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()) e1.I1() == e2.I1() || e1.I2() == e2.I2())
continue; continue;

10
libsrc/meshing/boundarylayer.cpp
View File

@ -87,7 +87,7 @@ Vec<3> CalcGrowthVector (FlatArray<Vec<3>> ns)
auto [maxpos1, maxpos2] = FindCloseVectors(ns); auto [maxpos1, maxpos2] = FindCloseVectors(ns);
Array<Vec<3>> new_normals; Array<Vec<3>> new_normals;
new_normals = ns; 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]); auto average = 0.5 * (ns[maxpos1] + ns[maxpos2]);
average.Normalize(); average.Normalize();
new_normals[maxpos1] = average; new_normals[maxpos1] = average;
@ -154,7 +154,7 @@ Vec<3> BoundaryLayerTool ::getEdgeTangent(PointIndex pi, int edgenr, FlatArray<S
auto& seg = *p_seg; auto& seg = *p_seg;
if (seg.edgenr != edgenr) if (seg.edgenr != edgenr)
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);
} }
@ -721,7 +721,7 @@ void BoundaryLayerTool ::InsertNewElements(
auto g1 = getGroups(p1, segj.si); auto g1 = getGroups(p1, segj.si);
if (g0.Size() == 1 && g1.Size() == 1) if (g0.Size() == 1 && g1.Size() == 1)
auto s = [[maybe_unused]] auto s =
addSegment(newPoint(p0, -1, g0[0]), newPoint(p1, -1, g1[0])); addSegment(newPoint(p0, -1, g0[0]), newPoint(p1, -1, g1[0]));
else else
{ {
@ -798,7 +798,7 @@ void BoundaryLayerTool ::InsertNewElements(
s3[0] = p3; s3[0] = p3;
s3[1] = p4; s3[1] = p4;
s3[2] = PointIndex::INVALID; s3[2] = PointIndex::INVALID;
auto pair = p3 < p4 ? make_pair(p3, p4) : make_pair(p4, p3); // auto pair = p3 < p4 ? make_pair(p3, p4) : make_pair(p4, p3);
s3.edgenr = getEdgeNr(segj.edgenr); s3.edgenr = getEdgeNr(segj.edgenr);
s3.si = segj.si; s3.si = segj.si;
new_segments.Append(s3); new_segments.Append(s3);
@ -853,7 +853,7 @@ void BoundaryLayerTool ::InsertNewElements(
if (groups.Size() == 1) if (groups.Size() == 1)
return groups[0]; return groups[0];
auto& growth_groups = special_boundary_points[pi].growth_groups; // auto& growth_groups = special_boundary_points[pi].growth_groups;
auto vdir = Center(mesh[sel[0]], mesh[sel[1]], mesh[sel[2]]) - mesh[pi]; auto vdir = Center(mesh[sel[0]], mesh[sel[1]], mesh[sel[2]]) - mesh[pi];
auto dot = vdir * special_boundary_points[pi].separating_direction; auto dot = vdir * special_boundary_points[pi].separating_direction;

6
libsrc/meshing/boundarylayer2d.cpp
View File

@ -329,7 +329,7 @@ namespace netgen
auto current_si = si; auto current_si = si;
auto first = current_seg[0]; auto first = current_seg[0];
auto current = -1; PointIndex current(PointIndex::INVALID);
auto next = current_seg[1]; auto next = current_seg[1];
if(points_done.Test(first)) if(points_done.Test(first))
@ -354,7 +354,7 @@ namespace netgen
current_si = sj; current_si = sj;
current_seg = mesh[sj]; current_seg = mesh[sj];
next = current_seg[0] + current_seg[1] - current; next = current_seg[0]-current + current_seg[1];
break; break;
} }
} }
@ -493,7 +493,7 @@ namespace netgen
if(growthvectors[pi].Length2() == 0.0) if(growthvectors[pi].Length2() == 0.0)
continue; continue;
PointIndex pi1 = seg0[0] + seg0[1] - pi; PointIndex pi1 = seg0[0] - pi + seg0[1];
auto p1 = mesh[pi1]; auto p1 = mesh[pi1];
auto p = mesh[pi]; auto p = mesh[pi];

2
libsrc/meshing/boundarylayer_interpolate.cpp
View File

@ -283,7 +283,7 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors()
auto np = mesh.GetNP(); auto np = mesh.GetNP();
auto hasMoved = [&] (PointIndex pi) { auto hasMoved = [&] (PointIndex pi) {
return (pi - PointIndex::BASE >= np_old) || mapto[pi].Size() > 0 || special_boundary_points.count(pi); return (pi - IndexBASE<PointIndex>() >= np_old) || mapto[pi].Size() > 0 || special_boundary_points.count(pi);
}; };
std::set<PointIndex> points_set; std::set<PointIndex> points_set;

12
libsrc/meshing/boundarylayer_limiter.hpp
View File

@ -288,7 +288,7 @@ struct GrowthVectorLimiter
return; return;
for (PointIndex pi : IntRange(tool.np, mesh.GetNP())) for (PointIndex pi : IntRange(tool.np, mesh.GetNP()))
{ {
auto pi_from = map_from[pi]; // auto pi_from = map_from[pi];
std::set<PointIndex> pis; std::set<PointIndex> pis;
for (auto sei : p2sel[pi]) for (auto sei : p2sel[pi])
for (auto pi_ : tool.new_sels[sei].PNums()) for (auto pi_ : tool.new_sels[sei].PNums())
@ -360,7 +360,7 @@ struct GrowthVectorLimiter
if (sel.GetNP() == 4) if (sel.GetNP() == 4)
continue; continue;
const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex()); // const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex());
auto np = sel.GetNP(); auto np = sel.GetNP();
double shift = 1.0; double shift = 1.0;
@ -447,7 +447,7 @@ struct GrowthVectorLimiter
for (auto sei : SurfaceElementsRange()) for (auto sei : SurfaceElementsRange())
{ {
const auto& sel = Get(sei); const auto& sel = Get(sei);
auto sel_index = sel.GetIndex(); // auto sel_index = sel.GetIndex();
Box<3> box(Box<3>::EMPTY_BOX); Box<3> box(Box<3>::EMPTY_BOX);
for (auto pi : sel.PNums()) for (auto pi : sel.PNums())
@ -466,7 +466,7 @@ struct GrowthVectorLimiter
RegionTimer rt(t); RegionTimer rt(t);
BuildSearchTree(trig_shift); BuildSearchTree(trig_shift);
auto np_new = mesh.Points().Size(); auto np_new = mesh.Points().Size();
int counter = 0; // int counter = 0;
for (auto i : IntRange(tool.np, np_new)) for (auto i : IntRange(tool.np, np_new))
{ {
PointIndex pi_to = i + PointIndex::BASE; PointIndex pi_to = i + PointIndex::BASE;
@ -478,7 +478,7 @@ struct GrowthVectorLimiter
continue; continue;
Box<3> box(Box<3>::EMPTY_BOX); 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, 0));
box.Add(GetPoint(pi_to, GetLimit(pi_from))); box.Add(GetPoint(pi_to, GetLimit(pi_from)));
@ -488,7 +488,7 @@ struct GrowthVectorLimiter
return false; return false;
if (sel.PNums().Contains(pi_to)) if (sel.PNums().Contains(pi_to))
return false; return false;
counter++; // counter++;
f(pi_to, sei); f(pi_to, sei);
return false; return false;
}); });

23
libsrc/meshing/meshfunc.cpp
View File

@ -59,6 +59,11 @@ namespace netgen
auto num_points = mesh.GetNP(); auto num_points = mesh.GetNP();
auto num_facedescriptors = mesh.GetNFD(); auto num_facedescriptors = mesh.GetNFD();
constexpr PointIndex state0 = IndexBASE<PointIndex>()-1;
constexpr PointIndex state1 = state0+1;
constexpr PointIndex state2 = state0+2;
for(auto i : Range(ret)) for(auto i : Range(ret))
{ {
auto & md = ret[i]; auto & md = ret[i];
@ -73,7 +78,7 @@ namespace netgen
m.SetLocalH(mesh.GetLocalH()); m.SetLocalH(mesh.GetLocalH());
ipmap[i].SetSize(num_points); ipmap[i].SetSize(num_points);
ipmap[i] = 0; // PointIndex::INVALID; ipmap[i] = state0; // 0; // PointIndex::INVALID;
m.SetDimension( mesh.GetDimension() ); m.SetDimension( mesh.GetDimension() );
m.SetGeometry( mesh.GetGeometry() ); m.SetGeometry( mesh.GetGeometry() );
@ -86,8 +91,8 @@ namespace netgen
{ {
if(seg.domin > 0 && seg.domin == seg.domout) if(seg.domin > 0 && seg.domin == seg.domout)
{ {
ipmap[seg.domin-1][seg[0]] = 1; ipmap[seg.domin-1][seg[0]] = state1; // 1;
ipmap[seg.domin-1][seg[1]] = 1; ipmap[seg.domin-1][seg[1]] = state1; // 1;
} }
} }
@ -105,7 +110,7 @@ namespace netgen
auto & sels = ret[dom-1].mesh->SurfaceElements(); auto & sels = ret[dom-1].mesh->SurfaceElements();
for(auto pi : sel.PNums()) for(auto pi : sel.PNums())
ipmap[dom-1][pi] = 1; ipmap[dom-1][pi] = state1; // 1;
sels.Append(sel); sels.Append(sel);
} }
} }
@ -114,17 +119,17 @@ namespace netgen
for(const auto & el : mesh.VolumeElements()) for(const auto & el : mesh.VolumeElements())
{ {
auto dom = el.GetIndex(); auto dom = el.GetIndex();
auto & els = ret[dom-1].mesh->VolumeElements(); auto & els = ret[dom-1].mesh->VolumeElements();
for(auto pi : el.PNums()) for(auto pi : el.PNums())
ipmap[dom-1][pi] = 1; ipmap[dom-1][pi] = state1; // 1;
els.Append(el); els.Append(el);
} }
// mark locked/fixed points for each domain TODO: domain bounding box to add only relevant points? // mark locked/fixed points for each domain TODO: domain bounding box to add only relevant points?
for(auto pi : mesh.LockedPoints()) for(auto pi : mesh.LockedPoints())
for(auto i : Range(ret)) for(auto i : Range(ret))
ipmap[i][pi] = 2; ipmap[i][pi] = state2; // 2;
// add used points to domain mesh, build point mapping // add used points to domain mesh, build point mapping
for(auto i : Range(ret)) for(auto i : Range(ret))
@ -133,11 +138,11 @@ namespace netgen
auto & pmap = ret[i].pmap; auto & pmap = ret[i].pmap;
for(auto pi : Range(ipmap[i])) for(auto pi : Range(ipmap[i]))
if(ipmap[i][pi]) if(ipmap[i][pi] != state0)
{ {
const auto& mp = mesh[pi]; const auto& mp = mesh[pi];
auto pi_new = m.AddPoint( mp, mp.GetLayer(), mp.Type() ); 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); mesh.AddLockedPoint(pi_new);
ipmap[i][pi] = pi_new; ipmap[i][pi] = pi_new;
pmap.Append( pi ); pmap.Append( pi );

60
libsrc/meshing/meshtype.hpp
View File

@ -223,18 +223,18 @@ namespace netgen
friend constexpr netgen::PointIndex ngcore::IndexBASE<netgen::PointIndex> (); friend constexpr netgen::PointIndex ngcore::IndexBASE<netgen::PointIndex> ();
template <int N> friend class PointIndices; template <int N> 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 bool operator>= (PointIndex a, PointIndex b);
friend bool operator== (PointIndex a, PointIndex b); friend constexpr bool operator== (PointIndex a, PointIndex b);
friend bool operator!= (PointIndex a, PointIndex b); friend constexpr bool operator!= (PointIndex a, PointIndex b);
*/ */
public: public:
@ -261,18 +261,19 @@ namespace netgen
void DoArchive (Archive & ar) { ar & i; } 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 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 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) { ; } 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) const { return PointIndex(INDEX_2::operator[](i)); }
PointIndex & operator[] (int i) { return reinterpret_cast<PointIndex&>(INDEX_2::operator[](i)); } PointIndex & operator[] (int i) { return reinterpret_cast<PointIndex&>(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; using INDEX_2::Sort;
static PointIndices Sort(PointIndex i1, PointIndex i2) { return INDEX_2::Sort(i1, i2); } static PointIndices Sort(PointIndex i1, PointIndex i2) { return INDEX_2::Sort(i1, i2); }
template <size_t J> template <size_t J>
@ -333,6 +340,13 @@ namespace netgen
constexpr PointIndices (PointIndex i1, PointIndex i2, PointIndex i3) : INDEX_3(i1,i2,i3) { ; } 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) const { return PointIndex(INDEX_3::operator[](i)); }
PointIndex & operator[] (int i) { return reinterpret_cast<PointIndex&>(INDEX_3::operator[](i)); } PointIndex & operator[] (int i) { return reinterpret_cast<PointIndex&>(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; using INDEX_3::Sort;
static PointIndices Sort(PointIndex i1, PointIndex i2, PointIndex i3) { return INDEX_3::Sort(i1, i2, i3); } static PointIndices Sort(PointIndex i1, PointIndex i2, PointIndex i3) { return INDEX_3::Sort(i1, i2, i3); }
template <size_t J> template <size_t J>
@ -347,6 +361,16 @@ namespace netgen
PointIndices (PointIndex i1, PointIndex i2, PointIndex i3, PointIndex i4) : INDEX_4(i1,i2,i3,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) const { return PointIndex(INDEX_4::operator[](i)); }
PointIndex & operator[] (int i) { return reinterpret_cast<PointIndex&>(INDEX_4::operator[](i)); } PointIndex & operator[] (int i) { return reinterpret_cast<PointIndex&>(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; using INDEX_4::Sort;
// static PointIndices Sort(PointIndex i1, PointIndex i2, PointIndex i3, PointIndex i4) { return INDEX_4::Sort(i1, i2, i3, i4); } // static PointIndices Sort(PointIndex i1, PointIndex i2, PointIndex i3, PointIndex i4) { return INDEX_4::Sort(i1, i2, i3, i4); }
template <size_t J> template <size_t J>
@ -1646,9 +1670,9 @@ namespace netgen
/// ///
int haltnode; int haltnode;
/// ///
int haltsegmentp1; PointIndex haltsegmentp1;
/// ///
int haltsegmentp2; PointIndex haltsegmentp2;
/// ///
int haltexistingline; int haltexistingline;
/// ///

4
libsrc/meshing/parallelmesh.cpp
View File

@ -329,7 +329,9 @@ namespace netgen
/** The same table as per_verts, but TRANSITIVE!! **/ /** The same table as per_verts, but TRANSITIVE!! **/
auto iterate_per_verts_trans = [&](auto f){ auto iterate_per_verts_trans = [&](auto f){
NgArray<int> allvs; NgArray<int> 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<PointIndex>();
k < GetNV()+IndexBASE<PointIndex>(); k++)
{ {
allvs.SetSize(0); allvs.SetSize(0);
allvs.Append(per_verts[k]); allvs.Append(per_verts[k]);

16
libsrc/meshing/parser3.cpp
View File

@ -879,16 +879,16 @@ void vnetrule :: LoadRule (istream & ist)
// NgArray<int> & freeset = *freesets.Get(fs); // NgArray<int> & freeset = *freesets.Get(fs);
NgArray<twoint> & freesetedges = *freeedges.Last(); NgArray<twoint> & freesetedges = *freeedges.Last();
NgArray<threeint> & freesetfaces = *freefaces.Get(fs); NgArray<threeint> & freesetfaces = *freefaces.Get(fs);
int k,l; // int k,l;
INDEX ind; // INDEX ind;
for (k = 1; k <= freesetfaces.Size(); k++) for (int k = 1; k <= freesetfaces.Size(); k++)
{ {
// threeint tr = freesetfaces.Get(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; if (!ind) continue;
INDEX_3 f1(freesetfaces.Get(k).i1, INDEX_3 f1(freesetfaces.Get(k).i1,
@ -897,7 +897,7 @@ void vnetrule :: LoadRule (istream & ist)
INDEX_3 f2(freesetfaces.Get(l).i1, INDEX_3 f2(freesetfaces.Get(l).i1,
freesetfaces.Get(l).i2, freesetfaces.Get(l).i2,
freesetfaces.Get(l).i3); 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 f11 = 1; f11 <= 3; f11++)
for (int f12 = 1; f12 <= 3; f12++) for (int f12 = 1; f12 <= 3; f12++)
if (f11 != f12) if (f11 != f12)
@ -916,8 +916,8 @@ void vnetrule :: LoadRule (istream & ist)
{ {
for (int elr = 1; elr <= 4; elr++) for (int elr = 1; elr <= 4; elr++)
{ {
if (GetPointNrMod (eli, elr) == ed.I(1) && if (GetPointNrMod (eli, elr) == ed[0] &&
GetPointNrMod (eli, elr+2) == ed.I(2)) GetPointNrMod (eli, elr+2) == ed[1])
{ {
/* /*
(*testout) << "ed is diagonal of rectangle" << endl; (*testout) << "ed is diagonal of rectangle" << endl;

27
libsrc/meshing/smoothing3.cpp
View File

@ -1731,26 +1731,31 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp,
pf.SetPointIndex (pi); pf.SetPointIndex (pi);
PointIndex brother (-1); constexpr PointIndex state0(PointIndex::INVALID);
constexpr PointIndex statem1 = state0-1;
PointIndex brother = statem1;
if(usesum) if(usesum)
{ {
for(int j=0; brother == -1 && j<used_idmaps->Size(); j++) for(int j=0; brother == statem1 && j<used_idmaps->Size(); j++)
{ {
if(pi < (*used_idmaps)[j]->Size() + IndexBASE<PointIndex>()) if(pi < (*used_idmaps)[j]->Size() + IndexBASE<PointIndex>())
{ {
brother = (*(*used_idmaps)[j])[pi]; brother = (*(*used_idmaps)[j])[pi];
if(brother == pi || brother == 0) if(brother == pi || brother == state0)
brother = -1; brother = statem1;
} }
} }
if(brother >= pi) // if(brother >= pi)
if(brother-pi >= 0)
{ {
pf2ptr->SetPointIndex(brother); pf2ptr->SetPointIndex(brother);
pf2ptr->SetNV(*nv[brother-1]); pf2ptr->SetNV(*nv[brother-1]);
} }
} }
if(usesum && brother < pi) // if(usesum && brother < pi)
if(usesum && (brother-pi < 0))
continue; continue;
//pf.UnSetNV(); x = 0; //pf.UnSetNV(); x = 0;
@ -1759,12 +1764,12 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp,
pf.SetNV(*nv[pi-1]); pf.SetNV(*nv[pi-1]);
x = 0; 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 (pok)
{ {
if(brother == -1) if(brother == statem1)
BFGS (x, pf, par); BFGS (x, pf, par);
else else
BFGS (x, pf_sum, par); BFGS (x, pf_sum, par);
@ -1773,7 +1778,7 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp,
for(int j=0; j<3; j++) for(int j=0; j<3; j++)
points[pi](j) += x(j);// - scal*nv[i-1].X(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++) for(int j=0; j<3; j++)
points[brother](j) += x(j);// - scal*nv[brother-1].X(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; cout << "el not ok" << endl;
(*testout) << "el not ok" << endl (*testout) << "el not ok" << endl
<< " func " << ((brother == -1) ? pf.Func(x) : pf_sum.Func (x)) << endl; << " func " << ((brother == statem1) ? pf.Func(x) : pf_sum.Func (x)) << endl;
if(brother != -1) if(brother != statem1)
(*testout) << " func1 " << pf.Func(x) << endl (*testout) << " func1 " << pf.Func(x) << endl
<< " func2 " << pf2ptr->Func(x) << endl; << " func2 " << pf2ptr->Func(x) << endl;
} }

87
libsrc/meshing/topology.cpp
View File

@ -102,11 +102,11 @@ namespace netgen
auto eledges = MeshTopology::GetEdges (el.GetType()); auto eledges = MeshTopology::GetEdges (el.GetType());
for (int k = 0; k < eledges.Size(); k++) 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()); bool edgedir = edge[0] > edge[1];
if (edgedir) swap (edge.I1(), edge.I2()); if (edgedir) swap (edge[0], edge[1]);
if (edge.I1() != v) continue; if (edge[0] != v) continue;
func (edge, elnr, k, 3); func (edge, elnr, k, 3);
} }
@ -119,7 +119,7 @@ namespace netgen
auto eledges = MeshTopology::GetEdges (el.GetType()); auto eledges = MeshTopology::GetEdges (el.GetType());
for (int k = 0; k < eledges.Size(); k++) 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()); int edgedir = (edge.I1() > edge.I2());
if (edgedir) swap (edge.I1(), edge.I2()); if (edgedir) swap (edge.I1(), edge.I2());
@ -133,7 +133,7 @@ namespace netgen
for (SegmentIndex elnr : top.GetVertexSegments(v)) for (SegmentIndex elnr : top.GetVertexSegments(v))
{ {
const Segment & el = mesh[elnr]; 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()); int edgedir = (edge.I1() > edge.I2());
if (edgedir) swap (edge.I1(), edge.I2()); if (edgedir) swap (edge.I1(), edge.I2());
@ -159,8 +159,8 @@ namespace netgen
if (elfaces[j][3] < 0) if (elfaces[j][3] < 0)
{ // triangle { // triangle
INDEX_4 face(el[elfaces[j][0]], el[elfaces[j][1]], PointIndices<4> face(el[elfaces[j][0]], el[elfaces[j][1]],
el[elfaces[j][2]], 0); el[elfaces[j][2]], 0);
[[maybe_unused]] int facedir = 0; [[maybe_unused]] int facedir = 0;
@ -197,8 +197,8 @@ namespace netgen
{ {
// quad // quad
// int facenum; // int facenum;
INDEX_4 face4(el[elfaces[j][0]], el[elfaces[j][1]], PointIndices<4> face4(el[elfaces[j][0]], el[elfaces[j][1]],
el[elfaces[j][2]], el[elfaces[j][3]]); el[elfaces[j][2]], el[elfaces[j][3]]);
// int facedir = 0; // int facedir = 0;
if (min2 (face4.I1(), face4.I2()) > if (min2 (face4.I1(), face4.I2()) >
@ -264,24 +264,24 @@ namespace netgen
// int facenum; // int facenum;
// int facedir; // int facedir;
INDEX_4 face(el.PNum(elfaces[0][0]), PointIndices<4> face(el.PNum(elfaces[0][0]),
el.PNum(elfaces[0][1]), el.PNum(elfaces[0][1]),
el.PNum(elfaces[0][2]),0); el.PNum(elfaces[0][2]),0);
// facedir = 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; // facedir += 1;
} }
if (face.I2() > face.I3()) if (face[1] > face[2])
{ {
swap (face.I2(), face.I3()); swap (face[1], face[2]);
// facedir += 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; // facedir += 4;
} }
@ -313,10 +313,10 @@ namespace netgen
// int facenum; // int facenum;
// int facedir; // int facedir;
INDEX_4 face4(el.PNum(elfaces[0][0]), PointIndices<4> face4(el.PNum(elfaces[0][0]),
el.PNum(elfaces[0][1]), el.PNum(elfaces[0][1]),
el.PNum(elfaces[0][2]), el.PNum(elfaces[0][2]),
el.PNum(elfaces[0][3])); el.PNum(elfaces[0][3]));
// facedir = 0; // facedir = 0;
if (min2 (face4.I1(), face4.I2()) > if (min2 (face4.I1(), face4.I2()) >
@ -698,7 +698,7 @@ namespace netgen
// edge is splitting edge in middle of triangle: // edge is splitting edge in middle of triangle:
for (int j = 1; j <= 2; j++) for (int j = 1; j <= 2; j++)
{ {
IVec<2> paedge1, paedge2, paedge3; IVec<2,PointIndex> paedge1, paedge2, paedge3;
int orient_inner = 0; int orient_inner = 0;
if (j == 1) if (j == 1)
{ {
@ -722,7 +722,7 @@ namespace netgen
Swap (paedge3[0], paedge3[1]); Swap (paedge3[0], paedge3[1]);
// if first vertex number is -1, then don't try to find entry in node2edge hash table // 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; continue;
int paedgenr1=-1, paedgenr2=-1, paedgenr3=-1, orient1 = 0, orient2 = 0; 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) 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 orient1 = 0, orient2 = 0, orient3=0;
// int orient_inner = 0; // int orient_inner = 0;
paedge1 = IVec<2> (pa0[0], pa0[1]); paedge1 = IVec<2,PointIndex> (pa0[0], pa0[1]);
paedge2 = IVec<2> (pa1[0], pa1[1]); paedge2 = IVec<2,PointIndex> (pa1[0], pa1[1]);
// find common vertex and the third pa edge // find common vertex and the third pa edge
if (pa0[0]==pa1[0]){// 00 if (pa0[0]==pa1[0]){// 00
//orient1 = 0; //orient1 = 0;
orient2 = 1; orient2 = 1;
if (pa0[1]<pa1[1]){ if (pa0[1]<pa1[1]){
orient3 = 1; orient3 = 1;
paedge3 = IVec<2> (pa0[1], pa1[1]); paedge3 = IVec<2,PointIndex> (pa0[1], pa1[1]);
}else{ }else{
//orient3 = 0; //orient3 = 0;
paedge3 = IVec<2> (pa1[1], pa0[1]); paedge3 = IVec<2,PointIndex> (pa1[1], pa0[1]);
} }
} }
else if (pa0[0]==pa1[1]){//01 else if (pa0[0]==pa1[1]){//01
@ -773,10 +773,10 @@ namespace netgen
//orient2 = 0; //orient2 = 0;
if (pa0[1]<pa1[0]){ if (pa0[1]<pa1[0]){
orient3 = 1; orient3 = 1;
paedge3 = IVec<2> (pa0[1], pa1[0]); paedge3 = IVec<2,PointIndex> (pa0[1], pa1[0]);
}else{ }else{
//orient3 = 0; //orient3 = 0;
paedge3 = IVec<2> (pa1[0], pa0[1]); paedge3 = IVec<2,PointIndex> (pa1[0], pa0[1]);
} }
} }
else if (pa0[1]==pa1[0]){//10 else if (pa0[1]==pa1[0]){//10
@ -784,10 +784,10 @@ namespace netgen
orient2 = 1; orient2 = 1;
if (pa0[0]<pa1[1]){ if (pa0[0]<pa1[1]){
orient3 = 1; orient3 = 1;
paedge3 = IVec<2> (pa0[0], pa1[1]); paedge3 = IVec<2,PointIndex> (pa0[0], pa1[1]);
}else{ }else{
//orient3 = 0; //orient3 = 0;
paedge3 = IVec<2> (pa1[1], pa0[0]); paedge3 = IVec<2,PointIndex> (pa1[1], pa0[0]);
} }
} }
else if (pa0[1]==pa1[1]){//11 else if (pa0[1]==pa1[1]){//11
@ -795,10 +795,10 @@ namespace netgen
//orient2 = 0; //orient2 = 0;
if (pa0[0]<pa1[0]){ if (pa0[0]<pa1[0]){
orient3 = 1; orient3 = 1;
paedge3 = IVec<2> (pa0[0], pa1[0]); paedge3 = IVec<2,PointIndex> (pa0[0], pa1[0]);
}else{ }else{
//orient3 = 0; //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) for(auto & face : elfaces)
{ {
auto v = face2vert[face-1]; auto v = face2vert[face-1];
if(v[3]!=0) if(v[3].IsValid())
cerr << "GetElementFaces with orientation currently not supported for quads" << endl; cerr << "GetElementFaces with orientation currently not supported for quads" << endl;
int classnr = 0; int classnr = 0;
@ -2270,12 +2270,13 @@ namespace netgen
void MeshTopology :: GetFaceEdges (int fnr, NgArray<int> & fedges, bool withorientation) const void MeshTopology :: GetFaceEdges (int fnr, NgArray<int> & fedges, bool withorientation) const
{ {
NgArrayMem<int,4> pi(4); // NgArrayMem<int,4> pi(4);
// NgArrayMem<int,12> eledges; // NgArrayMem<int,12> eledges;
fedges.SetSize (0); fedges.SetSize (0);
GetFaceVertices(fnr, pi); // GetFaceVertices(fnr, pi);
auto pi = GetFaceVertices(fnr-1);
// Sort Edges according to global vertex numbers // Sort Edges according to global vertex numbers
// e1 = fmax, f2 // e1 = fmax, f2
// e2 = fmax, f1 // e2 = fmax, f1
@ -2345,8 +2346,8 @@ namespace netgen
// fedges.Append (eledges[j]); // fedges.Append (eledges[j]);
for(int k=0;k<nfa_ref_edges;k++) for(int k=0;k<nfa_ref_edges;k++)
{ {
int w1 = el[ref_faces[fa][fa_ref_edges[k][0]-1]-1]; PointIndex w1 = el[ref_faces[fa][fa_ref_edges[k][0]-1]-1];
int w2 = el[ref_faces[fa][fa_ref_edges[k][1]-1]-1]; PointIndex w2 = el[ref_faces[fa][fa_ref_edges[k][1]-1]-1];
if(withorientation) if(withorientation)
{ {
@ -2392,7 +2393,7 @@ namespace netgen
} }
int MeshTopology :: GetVerticesEdge ( int v1, int v2 ) const int MeshTopology :: GetVerticesEdge ( PointIndex v1, PointIndex v2 ) const
{ {
/* /*
if (vert2element.Size() > 0) if (vert2element.Size() > 0)

2
libsrc/meshing/topology.hpp
View File

@ -207,7 +207,7 @@ public:
FlatArray<int> GetVertexPointElements (PointIndex vnr) const FlatArray<int> GetVertexPointElements (PointIndex vnr) const
{ return vert2pointelement[vnr]; } { 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<ElementIndex> & els ) const; void GetSegmentVolumeElements ( int segnr, NgArray<ElementIndex> & els ) const;
void GetSegmentSurfaceElements ( int segnr, NgArray<SurfaceElementIndex> & els ) const; void GetSegmentSurfaceElements ( int segnr, NgArray<SurfaceElementIndex> & els ) const;

4
libsrc/visualization/vssolution.cpp
View File

@ -248,7 +248,7 @@ namespace netgen
const Element2d & el = (*mesh)[sei]; const Element2d & el = (*mesh)[sei];
surf_ost << el.GetNP(); surf_ost << el.GetNP();
for (int j = 0; j < el.GetNP(); j++) for (int j = 0; j < el.GetNP(); j++)
surf_ost << " " << el[j] - PointIndex::BASE; surf_ost << " " << el[j] - IndexBASE<PointIndex>();
surf_ost << "\n"; surf_ost << "\n";
} }
surf_ost << "\nCELL_TYPES " << mesh->GetNSE() << "\n"; surf_ost << "\nCELL_TYPES " << mesh->GetNSE() << "\n";
@ -291,7 +291,7 @@ namespace netgen
const Element & el = (*mesh)[ei]; const Element & el = (*mesh)[ei];
ost << el.GetNP(); ost << el.GetNP();
for (int j = 0; j < el.GetNP(); j++) for (int j = 0; j < el.GetNP(); j++)
ost << " " << el[j] - PointIndex::BASE; ost << " " << el[j] - IndexBASE<PointIndex>();
ost << "\n"; ost << "\n";
} }
ost << "\nCELL_TYPES " << mesh->GetNE() << "\n"; ost << "\nCELL_TYPES " << mesh->GetNE() << "\n";