diff --git a/libsrc/core/array.hpp b/libsrc/core/array.hpp index b001acec..4418ef6c 100644 --- a/libsrc/core/array.hpp +++ b/libsrc/core/array.hpp @@ -554,6 +554,13 @@ namespace ngcore // { return CArray (data+pos); } NETGEN_INLINE T * operator+ (size_t pos) const { return data+pos; } + /// access first element. check by macro NETGEN_CHECK_RANGE + T & First () const + { + NETGEN_CHECK_RANGE(0,0,size); + return data[0]; + } + /// access last element. check by macro NETGEN_CHECK_RANGE T & Last () const { diff --git a/libsrc/meshing/improve2gen.cpp b/libsrc/meshing/improve2gen.cpp index 75e3fb14..db7485fc 100644 --- a/libsrc/meshing/improve2gen.cpp +++ b/libsrc/meshing/improve2gen.cpp @@ -11,9 +11,9 @@ namespace netgen public: NgArray oldels; NgArray newels; - NgArray deledges; - NgArray incelsonnode; - NgArray reused; + NgArray> deledges; + Array incelsonnode; + Array reused; int bonus; int onp; }; @@ -55,128 +55,137 @@ namespace netgen NgArray rules; NgArray elmap; NgArray elrot; - NgArray pmap; - NgArray pgi; + Array pmap; + Array pgi; int surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr(); ImprovementRule * r1; + PointIndex p1 = IndexBASE(); + PointIndex p2 = p1+1; + PointIndex p3 = p2+1; + PointIndex p4 = p3+1; + PointIndex p5 = p4+1; + PointIndex p6 = p5+1; + PointIndex p7 = p6+1; + + // 2 triangles to quad r1 = new ImprovementRule; - r1->oldels.Append (Element2d (1, 2, 3)); - r1->oldels.Append (Element2d (3, 2, 4)); - r1->newels.Append (Element2d (1, 2, 4, 3)); - r1->deledges.Append (INDEX_2 (2,3)); + r1->oldels.Append (Element2d (p1, p2, p3)); + r1->oldels.Append (Element2d (p3, p2, p4)); + r1->newels.Append (Element2d (p1, p2, p4, p3)); + r1->deledges.Append ( { p2, p3 } ); r1->onp = 4; r1->bonus = 2; rules.Append (r1); // 2 quad to 1 quad r1 = new ImprovementRule; - r1->oldels.Append (Element2d (1, 2, 3, 4)); - r1->oldels.Append (Element2d (4, 3, 2, 5)); - r1->newels.Append (Element2d (1, 2, 5, 4)); - r1->deledges.Append (INDEX_2 (2, 3)); - r1->deledges.Append (INDEX_2 (3, 4)); + r1->oldels.Append (Element2d (p1, p2, p3, p4)); + r1->oldels.Append (Element2d (p4, p3, p2, p5)); + r1->newels.Append (Element2d (p1, p2, p5, p4)); + r1->deledges.Append ( { p2, p3 } ); + r1->deledges.Append ( { p3, p4 } ); r1->onp = 5; r1->bonus = 0; rules.Append (r1); // swap quads r1 = new ImprovementRule; - r1->oldels.Append (Element2d (1, 2, 3, 4)); - r1->oldels.Append (Element2d (3, 2, 5, 6)); - r1->newels.Append (Element2d (1, 6, 3, 4)); - r1->newels.Append (Element2d (1, 2, 5, 6)); - r1->deledges.Append (INDEX_2 (2, 3)); + r1->oldels.Append (Element2d (p1, p2, p3, p4)); + r1->oldels.Append (Element2d (p3, p2, p5, p6)); + r1->newels.Append (Element2d (p1, p6, p3, p4)); + r1->newels.Append (Element2d (p1, p2, p5, p6)); + r1->deledges.Append ( { p2, p3 } ); r1->onp = 6; r1->bonus = 0; rules.Append (r1); // three quads to 2 r1 = new ImprovementRule; - r1->oldels.Append (Element2d (1, 2, 3, 4)); - r1->oldels.Append (Element2d (2, 5, 6, 3)); - r1->oldels.Append (Element2d (3, 6, 7, 4)); - r1->newels.Append (Element2d (1, 2, 5, 4)); - r1->newels.Append (Element2d (4, 5, 6, 7)); - r1->deledges.Append (INDEX_2 (2, 3)); - r1->deledges.Append (INDEX_2 (3, 4)); - r1->deledges.Append (INDEX_2 (3, 6)); + r1->oldels.Append (Element2d (p1, p2, p3, p4)); + r1->oldels.Append (Element2d (p2, p5, p6, p3)); + r1->oldels.Append (Element2d (p3, p6, p7, p4)); + r1->newels.Append (Element2d (p1, p2, p5, p4)); + r1->newels.Append (Element2d (p4, p5, p6, p7)); + r1->deledges.Append ( { p2, p3 } ); + r1->deledges.Append ( { p3, p4 } ); + r1->deledges.Append ( { p3, p6 } ); r1->onp = 7; r1->bonus = -1; rules.Append (r1); // quad + 2 connected trigs to quad r1 = new ImprovementRule; - r1->oldels.Append (Element2d (1, 2, 3, 4)); - r1->oldels.Append (Element2d (2, 5, 3)); - r1->oldels.Append (Element2d (3, 5, 4)); - r1->newels.Append (Element2d (1, 2, 5, 4)); - r1->deledges.Append (INDEX_2 (2, 3)); - r1->deledges.Append (INDEX_2 (3, 4)); - r1->deledges.Append (INDEX_2 (3, 5)); + r1->oldels.Append (Element2d (p1, p2, p3, p4)); + r1->oldels.Append (Element2d (p2, p5, p3)); + r1->oldels.Append (Element2d (p3, p5, p4)); + r1->newels.Append (Element2d (p1, p2, p5, p4)); + r1->deledges.Append ( { p2, p3 } ); + r1->deledges.Append ( { p3, p4 } ); + r1->deledges.Append ( { p3, p5 } ); r1->onp = 5; r1->bonus = 0; rules.Append (r1); // quad + 2 non-connected trigs to quad (a and b) r1 = new ImprovementRule; - r1->oldels.Append (Element2d (1, 2, 3, 4)); - r1->oldels.Append (Element2d (2, 6, 3)); - r1->oldels.Append (Element2d (1, 4, 5)); - r1->newels.Append (Element2d (1, 3, 4, 5)); - r1->newels.Append (Element2d (1, 2, 6, 3)); - r1->deledges.Append (INDEX_2 (1, 4)); - r1->deledges.Append (INDEX_2 (2, 3)); + r1->oldels.Append (Element2d (p1, p2, p3, p4)); + r1->oldels.Append (Element2d (p2, p6, p3)); + r1->oldels.Append (Element2d (p1, p4, p5)); + r1->newels.Append (Element2d (p1, p3, p4, p5)); + r1->newels.Append (Element2d (p1, p2, p6, p3)); + r1->deledges.Append ( { p1, p4 } ); + r1->deledges.Append ( { p2, p3 } ); r1->onp = 6; r1->bonus = 0; rules.Append (r1); r1 = new ImprovementRule; - r1->oldels.Append (Element2d (1, 2, 3, 4)); - r1->oldels.Append (Element2d (2, 6, 3)); - r1->oldels.Append (Element2d (1, 4, 5)); - r1->newels.Append (Element2d (1, 2, 4, 5)); - r1->newels.Append (Element2d (4, 2, 6, 3)); - r1->deledges.Append (INDEX_2 (1, 4)); - r1->deledges.Append (INDEX_2 (2, 3)); + r1->oldels.Append (Element2d (p1, p2, p3, p4)); + r1->oldels.Append (Element2d (p2, p6, p3)); + r1->oldels.Append (Element2d (p1, p4, p5)); + r1->newels.Append (Element2d (p1, p2, p4, p5)); + r1->newels.Append (Element2d (p4, p2, p6, p3)); + r1->deledges.Append ( { p1, p4 } ); + r1->deledges.Append ( { p2, p3 } ); r1->onp = 6; r1->bonus = 0; rules.Append (r1); // two quad + trig -> one quad + trig r1 = new ImprovementRule; - r1->oldels.Append (Element2d (1, 2, 3, 4)); - r1->oldels.Append (Element2d (2, 5, 6, 3)); - r1->oldels.Append (Element2d (4, 3, 6)); - r1->newels.Append (Element2d (1, 2, 6, 4)); - r1->newels.Append (Element2d (2, 5, 6)); - r1->deledges.Append (INDEX_2 (2, 3)); - r1->deledges.Append (INDEX_2 (3, 4)); - r1->deledges.Append (INDEX_2 (3, 6)); + r1->oldels.Append (Element2d (p1, p2, p3, p4)); + r1->oldels.Append (Element2d (p2, p5, p6, p3)); + r1->oldels.Append (Element2d (p4, p3, p6)); + r1->newels.Append (Element2d (p1, p2, p6, p4)); + r1->newels.Append (Element2d (p2, p5, p6)); + r1->deledges.Append ( { p2, p3 } ); + r1->deledges.Append ( { p3, p4 } ); + r1->deledges.Append ( { p3, p6 } ); r1->onp = 6; r1->bonus = -1; rules.Append (r1); // swap quad + trig (a and b) r1 = new ImprovementRule; - r1->oldels.Append (Element2d (1, 2, 3, 4)); - r1->oldels.Append (Element2d (2, 5, 3)); - r1->newels.Append (Element2d (2, 5, 3, 4)); - r1->newels.Append (Element2d (1, 2, 4)); - r1->deledges.Append (INDEX_2 (2, 3)); + r1->oldels.Append (Element2d (p1, p2, p3, p4)); + r1->oldels.Append (Element2d (p2, p5, p3)); + r1->newels.Append (Element2d (p2, p5, p3, p4)); + r1->newels.Append (Element2d (p1, p2, p4)); + r1->deledges.Append ( { p2, p3 } ); r1->onp = 5; r1->bonus = 0; rules.Append (r1); r1 = new ImprovementRule; - r1->oldels.Append (Element2d (1, 2, 3, 4)); - r1->oldels.Append (Element2d (2, 5, 3)); - r1->newels.Append (Element2d (1, 2, 5, 3)); - r1->newels.Append (Element2d (1, 3, 4)); - r1->deledges.Append (INDEX_2 (2, 3)); + r1->oldels.Append (Element2d (p1, p2, p3, p4)); + r1->oldels.Append (Element2d (p2, p5, p3)); + r1->newels.Append (Element2d (p1, p2, p5, p3)); + r1->newels.Append (Element2d (p1, p3, p4)); + r1->deledges.Append ( { p2, p3 } ); r1->onp = 5; r1->bonus = 0; rules.Append (r1); @@ -184,12 +193,12 @@ namespace netgen // 2 quads to quad + 2 trigs r1 = new ImprovementRule; - r1->oldels.Append (Element2d (1, 2, 3, 4)); - r1->oldels.Append (Element2d (3, 2, 5, 6)); - r1->newels.Append (Element2d (1, 5, 6, 4)); - r1->newels.Append (Element2d (1, 2, 5)); - r1->newels.Append (Element2d (4, 6, 3)); - r1->deledges.Append (INDEX_2 (2, 3)); + r1->oldels.Append (Element2d (p1, p2, p3, p4)); + r1->oldels.Append (Element2d (p3, p2, p5, p6)); + r1->newels.Append (Element2d (p1, p5, p6, p4)); + r1->newels.Append (Element2d (p1, p2, p5)); + r1->newels.Append (Element2d (p4, p6, p3)); + r1->deledges.Append ( { p2, p3 } ); r1->onp = 6; r1->bonus = 0; // rules.Append (r1); @@ -203,18 +212,21 @@ namespace netgen mapped = 0; - for (int ri = 0; ri < rules.Size(); ri++) { ImprovementRule & rule = *rules[ri]; rule.incelsonnode.SetSize (rule.onp); rule.reused.SetSize (rule.onp); + /* for (int j = 0; j < rule.onp; j++) { rule.incelsonnode[j] = 0; rule.reused[j] = 0; } + */ + rule.incelsonnode = 0; + rule.reused = 0; /* for (int j = 1; j <= rule.oldels.Size(); j++) @@ -226,15 +238,15 @@ namespace netgen */ for (const Element2d & el : rule.oldels) for (PointIndex pi : el.PNums()) - rule.incelsonnode[pi-IndexBASE()]--; + rule.incelsonnode[pi]--; for (int j = 1; j <= rule.newels.Size(); j++) { const Element2d & el = rule.newels.Elem(j); for (int k = 1; k <= el.GetNP(); k++) { - rule.incelsonnode.Elem(el.PNum(k))++; - rule.reused.Elem(el.PNum(k)) = 1; + rule.incelsonnode[el.PNum(k)]++; + rule.reused[el.PNum(k)] = 1; } } } @@ -242,8 +254,8 @@ namespace netgen - TABLE elonnode(np); - NgArray nelonnode(np); + DynamicTable elonnode(np); + Array nelonnode(np); TABLE nbels(ne); nelonnode = -4; @@ -314,13 +326,13 @@ namespace netgen if (el0.IsDeleted()) continue; if (el0.GetNP() != rel0.GetNP()) continue; - - pmap = PointIndex (-1); + // pmap = PointIndex (-1); + for (auto & p : pmap) p.Invalidate(); for (int k = 0; k < el0.GetNP(); k++) { - pmap.Elem(rel0[k]) = el0.PNumMod(k+elrot[0]+1); - pgi.Elem(rel0[k]) = el0.GeomInfoPiMod(k+elrot[0]+1); + pmap[rel0[k]] = el0.PNumMod(k+elrot[0]+1); + pgi[rel0[k]] = el0.GeomInfoPiMod(k+elrot[0]+1); } ok = 1; @@ -342,16 +354,16 @@ namespace netgen possible = 1; for (int k = 0; k < rel.GetNP(); k++) - if (pmap.Elem(rel[k]).IsValid() && - pmap.Elem(rel[k]) != el.PNumMod (k+elrot[i]+1)) + if (pmap[rel[k]].IsValid() && + pmap[rel[k]] != el.PNumMod (k+elrot[i]+1)) possible = 0; if (possible) { for (int k = 0; k < el.GetNP(); k++) { - pmap.Elem(rel[k]) = el.PNumMod(k+elrot[i]+1); - pgi.Elem(rel[k]) = el.GeomInfoPiMod(k+elrot[i]+1); + pmap[rel[k]] = el.PNumMod(k+elrot[i]+1); + pgi[rel[k]] = el.GeomInfoPiMod(k+elrot[i]+1); } break; } @@ -370,8 +382,8 @@ namespace netgen for(int i=0; ok && i olddef) continue; @@ -398,7 +410,8 @@ namespace netgen // calc metric badness double bad1 = 0, bad2 = 0; // SelectSurfaceOfPoint (mesh.Point(pmap.Get(1)), pgi.Get(1)); - auto n = geo.GetNormal(surfnr, mesh.Point(pmap.Get(1)), &pgi.Elem(1)); + // auto n = geo.GetNormal(surfnr, mesh.Point(pmap.Get(1)), &pgi.Elem(1)); + auto n = geo.GetNormal(surfnr, mesh.Point(pmap.First()), pgi.Data()); for (int j = 0; j < rule.oldels.Size(); j++) bad1 += mesh[elmap[j]].CalcJacobianBadness (mesh.Points(), n); @@ -409,7 +422,7 @@ namespace netgen const Element2d & rnel = rule.newels.Get(j); Element2d nel(rnel.GetNP()); for (int k = 1; k <= rnel.GetNP(); k++) - nel.PNum(k) = pmap.Get(rnel.PNum(k)); + nel.PNum(k) = pmap[rnel.PNum(k)]; bad2 += nel.CalcJacobianBadness (mesh.Points(), n); } @@ -427,8 +440,8 @@ namespace netgen nel.SetIndex (faceindex); for (int k = 1; k <= rnel.GetNP(); k++) { - nel.PNum(k) = pmap.Get(rnel.PNum(k)); - nel.GeomInfoPi(k) = pgi.Get(rnel.PNum(k)); + nel.PNum(k) = pmap[rnel.PNum(k)]; + nel.GeomInfoPi(k) = pgi[rnel.PNum(k)]; } mesh.AddSurfaceElement(nel); @@ -437,8 +450,8 @@ namespace netgen for (int j = 0; j < rule.oldels.Size(); j++) mesh.Delete (elmap[j]); - for (int j = 1; j <= pmap.Size(); j++) - nelonnode[pmap.Get(j)] += rule.incelsonnode.Get(j); + for (PointIndex j : pmap.Range()) + nelonnode[pmap[j]] += rule.incelsonnode[j]; used[ri]++; } diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 2b64bea7..15d71354 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -73,7 +73,7 @@ namespace netgen m.SetLocalH(mesh.GetLocalH()); ipmap[i].SetSize(num_points); - ipmap[i] = PointIndex::INVALID; + ipmap[i] = 0; // PointIndex::INVALID; m.SetDimension( mesh.GetDimension() ); m.SetGeometry( mesh.GetGeometry() ); diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 97848e86..2ed0e481 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -219,10 +219,11 @@ namespace netgen // throw Exception("illegal PointIndex, use PointIndex::INVALID instead"); #endif } + /* friend constexpr netgen::PointIndex ngcore::IndexBASE (); - friend istream & operator>> (istream &, PointIndex &); - friend ostream & operator<< (ostream &, const PointIndex &); + // friend istream & operator>> (istream &, PointIndex &); + // friend ostream & operator<< (ostream &, const PointIndex &); template friend class PointIndices; friend PointIndex operator+ (PointIndex, int); friend PointIndex operator+ (PointIndex, size_t); @@ -259,7 +260,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); } @@ -284,15 +284,21 @@ namespace ngcore namespace netgen { - + // input-output is 1-based inline istream & operator>> (istream & ist, PointIndex & pi) { - int i; ist >> i; pi = PointIndex(i); return ist; + // int i; ist >> i; pi = PointIndex(i); return ist; + int i; ist >> i; + // pi = PointIndex(i); + pi = IndexBASE()+i-1; + return ist; } inline ostream & operator<< (ostream & ost, const PointIndex & pi) { - return (ost << int(pi)); + int intpi = pi - IndexBASE() + 1; + return (ost << intpi); + // return (ost << int(pi)); } template class PointIndices;