diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index 5f63d35e..cbd2a73b 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -6,7 +6,6 @@ namespace netgen { - static const int deltetfaces[][3] = { { 1, 2, 3 }, { 2, 0, 3 }, @@ -14,10 +13,6 @@ namespace netgen { 1, 0, 2 } }; - - - - class DelaunayTet { PointIndex pnums[4]; @@ -41,9 +36,6 @@ namespace netgen PointIndex & operator[] (int i) { return pnums[i]; } PointIndex operator[] (int i) const { return pnums[i]; } - // int & NB1(int i) { return nb[i-1]; } - // int NB1(int i) const { return nb[i-1]; } - int & NB(int i) { return nb[i]; } int NB(int i) const { return nb[i]; } @@ -58,29 +50,6 @@ namespace netgen return 3; } - /* - void GetFace1 (int i, INDEX_3 & face) const - { - face.I(1) = pnums[deltetfaces[i-1][0]]; - face.I(2) = pnums[deltetfaces[i-1][1]]; - face.I(3) = pnums[deltetfaces[i-1][2]]; - } - */ - void GetFace (int i, INDEX_3 & face) const - { - face.I(1) = pnums[deltetfaces[i][0]]; - face.I(2) = pnums[deltetfaces[i][1]]; - face.I(3) = pnums[deltetfaces[i][2]]; - } - - /* - INDEX_3 GetFace1 (int i) const - { - return INDEX_3 (pnums[deltetfaces[i-1][0]], - pnums[deltetfaces[i-1][1]], - pnums[deltetfaces[i-1][2]]); - } - */ INDEX_3 GetFace (int i) const { return INDEX_3 (pnums[deltetfaces[i][0]], @@ -88,15 +57,6 @@ namespace netgen pnums[deltetfaces[i][2]]); } - /* - void GetFace1 (int i, Element2d & face) const - { - // face.SetType(TRIG); - face[0] = pnums[deltetfaces[i-1][0]]; - face[1] = pnums[deltetfaces[i-1][1]]; - face[2] = pnums[deltetfaces[i-1][2]]; - } - */ void GetFace (int i, Element2d & face) const { // face.SetType(TRIG); @@ -104,7 +64,6 @@ namespace netgen face[1] = pnums[deltetfaces[i][1]]; face[2] = pnums[deltetfaces[i][2]]; } - }; @@ -113,8 +72,6 @@ namespace netgen - - /* Table to maintain neighbour elements */ @@ -147,7 +104,7 @@ namespace netgen // get neighbour of element elnr in direction fnr int GetNB (int elnr, int fnr) { - return tets.Get(elnr).NB(fnr-1); + return tets.Get(elnr).NB(fnr); } // @@ -193,7 +150,6 @@ namespace netgen - /* connected lists of cosphereical elements */ @@ -265,24 +221,23 @@ namespace netgen Array & freelist, SphereList & list, IndexSet & insphere, IndexSet & closesphere) { + // static Timer t("Meshing3::AddDelaunayPoint"); RegionTimer reg(t); /* find any sphere, such that newp is contained in */ - DelaunayTet el; int cfelind = -1; const Point<3> * pp[4]; Point<3> pc; - double r2; Point3d tpmin, tpmax; tettree.GetIntersecting (newp, newp, treesearch); + double quot,minquot(1e20); - for (int j = 0; j < treesearch.Size(); j++) + for (auto jjj : treesearch) { - int jjj = treesearch[j]; quot = Dist2 (centers.Get(jjj), newp) / radi2.Get(jjj); if((cfelind == -1 || quot < 0.99*minquot) && quot < 1) @@ -296,46 +251,12 @@ namespace netgen } - /* - int i, j, k, l; - if (!felind) - { - cerr << "not in any sphere, 1" << endl; - // old, non tree search - - double mindist = 1e10; - for (j = 1; j <= tempels.Size(); j++) - { - if (tempels.Get(j).PNum(1)) - { - double toofar = - Dist2 (centers.Get(j), newp) - radi2.Get(j); - if (toofar < mindist || toofar < 1e-7) - { - mindist = toofar; - cout << " dist2 = " << Dist2 (centers.Get(j), newp) - << " radi2 = " << radi2.Get(j) << endl; - } - if (toofar < 0) - { - el = tempels.Get(j); - felind = j; - cout << "sphere found !" << endl; - break; - } - } - } - cout << "point is too far from sheres: " << mindist << endl; - } - */ - if (cfelind == -1) { PrintWarning ("Delaunay, point not in any sphere"); return; } - /* insphere: point is in sphere -> delete element closesphere: point is close to sphere -> considered for same center @@ -388,7 +309,7 @@ namespace netgen // check neighbour-tets for (int j = starti; j < nstarti; j++) - for (int k = 1; k <= 4; k++) + for (int k = 0; k < 4; k++) { int helind = insphere.GetArray().Get(j); int nbind = meshnb.GetNB (helind, k); @@ -416,61 +337,48 @@ namespace netgen } else { - /* - Element2d face; - tempels.Get(helind).GetFace (k, face); + INDEX_3 i3 = tempels.Get(helind).GetFace (k); - const Point3d & p1 = mesh.Point (face.PNum(1)); - const Point3d & p2 = mesh.Point (face[1]); - const Point3d & p3 = mesh.Point (face[2]); - */ + const Point<3> & p1 = mesh.Point ( PointIndex (i3.I1()) ); + const Point<3> & p2 = mesh.Point ( PointIndex (i3.I2()) ); + const Point<3> & p3 = mesh.Point ( PointIndex (i3.I3()) ); - INDEX_3 i3 = tempels.Get(helind).GetFace (k-1); - - const Point3d & p1 = mesh.Point ( PointIndex (i3.I1()) ); - const Point3d & p2 = mesh.Point ( PointIndex (i3.I2()) ); - const Point3d & p3 = mesh.Point ( PointIndex (i3.I3()) ); - - - Vec3d v1(p1, p2); - Vec3d v2(p1, p3); - Vec3d n = Cross (v1, v2); + Vec<3> v1 = p2-p1; + Vec<3> v2 = p3-p1; + Vec<3> n = Cross (v1, v2); n /= n.Length(); - if (n * Vec3d (p1, mesh.Point (tempels.Get(helind)[k-1])) > 0) + if (n * Vec3d (p1, mesh.Point (tempels.Get(helind)[k])) > 0) n *= -1; double dist = n * Vec3d (p1, newp); - if (dist > -1e-10) // 1e-10 { insphere.Add (nbind); changed = 1; } - - } } } } // while (changed) - // (*testout) << "newels: " << endl; - Array newels; + + // Array newels; + Array newels; Element2d face(TRIG); - for (int j = 1; j <= insphere.GetArray().Size(); j++) - for (int k = 1; k <= 4; k++) + for (int celind : insphere.GetArray()) + for (int k = 0; k < 4; k++) { - // int elind = insphere.GetArray().Get(j); - int celind = insphere.GetArray().Get(j); int nbind = meshnb.GetNB (celind, k); if (!nbind || !insphere.IsIn (nbind)) { - tempels.Get (celind).GetFace (k-1, face); + tempels.Get (celind).GetFace (k, face); - Element newel(TET); + // Element newel(TET); + DelaunayTet newel; for (int l = 0; l < 3; l++) newel[l] = face[l]; newel[3] = newpi; @@ -483,7 +391,7 @@ namespace netgen n.Normalize(); if (n * Vec3d(mesh.Point (face[0]), - mesh.Point (tempels.Get(insphere.GetArray().Get(j))[k-1])) + mesh.Point (tempels.Get(celind)[k])) > 0) n *= -1; @@ -504,38 +412,36 @@ namespace netgen meshnb.ResetFaceHT (10*insphere.GetArray().Size()+1); - for (int j = 1; j <= insphere.GetArray().Size(); j++) + for (auto celind : insphere.GetArray()) { - // int elind = - int celind = insphere.GetArray().Get(j); - meshnb.Delete (celind); list.DeleteElement (celind); for (int k = 0; k < 4; k++) tempels.Elem(celind)[k] = -1; - // ((ADTree6&)tettree.Tree()).DeleteElement (celind); tettree.DeleteElement (celind); freelist.Append (celind); } - int hasclose = 0; - for (int j = 1; j <= closesphere.GetArray().Size(); j++) + + + bool hasclose = false; + for (int ind : closesphere.GetArray()) { - int ind = closesphere.GetArray().Get(j); if (!insphere.IsIn(ind) && fabs (Dist2 (centers.Get (ind), newp) - radi2.Get(ind)) < 1e-8 ) - hasclose = 1; + hasclose = true; } for (int j = 1; j <= newels.Size(); j++) { + const auto & newel = newels.Get(j); int nelind; if (!freelist.Size()) { - tempels.Append (newels.Get(j)); + tempels.Append (newel); nelind = tempels.Size(); } else @@ -543,29 +449,29 @@ namespace netgen nelind = freelist.Last(); freelist.DeleteLast(); - tempels.Elem(nelind) = newels.Get(j); + tempels.Elem(nelind) = newel; } meshnb.Add (nelind); list.AddElement (nelind); for (int k = 0; k < 4; k++) - pp[k] = &mesh.Point (newels.Get(j)[k]); + pp[k] = &mesh.Point (newel[k]); if (CalcSphereCenter (&pp[0], pc) ) { PrintSysError ("Delaunay: New tet is flat"); (*testout) << "new tet is flat" << endl; - for (int k = 1; k <= 4; k++) - (*testout) << newels.Get(j).PNum(k) << " "; + for (int k = 0; k < 4; k++) + (*testout) << newel[k] << " "; (*testout) << endl; - for (int k = 1; k <= 4; k++) + for (int k = 0; k < 4; k++) (*testout) << *pp[k-1] << " "; (*testout) << endl; } - r2 = Dist2 (*pp[0], pc); + double r2 = Dist2 (*pp[0], pc); if (hasclose) for (int k = 1; k <= closesphere.GetArray().Size(); k++) { @@ -615,33 +521,11 @@ namespace netgen int oldnp, DelaunayTet & startel, Point3d & pmin, Point3d & pmax) { static Timer t("Meshing3::Delaunay1"); RegionTimer reg(t); + static Timer tloop("Meshing3::Delaunay1 loop"); Array> centers; Array radi2; - Point3d tpmin, tpmax; - - - // new: local box - /* - mesh.GetBox (pmax, pmin); // lower bound for pmax, upper for pmin - for (int i = 1; i <= adfront->GetNF(); i++) - { - const MiniElement2d & face = adfront->GetFace(i); - for (PointIndex pi : face.PNums()) - { - pmin.SetToMin (mesh.Point (pi)); - pmax.SetToMax (mesh.Point (pi)); - } - } - - for (PointIndex pi : mesh.LockedPoints()) - { - pmin.SetToMin (mesh.Point (pi)); - pmax.SetToMax (mesh.Point (pi)); - } - cout << "pmin = " << pmin << ", pmax = " << pmax << endl; - */ Box<3> bbox(Box<3>::EMPTY_BOX); for (auto & face : adfront->Faces()) @@ -653,26 +537,24 @@ namespace netgen pmin = bbox.PMin(); pmax = bbox.PMax(); - // cout << "bbox = " << bbox << endl; - Vec3d vdiag(pmin, pmax); + + Vec<3> vdiag = pmax-pmin; // double r1 = vdiag.Length(); - double r1 = sqrt (3.0) * max3(vdiag.X(), vdiag.Y(), vdiag.Z()); - vdiag = Vec3d (r1, r1, r1); + double r1 = sqrt (3.0) * max3(vdiag(0), vdiag(1), vdiag(2)); + vdiag = Vec<3> (r1, r1, r1); //double r2; - Point3d pmin2 = pmin - 8 * vdiag; - Point3d pmax2 = pmax + 8 * vdiag; + Point<3> pmin2 = pmin - 8 * vdiag; + Point<3> pmax2 = pmax + 8 * vdiag; - Point3d cp1(pmin2), cp2(pmax2), cp3(pmax2), cp4(pmax2); - cp2.X() = pmin2.X(); - cp3.Y() = pmin2.Y(); - cp4.Z() = pmin2.Z(); + Point<3> cp1(pmin2), cp2(pmax2), cp3(pmax2), cp4(pmax2); + cp2(0) = pmin2(0); + cp3(1) = pmin2(1); + cp4(2) = pmin2(2); - - - int np = mesh.GetNP(); + size_t np = mesh.GetNP(); startel[0] = mesh.AddPoint (cp1); startel[1] = mesh.AddPoint (cp2); @@ -682,24 +564,21 @@ namespace netgen // flag points to use for Delaunay: BitArrayChar usep(np); usep.Clear(); - for (int i = 1; i <= adfront->GetNF(); i++) - { - const MiniElement2d & face = adfront->GetFace(i); - for (int j = 0; j < face.GetNP(); j++) - usep.Set (face[j]); - } - for (int i = oldnp + PointIndex::BASE; + for (auto & face : adfront->Faces()) + for (PointIndex pi : face.Face().PNums()) + usep.Set (pi); + + for (size_t i = oldnp + PointIndex::BASE; i < np + PointIndex::BASE; i++) usep.Set (i); - for (int i = 0; i < mesh.LockedPoints().Size(); i++) - usep.Set (mesh.LockedPoints()[i]); - + for (PointIndex pi : mesh.LockedPoints()) + usep.Set (pi); + Array freelist; - int cntp = 0; MeshNB meshnb (tempels, mesh.GetNP() + 5); @@ -716,13 +595,12 @@ namespace netgen list.AddElement (1); Array connected, treesearch; - - tpmin = tpmax = mesh.Point(startel[0]); - for (int k = 1; k < 4; k++) - { - tpmin.SetToMin (mesh.Point (startel[k])); - tpmax.SetToMax (mesh.Point (startel[k])); - } + Box<3> tbox(Box<3>::EMPTY_BOX); + for (size_t k = 0; k < 4; k++) + tbox.Add (mesh.Point(startel[k])); + Point<3> tpmin = tbox.PMin(); + Point<3> tpmax = tbox.PMax(); + tpmax = tpmax + 0.01 * (tpmax - tpmin); tettree.Insert (tpmin, tpmax, 1); @@ -754,6 +632,7 @@ namespace netgen for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End()-4; pi++) mixed[pi] = PointIndex ( (prim * pi) % np + PointIndex::BASE ); + tloop.Start(); for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End()-4; pi++) { if (pi % 1000 == 0) @@ -782,7 +661,8 @@ namespace netgen connected, treesearch, freelist, list, insphere, closesphere); } - + tloop.Stop(); + for (int i = tempels.Size(); i >= 1; i--) if (tempels.Get(i)[0] <= 0) tempels.DeleteElement (i); @@ -839,26 +719,27 @@ namespace netgen // improve delaunay - mesh by swapping !!!! Mesh tempmesh; - for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++) - tempmesh.AddPoint (mesh[pi]); + + for (auto & meshpoint : mesh.Points()) + tempmesh.AddPoint (meshpoint); - for (int i = 1; i <= tempels.Size(); i++) + for (auto & tempel : tempels) { Element el(4); for (int j = 0; j < 4; j++) - el[j] = tempels.Elem(i)[j]; + el[j] = tempel[j]; el.SetIndex (1); - const Point3d & lp1 = mesh.Point (el[0]); - const Point3d & lp2 = mesh.Point (el[1]); - const Point3d & lp3 = mesh.Point (el[2]); - const Point3d & lp4 = mesh.Point (el[3]); - Vec3d v1(lp1, lp2); - Vec3d v2(lp1, lp3); - Vec3d v3(lp1, lp4); + const Point<3> & lp1 = mesh.Point (el[0]); + const Point<3> & lp2 = mesh.Point (el[1]); + const Point<3> & lp3 = mesh.Point (el[2]); + const Point<3> & lp4 = mesh.Point (el[3]); + Vec<3> v1 = lp2-lp1; + Vec<3> v2 = lp3-lp1; + Vec<3> v3 = lp4-lp1; - Vec3d n = Cross (v1, v2); + Vec<3> n = Cross (v1, v2); double vol = n * v3; if (vol > 0) swap (el[2], el[3]); @@ -1261,6 +1142,7 @@ namespace netgen INDEX_3_HASHTABLE boundaryfaces(mesh.GetNOpenElements()/3+1); + /* for (int i = 1; i <= mesh.GetNOpenElements(); i++) { const Element2d & tri = mesh.OpenElement(i); @@ -1268,6 +1150,13 @@ namespace netgen i3.Sort(); boundaryfaces.PrepareSet (i3); } + */ + for (const Element2d & tri : mesh.OpenElements()) + { + INDEX_3 i3 (tri[0], tri[1], tri[2]); + i3.Sort(); + boundaryfaces.PrepareSet (i3); + } boundaryfaces.AllocateElements(); for (int i = 1; i <= mesh.GetNOpenElements(); i++) { @@ -1277,14 +1166,23 @@ namespace netgen boundaryfaces.Set (i3, 1); } + /* for (int i = 0; i < tempels.Size(); i++) for (int j = 0; j < 4; j++) tempels[i].NB(j) = 0; - + */ + for (auto & el : tempels) + for (int j = 0; j < 4; j++) + el.NB(j) = 0; + TABLE elsonpoint(mesh.GetNP()); + /* for (int i = 0; i < tempels.Size(); i++) { const DelaunayTet & el = tempels[i]; + */ + for (const DelaunayTet & el : tempels) + { INDEX_4 i4(el[0], el[1], el[2], el[3]); i4.Sort(); elsonpoint.IncSizePrepare (i4.I1()); @@ -1308,7 +1206,8 @@ namespace netgen INDEX_3_CLOSED_HASHTABLE faceht(100); Element2d hel(TRIG); - for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++) + // for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++) + for (PointIndex pi : mesh.Points().Range()) { faceht.SetSize (4 * elsonpoint[pi].Size()); for (int ii = 0; ii < elsonpoint[pi].Size(); ii++) diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index e3fe5025..d1f09b91 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -454,7 +454,8 @@ namespace netgen const Element2d & OpenElement(int i) const { return openelements.Get(i); } - + auto & OpenElements() const { return openelements; } + /// are also quads open elements bool HasOpenQuads () const;