diff --git a/CMakeLists.txt b/CMakeLists.txt index d2766638..e6df0df2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -172,7 +172,6 @@ if (USE_GUI) get_filename_component(MY_LIB_DIR ${TK_LIBRARY} DIRECTORY) install( DIRECTORY "${MY_LIB_DIR}/tcl8.6" DESTINATION lib COMPONENT netgen ) install( DIRECTORY "${MY_LIB_DIR}/tk8.6" DESTINATION lib COMPONENT netgen ) - install( DIRECTORY "${MY_LIB_DIR}/tix8.4.3" DESTINATION lib COMPONENT netgen ) install( DIRECTORY "${MY_LIB_DIR}/../bin" DESTINATION . COMPONENT netgen ) add_definitions(-DTOGL_WGL) else(WIN32) diff --git a/libsrc/csg/csgeom.cpp b/libsrc/csg/csgeom.cpp index a1abe6d6..0cb2c8c2 100644 --- a/libsrc/csg/csgeom.cpp +++ b/libsrc/csg/csgeom.cpp @@ -648,7 +648,8 @@ namespace netgen int CSGeometry :: SetTopLevelObject (Solid * sol, Surface * surf) { - return toplevelobjects.Append (new TopLevelObject (sol, surf)) - 1; + toplevelobjects.Append (new TopLevelObject (sol, surf)); + return toplevelobjects.Size()-1; } TopLevelObject * CSGeometry :: diff --git a/libsrc/csg/edgeflw.cpp b/libsrc/csg/edgeflw.cpp index 6afa4957..68cac826 100644 --- a/libsrc/csg/edgeflw.cpp +++ b/libsrc/csg/edgeflw.cpp @@ -617,7 +617,9 @@ namespace netgen (*testout) << "Point on edge" << endl << "seg = " << i2 << ", p = " << pi << endl << "pos = " << p << ", projected = " << hp << endl - << "seg is = " << mesh.Point(i2.I1()) << " - " << mesh.Point(i2.I2()) << endl; + << "seg is = " + << mesh.Point(PointIndex(i2.I1())) << " - " + << mesh.Point(PointIndex(i2.I2())) << endl; } } } @@ -654,9 +656,10 @@ namespace netgen } - for (int i = 1; i <= nseg; i++) + // for (int i = 1; i <= nseg; i++) + for (Segment & seg : mesh.LineSegments()) { - Segment & seg = mesh.LineSegment (i); + // Segment & seg = mesh.LineSegment (i); if (seg.edgenr >= 1 && seg.edgenr <= cntedge) { if (osedges.Get(seg.edgenr) != -1) @@ -1156,7 +1159,8 @@ namespace netgen //seg.surfnr2 = s2_rep; seg.surfnr1 = s1; seg.surfnr2 = s2; - hi = refedges.Append (seg); + refedges.Append (seg); + hi = refedges.Size(); refedgesinv.Append (edgeinv); edges_priority.Append((pre_ok[k-1]) ? 1 : 0); } diff --git a/libsrc/csg/genmesh.cpp b/libsrc/csg/genmesh.cpp index 566f94a1..fb4a8f14 100644 --- a/libsrc/csg/genmesh.cpp +++ b/libsrc/csg/genmesh.cpp @@ -759,8 +759,9 @@ namespace netgen << "time = " << GetTime() << " sec" << endl << "points: " << mesh->GetNP() << endl; #endif - - if (mparam.uselocalh && 0) + + /* + if (mparam.uselocalh) { mesh->CalcLocalH(mparam.grading); mesh->DeleteMesh(); @@ -773,6 +774,7 @@ namespace netgen MeshSurface (geom, *mesh, mparam); if (multithread.terminate) return TCL_OK; } + */ #ifdef LOG_STREAM (*logout) << "Surfaces remeshed" << endl diff --git a/libsrc/csg/polyhedra.cpp b/libsrc/csg/polyhedra.cpp index a09d7433..1bb1d1d2 100644 --- a/libsrc/csg/polyhedra.cpp +++ b/libsrc/csg/polyhedra.cpp @@ -499,7 +499,8 @@ int Polyhedra :: AddPoint (const Point<3> & p) else poly_bbox.Add(p); - return points.Append (p); + points.Append (p); + return points.Size(); } int Polyhedra :: AddFace (int pi1, int pi2, int pi3, int inputnum) diff --git a/libsrc/csg/specpoin.cpp b/libsrc/csg/specpoin.cpp index 10da4bb9..15d41d39 100644 --- a/libsrc/csg/specpoin.cpp +++ b/libsrc/csg/specpoin.cpp @@ -2027,7 +2027,8 @@ namespace netgen if (spi == -1) { - spi = specpoints.Append (SpecialPoint()) - 1; + specpoints.Append (SpecialPoint()); + spi = specpoints.Size()-1; specpoint2point.Append (i); specpoints.Last().unconditional = 0; searchtree.Insert (apoints[i], spi); diff --git a/libsrc/csg/vscsg.cpp b/libsrc/csg/vscsg.cpp index b35b4b10..8fab26bd 100644 --- a/libsrc/csg/vscsg.cpp +++ b/libsrc/csg/vscsg.cpp @@ -130,18 +130,18 @@ namespace netgen int hasp = 0; for (int i = 0; i < geometry->GetNTopLevelObjects(); i++) { - const TriangleApproximation & ta = - *geometry->GetTriApprox(i); - if (!&ta) continue; + const TriangleApproximation * ta = + geometry->GetTriApprox(i); + if (!ta) continue; - for (int j = 0; j < ta.GetNP(); j++) + for (int j = 0; j < ta->GetNP(); j++) { if (hasp) - box.Add (ta.GetPoint(j)); + box.Add (ta->GetPoint(j)); else { hasp = 1; - box.Set (ta.GetPoint(j)); + box.Set (ta->GetPoint(j)); } } } @@ -167,18 +167,18 @@ namespace netgen trilists.Append (glGenLists (1)); glNewList (trilists.Last(), GL_COMPILE); glEnable (GL_NORMALIZE); - const TriangleApproximation & ta = - *geometry->GetTriApprox(i); - if (&ta) + const TriangleApproximation * ta = + geometry->GetTriApprox(i); + if (ta) { glEnableClientState(GL_VERTEX_ARRAY); - glVertexPointer(3, GL_DOUBLE, 0, &ta.GetPoint(0)(0)); + glVertexPointer(3, GL_DOUBLE, 0, &ta->GetPoint(0)(0)); glEnableClientState(GL_NORMAL_ARRAY); - glNormalPointer(GL_DOUBLE, 0, &ta.GetNormal(0)(0)); + glNormalPointer(GL_DOUBLE, 0, &ta->GetNormal(0)(0)); - for (int j = 0; j < ta.GetNT(); j++) - glDrawElements(GL_TRIANGLES, 3, GL_UNSIGNED_INT, & (ta.GetTriangle(j)[0])); + for (int j = 0; j < ta->GetNT(); j++) + glDrawElements(GL_TRIANGLES, 3, GL_UNSIGNED_INT, & (ta->GetTriangle(j)[0])); glDisableClientState(GL_VERTEX_ARRAY); glDisableClientState(GL_NORMAL_ARRAY); @@ -382,13 +382,20 @@ namespace netgen glDisable (GL_COLOR_MATERIAL); glDisable (GL_LIGHTING); glDisable (GL_CLIP_PLANE0); - + + /* for (int i = 1; i <= mesh -> GetNP(); i++) { const Point3d & p = mesh -> Point(i); glRasterPos3d (p.X(), p.Y(), p.Z()); glBitmap (7, 7, 3, 3, 0, 0, &knoedel[0]); } + */ + for (const Point3d & p : mesh->Points()) + { + glRasterPos3d (p.X(), p.Y(), p.Z()); + glBitmap (7, 7, 3, 3, 0, 0, &knoedel[0]); + } } if (vispar.drawedpointnrs) @@ -403,12 +410,13 @@ namespace netgen // glListBase (fontbase); char buf[20]; - for (int i = 1; i <= mesh->GetNP(); i++) + // for (int i = 1; i <= mesh->GetNP(); i++) + for (auto i : mesh->Points().Range()) { const Point3d & p = mesh->Point(i); glRasterPos3d (p.X(), p.Y(), p.Z()); - sprintf (buf, "%d", i); + sprintf (buf, "%d", int(i)); // glCallLists (GLsizei(strlen (buf)), GL_UNSIGNED_BYTE, buf); MyOpenGLText (buf); } diff --git a/libsrc/general/array.hpp b/libsrc/general/array.hpp index 69cf2f30..9c6bcbe5 100644 --- a/libsrc/general/array.hpp +++ b/libsrc/general/array.hpp @@ -293,13 +293,13 @@ namespace netgen /// Add element at end of array. reallocation if necessary. - int Append (const T & el) + void Append (const T & el) { if (size == allocsize) ReSize (size+1); data[size] = el; size++; - return size; + // return size; } template diff --git a/libsrc/include/mydefs.hpp b/libsrc/include/mydefs.hpp index f9118daf..bf3cb1e9 100644 --- a/libsrc/include/mydefs.hpp +++ b/libsrc/include/mydefs.hpp @@ -38,7 +38,8 @@ - +// #define BASE0 +// #define DEBUG #define noDEMOVERSION diff --git a/libsrc/include/nginterface_v2.hpp b/libsrc/include/nginterface_v2.hpp index 742e4258..5371a97e 100644 --- a/libsrc/include/nginterface_v2.hpp +++ b/libsrc/include/nginterface_v2.hpp @@ -14,15 +14,20 @@ namespace netgen { + + static constexpr int POINTINDEX_BASE = 1; + struct T_EDGE2 { - int orient:1; - int nr:31; // 0-based + // int orient:1; + // int nr:31; // 0-based + int nr; // 0-based }; struct T_FACE2 { - int orient:3; - int nr:29; // 0-based + // int orient:3; + // int nr:29; // 0-based + int nr; // 0-based }; class Ng_Element @@ -35,7 +40,7 @@ namespace netgen const int * ptr; int Size() const { return num; } - int operator[] (int i) const { return ptr[i]-1; } + int operator[] (int i) const { return ptr[i]-POINTINDEX_BASE; } }; @@ -46,7 +51,7 @@ namespace netgen const int * ptr; int Size() const { return num; } - int operator[] (int i) const { return ptr[i]-1; } + int operator[] (int i) const { return ptr[i]-POINTINDEX_BASE; } }; class Ng_Edges @@ -71,6 +76,17 @@ namespace netgen int operator[] (int i) const { return ptr[i].nr; } }; + class Ng_Facets + { + public: + int num; + const int * ptr; + + int Size() const { return num; } + int operator[] (int i) const { return ptr[i]; } + }; + + public: NG_ELEMENT_TYPE type; int index; // material / boundary condition @@ -81,6 +97,7 @@ namespace netgen Ng_Vertices vertices; Ng_Edges edges; Ng_Faces faces; + Ng_Facets facets; bool is_curved; }; @@ -131,7 +148,7 @@ namespace netgen const int * ptr; int Size() const { return 2; } - int operator[] (int i) const { return ptr[i]-1; } + int operator[] (int i) const { return ptr[i]-POINTINDEX_BASE; } }; @@ -151,7 +168,7 @@ namespace netgen const int * ptr; int Size() const { return nv; } - int operator[] (int i) const { return ptr[i]-1; } + int operator[] (int i) const { return ptr[i]-POINTINDEX_BASE; } }; class Ng_Edges diff --git a/libsrc/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp index dfdaf46b..78b829e2 100644 --- a/libsrc/include/nginterface_v2_impl.hpp +++ b/libsrc/include/nginterface_v2_impl.hpp @@ -1,6 +1,6 @@ NGX_INLINE DLL_HEADER Ng_Point Ngx_Mesh :: GetPoint (int nr) const { - return Ng_Point (&mesh->Point(nr + PointIndex::BASE)(0)); + return Ng_Point (&mesh->Point(PointIndex(nr+PointIndex::BASE))(0)); } @@ -86,6 +86,20 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (int nr) const ret.faces.num = 0; ret.faces.ptr = NULL; + if (mesh->GetDimension() == 2) + { + ret.facets.num = 1; + ret.facets.ptr = (int*)ret.edges.ptr; + } + else + { + ret.facets.num = 0; + ret.facets.ptr = nullptr; + // not working as long as vertices are 1-based + // ret.facets.num = 2; + // ret.facets.ptr = (int*)&(el[0]); + } + // ret.is_curved = mesh->GetCurvedElements().IsSegmentCurved(nr); ret.is_curved = el.IsCurved(); @@ -117,6 +131,16 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<2> (int nr) const ret.faces.num = MeshTopology::GetNFaces (el.GetType()); ret.faces.ptr = (T_FACE2*)mesh->GetTopology().GetSurfaceElementFacesPtr (nr); + if (mesh->GetDimension() == 3) + { + ret.facets.num = ret.faces.num; + ret.facets.ptr = (int*)ret.faces.ptr; + } + else + { + ret.facets.num = ret.edges.num; + ret.facets.ptr = (int*)ret.edges.ptr; + } ret.is_curved = el.IsCurved(); return ret; } @@ -142,6 +166,9 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<3> (int nr) const ret.faces.num = MeshTopology::GetNFaces (el.GetType()); ret.faces.ptr = (T_FACE2*)mesh->GetTopology().GetElementFacesPtr (nr); + ret.facets.num = ret.faces.num; + ret.facets.ptr = (int*)ret.faces.ptr; + ret.is_curved = el.IsCurved(); return ret; } diff --git a/libsrc/interface/writejcm.cpp b/libsrc/interface/writejcm.cpp index 0b2ced88..59b64101 100644 --- a/libsrc/interface/writejcm.cpp +++ b/libsrc/interface/writejcm.cpp @@ -33,7 +33,7 @@ void WriteJCMFormat (const Mesh & mesh, int np = mesh.GetNP(); // Identic points - Array identmap1, identmap2, identmap3; + Array identmap1, identmap2, identmap3; mesh.GetIdentifications().GetMap(1, identmap1); mesh.GetIdentifications().GetMap(2, identmap2); mesh.GetIdentifications().GetMap(3, identmap3); diff --git a/libsrc/interface/wuchemnitz.cpp b/libsrc/interface/wuchemnitz.cpp index 8641a8cd..46872fb0 100644 --- a/libsrc/interface/wuchemnitz.cpp +++ b/libsrc/interface/wuchemnitz.cpp @@ -202,7 +202,8 @@ namespace netgen fa.p1 = i3.I1(); fa.p2 = i3.I2(); fa.p3 = i3.I3(); - facei = faces.Append (fa); + faces.Append (fa); + facei = faces.Size(); faceindex.Set (i3, facei); } @@ -241,7 +242,8 @@ namespace netgen EDGE ed; ed.p1 = i2.I1(); ed.p2 = i2.I2(); - edgei = edges.Append (ed); + edges.Append (ed); + edgei = edges.Size(); edgeindex.Set (i2, edgei); } diff --git a/libsrc/meshing/adfront2.cpp b/libsrc/meshing/adfront2.cpp index 647de0b4..7eb426ad 100644 --- a/libsrc/meshing/adfront2.cpp +++ b/libsrc/meshing/adfront2.cpp @@ -88,7 +88,8 @@ namespace netgen } else { - pi = points.Append (FrontPoint2 (p, globind, mgi, pointonsurface)) - 1; + points.Append (FrontPoint2 (p, globind, mgi, pointonsurface)); + pi = points.Size()-1; } if (mgi) @@ -128,7 +129,8 @@ namespace netgen } else { - li = lines.Append(FrontLine (INDEX_2(pi1, pi2))) - 1; + lines.Append(FrontLine (INDEX_2(pi1, pi2))); + li = lines.Size()-1; } @@ -332,7 +334,8 @@ namespace netgen { pindex.Append (pi); invpindex[pi] = pindex.Size(); - loclines[i][j] = locpoints.Append (points[pi].P()); + locpoints.Append (points[pi].P()); + loclines[i][j] = locpoints.Size(); } else loclines[i][j] = invpindex[pi]; @@ -349,7 +352,8 @@ namespace netgen // Dist2 (points.Get(i).P(), p0) <= xh2 && invpindex[i] <= 0) { - invpindex[i] = locpoints.Append (points[i].P()); + locpoints.Append (points[i].P()); + invpindex[i] = locpoints.Size(); pindex.Append(i); } } diff --git a/libsrc/meshing/adfront2.hpp b/libsrc/meshing/adfront2.hpp index 4d76faaf..1bb21c40 100644 --- a/libsrc/meshing/adfront2.hpp +++ b/libsrc/meshing/adfront2.hpp @@ -35,7 +35,7 @@ /// FrontPoint2 () { - globalindex = -1; + globalindex.Invalidate(); // = -1; nlinetopoint = 0; frontnr = INT_MAX-10; // attention: overflow on calculating INT_MAX + 1 mgi = NULL; diff --git a/libsrc/meshing/adfront3.cpp b/libsrc/meshing/adfront3.cpp index 14d721d7..089b96ab 100644 --- a/libsrc/meshing/adfront3.cpp +++ b/libsrc/meshing/adfront3.cpp @@ -9,7 +9,7 @@ namespace netgen FrontPoint3 :: FrontPoint3 () { - globalindex = -1; + globalindex.Invalidate(); // = -1; nfacetopoint = 0; frontnr = 1000; cluster = 0; @@ -159,8 +159,9 @@ INDEX AdFront3 :: AddFace (const MiniElement2d & aface) for (i = 1; i <= aface.GetNP(); i++) points[aface.PNum(i)].DecFrontNr (minfn+1); - - int nfn = faces.Append(FrontFace (aface)); + + faces.Append(FrontFace (aface)); + int nfn = faces.Size(); faces.Elem(nfn).cluster = cluster; if (hashon && hashcreated) @@ -484,9 +485,9 @@ int AdFront3 :: SelectBaseElement () int AdFront3 :: GetLocals (int fstind, - Array & locpoints, + Array & locpoints, Array & locfaces, // local index - Array & pindex, + Array & pindex, Array & findex, INDEX_2_HASHTABLE & getconnectedpairs, float xh, @@ -600,12 +601,13 @@ int AdFront3 :: GetLocals (int fstind, if (invpindex[pi] == -1) { pindex.Append (pi); - invpindex[pi] = pindex.Size(); // -1+PointIndex::BASE; - locfaces.Elem(i).PNum(j) = locpoints.Append (points[pi].P()); - } - else - locfaces.Elem(i).PNum(j) = invpindex[pi]; - + locpoints.Append (points[pi].P()); + invpindex[pi] = pindex.Size()-1+PointIndex::BASE; + } + // locfaces.Elem(i).PNum(j) = locpoints.Append (points[pi].P()); + // } + // else + locfaces.Elem(i).PNum(j) = invpindex[pi]; } } @@ -655,9 +657,9 @@ int AdFront3 :: GetLocals (int fstind, // returns all points connected with fi void AdFront3 :: GetGroup (int fi, - Array & grouppoints, + Array & grouppoints, Array & groupelements, - Array & pindex, + Array & pindex, Array & findex) { // static Array pingroup; diff --git a/libsrc/meshing/adfront3.hpp b/libsrc/meshing/adfront3.hpp index 11993733..9243db42 100644 --- a/libsrc/meshing/adfront3.hpp +++ b/libsrc/meshing/adfront3.hpp @@ -17,15 +17,15 @@ class FrontPoint3 { /// coordinates -Point<3> p; + Point<3> p; /// global node index -PointIndex globalindex; + PointIndex globalindex; /// number of faces connected to point -int nfacetopoint; + int nfacetopoint; /// distance to original boundary -int frontnr; + int frontnr; /// -int cluster; + int cluster; public: /// FrontPoint3 (); @@ -43,7 +43,7 @@ public: void AddFace () { nfacetopoint++; } - /// + /// if last face is removed, then point is invalidated void RemoveFace() { nfacetopoint--; @@ -51,7 +51,7 @@ public: } /// - int Valid () const + bool Valid () const { return nfacetopoint >= 0; } /// @@ -74,7 +74,7 @@ class MiniElement2d { protected: int np; - PointIndex pnum[4]; + PointIndex pnum[4]; // can be global or local nums bool deleted; public: MiniElement2d () @@ -89,8 +89,8 @@ public: const PointIndex PNum (int i) const { return pnum[i-1]; } PointIndex & PNum (int i) { return pnum[i-1]; } const PointIndex PNumMod (int i) const { return pnum[(i-1)%np]; } - - void Delete () { deleted = 1; pnum[0] = pnum[1] = pnum[2] = pnum[3] = PointIndex::BASE-1; } + auto PNums() const { return FlatArray (np, &pnum[0]); } + void Delete () { deleted = true; for (PointIndex & p : pnum) p.Invalidate(); } bool IsDeleted () const { return deleted; } }; @@ -120,7 +120,7 @@ private: int hashvalue; /// int cluster; - + public: /// FrontFace (); @@ -178,43 +178,43 @@ class AdFront3 /// Array points; /// -Array faces; + Array faces; /// -Array delpointl; - + Array delpointl; + /// which points are connected to pi ? -TABLE * connectedpairs; + TABLE * connectedpairs; /// number of total front faces; -int nff; + int nff; /// number of quads in front -int nff4; + int nff4; /// -double vol; - + double vol; + /// -GeomSearch3d hashtable; - + GeomSearch3d hashtable; + /// -int hashon; + int hashon; /// -int hashcreated; - + int hashcreated; + /// counter for rebuilding internal tables -int rebuildcounter; + int rebuildcounter; /// last base element -int lasti; + int lasti; /// minimal selection-value of baseelements -int minval; - Array invpindex; + int minval; + Array invpindex; Array pingroup; /// -class Box3dTree * facetree; + class Box3dTree * facetree; public: - + /// AdFront3 (); /// @@ -260,9 +260,9 @@ public: /// int GetLocals (int baseelement, - Array & locpoints, + Array & locpoints, Array & locfaces, // local index - Array & pindex, + Array & pindex, Array & findex, INDEX_2_HASHTABLE & connectedpairs, float xh, @@ -271,9 +271,9 @@ public: /// void GetGroup (int fi, - Array & grouppoints, + Array & grouppoints, Array & groupelements, - Array & pindex, + Array & pindex, Array & findex); /// diff --git a/libsrc/meshing/clusters.cpp b/libsrc/meshing/clusters.cpp index 7e656b17..c19862b0 100644 --- a/libsrc/meshing/clusters.cpp +++ b/libsrc/meshing/clusters.cpp @@ -96,7 +96,7 @@ namespace netgen nnums.SetSize(elnv+elned+elnfa+1); for (int j = 1; j <= elnv; j++) - nnums.Elem(j) = el.PNum(j); + nnums.Elem(j) = el.PNum(j)+1-PointIndex::BASE; for (int j = 1; j <= elned; j++) nnums.Elem(elnv+j) = nv+ednums.Elem(j); for (int j = 1; j <= elnfa; j++) @@ -108,7 +108,6 @@ namespace netgen } }); - NgProfiler::StopTimer(timer1); NgProfiler::StartTimer(timer2); @@ -125,7 +124,7 @@ namespace netgen nnums.SetSize(elnv+elned+1); for (int j = 1; j <= elnv; j++) - nnums.Elem(j) = el.PNum(j); + nnums.Elem(j) = el.PNum(j)+1-PointIndex::BASE; for (int j = 1; j <= elned; j++) nnums.Elem(elnv+j) = nv+ednums.Elem(j); nnums.Elem(elnv+elned+1) = fanum; @@ -197,7 +196,7 @@ namespace netgen nnums.SetSize(elnv+elned+elnfa+1); for (int j = 1; j <= elnv; j++) - nnums.Elem(j) = el.PNum(j); + nnums.Elem(j) = el.PNum(j)+1-PointIndex::BASE; for (int j = 1; j <= elned; j++) nnums.Elem(elnv+j) = nv+ednums.Elem(j); for (int j = 1; j <= elnfa; j++) @@ -220,23 +219,23 @@ namespace netgen break; case TET: case TET10: - if (cluster_reps.Get(el.PNum(1)) == - cluster_reps.Get(el.PNum(2))) + if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) == + cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE)) clustertab = tet_cluster12; - else if (cluster_reps.Get(el.PNum(1)) == - cluster_reps.Get(el.PNum(3))) + else if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) == + cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE)) clustertab = tet_cluster13; - else if (cluster_reps.Get(el.PNum(1)) == - cluster_reps.Get(el.PNum(4))) + else if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) == + cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE)) clustertab = tet_cluster14; - else if (cluster_reps.Get(el.PNum(2)) == - cluster_reps.Get(el.PNum(3))) + else if (cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE) == + cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE)) clustertab = tet_cluster23; - else if (cluster_reps.Get(el.PNum(2)) == - cluster_reps.Get(el.PNum(4))) + else if (cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE) == + cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE)) clustertab = tet_cluster24; - else if (cluster_reps.Get(el.PNum(3)) == - cluster_reps.Get(el.PNum(4))) + else if (cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE) == + cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE)) clustertab = tet_cluster34; else diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index 14b0d4df..48538359 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -3873,7 +3873,7 @@ namespace netgen Mat _dxdxi; if (!EvaluateMapping (info, _xi, _x, _dxdxi)) { ok = false; break; } - // cout << "x = " << _x << ", dxdxi = " << _dxdxi << endl; + // *testout << "x = " << _x << ", dxdxi = " << _dxdxi << endl; if (x) for (int j = 0; j < DIM_SPACE; j++) x[i*sx+j] = _x[j]; diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index 9cb37644..69bffb0f 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -415,9 +415,9 @@ namespace netgen INDEX_3 i3 = tempels.Get(helind).GetFace (k-1); - const Point3d & p1 = mesh.Point ( i3.I1()); - const Point3d & p2 = mesh.Point ( i3.I2()); - const Point3d & p3 = mesh.Point ( i3.I3()); + 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); @@ -612,20 +612,20 @@ namespace netgen for (int i = 1; i <= adfront->GetNF(); i++) { const MiniElement2d & face = adfront->GetFace(i); - for (int j = 0; j < face.GetNP(); j++) + for (PointIndex pi : face.PNums()) { - pmin.SetToMin (mesh.Point (face[j])); - pmax.SetToMax (mesh.Point (face[j])); + pmin.SetToMin (mesh.Point (pi)); + pmax.SetToMax (mesh.Point (pi)); } } - for (int i = 0; i < mesh.LockedPoints().Size(); i++) + for (PointIndex pi : mesh.LockedPoints()) { - pmin.SetToMin (mesh.Point (mesh.LockedPoints()[i])); - pmax.SetToMax (mesh.Point (mesh.LockedPoints()[i])); + pmin.SetToMin (mesh.Point (pi)); + pmax.SetToMax (mesh.Point (pi)); } - - + + Vec3d vdiag(pmin, pmax); // double r1 = vdiag.Length(); diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp index 43b67947..7aa5b43d 100644 --- a/libsrc/meshing/improve3.cpp +++ b/libsrc/meshing/improve3.cpp @@ -274,16 +274,8 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh, void MeshOptimize3d :: SplitImprove (Mesh & mesh, OPTIMIZEGOAL goal) { - int j, k, l; - Point3d p1, p2, pnew; - - ElementIndex ei; - SurfaceElementIndex sei; - PointIndex pi1, pi2; - double bad1, bad2, badmax, badlimit; - int cnt = 0; int np = mesh.GetNP(); int ne = mesh.GetNE(); @@ -291,7 +283,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, TABLE elementsonnode(np); Array hasbothpoints; - BitArray origpoint(np), boundp(np); + BitArray origpoint(np+1), boundp(np+1); // big enough for 0 and 1-based origpoint.Set(); Array elerrs(ne); @@ -301,7 +293,6 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, const char * savetask = multithread.task; multithread.task = "Split Improve"; - PrintMessage (3, "SplitImprove"); (*testout) << "start SplitImprove" << "\n"; @@ -311,10 +302,11 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, bad1 = 0; badmax = 0; - for (ei = 0; ei < ne; ei++) + for (ElementIndex ei = 0; ei < ne; ei++) { - if(mesh.GetDimension()==3 && mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex()) + if(mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex()) continue; + elerrs[ei] = CalcBad (mesh.Points(), mesh[ei], 0); bad1 += elerrs[ei]; if (elerrs[ei] > badmax) badmax = elerrs[ei]; @@ -325,41 +317,40 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, boundp.Clear(); - for (sei = 0; sei < mesh.GetNSE(); sei++) - for (j = 0; j < 3; j++) - boundp.Set (mesh[sei][j]); - + for (auto & el : mesh.SurfaceElements()) + for (PointIndex pi : el.PNums()) + boundp.Set (pi); + if (goal == OPT_QUALITY) { bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements()); (*testout) << "Total badness = " << bad1 << endl; } - for (ei = 0; ei < ne; ei++) - for (j = 0; j < mesh[ei].GetNP(); j++) - elementsonnode.Add (mesh[ei][j], ei); - + for (ElementIndex ei : mesh.VolumeElements().Range()) + for (PointIndex pi : mesh[ei].PNums()) + elementsonnode.Add (pi, ei); mesh.MarkIllegalElements(); if (goal == OPT_QUALITY || goal == OPT_LEGAL) { int cntill = 0; - for (ei = 0; ei < ne; ei++) + for (ElementIndex ei = 0; ei < ne; ei++) { - // if (!LegalTet (volelements.Get(i))) - if (mesh[ei].flags.illegal) - { - cntill++; - illegaltet.Set (ei+1); - } + // if (!LegalTet (volelements.Get(i))) + if (mesh[ei].flags.illegal) + { + cntill++; + illegaltet.Set (ei); + } } - // (*mycout) << cntill << " illegal tets" << endl; } - - for (ei = 0; ei < ne; ei++) + for (ElementIndex ei : mesh.VolumeElements().Range()) { - if(mesh.GetDimension()==3 && mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex()) + Element & elem = mesh[ei]; + + if(mp.only3D_domain_nr && mp.only3D_domain_nr != elem.GetIndex()) continue; if (multithread.terminate) break; @@ -368,82 +359,74 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, bool ltestmode = 0; - - if (elerrs[ei] < badlimit && !illegaltet.Test(ei+1)) continue; + if (elerrs[ei] < badlimit && !illegaltet.Test(ei)) continue; if ((goal == OPT_LEGAL) && - !illegaltet.Test(ei+1) && - CalcBad (mesh.Points(), mesh[ei], 0) < 1e3) + !illegaltet.Test(ei) && + CalcBad (mesh.Points(), elem, 0) < 1e3) continue; - - Element & elem = mesh[ei]; - if (ltestmode) { (*testout) << "test el " << ei << endl; - for (j = 0; j < 4; j++) + for (int j = 0; j < 4; j++) (*testout) << elem[j] << " "; (*testout) << endl; } - - for (j = 0; j < 6; j++) + for (int j = 0; j < 6; j++) { - static const int tetedges[6][2] = { { 0, 1 }, { 0, 2 }, { 0, 3 }, { 1, 2 }, { 1, 3 }, { 2, 3 } }; - pi1 = elem[tetedges[j][0]]; - pi2 = elem[tetedges[j][1]]; + PointIndex pi1 = elem[tetedges[j][0]]; + PointIndex pi2 = elem[tetedges[j][1]]; if (pi2 < pi1) Swap (pi1, pi2); - if (pi2 > elementsonnode.Size()) continue; + if (pi2 >= elementsonnode.Size()+PointIndex::BASE) continue; // old number of points if (!origpoint.Test(pi1) || !origpoint.Test(pi2)) continue; - INDEX_2 i2(pi1, pi2); i2.Sort(); if (mesh.BoundaryEdge (pi1, pi2)) continue; - if (edgetested.Used (i2) && !illegaltet.Test(ei+1)) continue; + if (edgetested.Used (i2) && !illegaltet.Test(ei)) continue; edgetested.Set (i2, 1); hasbothpoints.SetSize (0); - for (k = 1; k <= elementsonnode.EntrySize(pi1); k++) + /* + for (int k = 1; k <= elementsonnode.EntrySize(pi1); k++) { - bool has1 = 0, has2 = 0; - ElementIndex elnr = elementsonnode.Get(pi1, k); - Element & el = mesh[elnr]; + */ + for (ElementIndex ei : elementsonnode[pi1]) + { + Element & el = mesh[ei]; + bool has1 = el.PNums().Contains(pi1); + bool has2 = el.PNums().Contains(pi2); - for (l = 0; l < el.GetNP(); l++) - { - if (el[l] == pi1) has1 = 1; - if (el[l] == pi2) has2 = 1; - } if (has1 && has2) - if (!hasbothpoints.Contains (elnr)) - hasbothpoints.Append (elnr); + if (!hasbothpoints.Contains (ei)) + hasbothpoints.Append (ei); } bad1 = 0; - for (k = 0; k < hasbothpoints.Size(); k++) - bad1 += CalcBad (mesh.Points(), mesh[hasbothpoints[k]], 0); - + for (ElementIndex ei : hasbothpoints) + bad1 += CalcBad (mesh.Points(), mesh[ei], 0); + bool puretet = 1; - for (k = 0; k < hasbothpoints.Size(); k++) - if (mesh[hasbothpoints[k]].GetType() != TET) + for (ElementIndex ei : hasbothpoints) + if (mesh[ei].GetType() != TET) puretet = 0; if (!puretet) continue; - p1 = mesh[pi1]; - p2 = mesh[pi2]; + Point3d p1 = mesh[pi1]; + Point3d p2 = mesh[pi2]; /* pnew = Center (p1, p2); @@ -463,13 +446,12 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, points.Elem(pi2) = p2; */ - locfaces.SetSize (0); - for (k = 0; k < hasbothpoints.Size(); k++) + for (ElementIndex ei : hasbothpoints) { - const Element & el = mesh[hasbothpoints[k]]; - - for (l = 0; l < 4; l++) + const Element & el = mesh[ei]; + + for (int l = 0; l < 4; l++) if (el[l] == pi1 || el[l] == pi2) { INDEX_3 i3; @@ -486,7 +468,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, par.maxit_linsearch = 50; par.maxit_bfgs = 20; - pnew = Center (p1, p2); + Point3d pnew = Center (p1, p2); Vector px(3); px(0) = pnew.X(); px(1) = pnew.Y(); @@ -501,11 +483,10 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, pnew.Y() = px(1); pnew.Z() = px(2); - - int hpinew = mesh.AddPoint (pnew); + PointIndex hpinew = mesh.AddPoint (pnew); // ptyps.Append (INNERPOINT); - for (k = 0; k < hasbothpoints.Size(); k++) + for (int k = 0; k < hasbothpoints.Size(); k++) { Element & oldel = mesh[hasbothpoints[k]]; Element newel1 = oldel; @@ -515,7 +496,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, newel1.flags.illegal_valid = 0; newel2.flags.illegal_valid = 0; - for (l = 0; l < 4; l++) + for (int l = 0; l < 4; l++) { if (newel1[l] == pi2) newel1[l] = hpinew; if (newel2[l] == pi1) newel2[l] = hpinew; @@ -536,15 +517,15 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, PointIndex pinew = mesh.AddPoint (pnew); - for (k = 0; k < hasbothpoints.Size(); k++) + for (ElementIndex ei : hasbothpoints) { - Element & oldel = mesh[hasbothpoints[k]]; + Element & oldel = mesh[ei]; Element newel = oldel; - + newel.flags.illegal_valid = 0; oldel.flags.illegal_valid = 0; - for (l = 0; l < 4; l++) + for (int l = 0; l < 4; l++) { origpoint.Clear (oldel[l]); @@ -554,7 +535,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, mesh.AddVolumeElement (newel); } - j = 10; + j = 10; // end j-loop } } } @@ -572,11 +553,9 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, int cntill = 0; ne = mesh.GetNE(); - for (ei = 0; ei < ne; ei++) - { - if (!mesh.LegalTet (mesh[ei])) - cntill++; - } + for (ElementIndex ei = 0; ei < ne; ei++) + if (!mesh.LegalTet (mesh[ei])) + cntill++; // cout << cntill << " illegal tets" << endl; } @@ -590,7 +569,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, const BitArray * working_elements) { - PointIndex pi1(0), pi2(0), pi3(0), pi4(0), pi5(0), pi6(0); + PointIndex pi3(0), pi4(0), pi5(0), pi6(0); int cnt = 0; Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET); @@ -636,8 +615,12 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, // find elements on node for (ElementIndex ei = 0; ei < ne; ei++) + for (PointIndex pi : mesh[ei].PNums()) + elementsonnode.Add (pi, ei); + /* for (int j = 0; j < mesh[ei].GetNP(); j++) elementsonnode.Add (mesh[ei][j], ei); + */ // INDEX_2_HASHTABLE edgeused(2 * ne + 5); @@ -648,7 +631,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, if (multithread.terminate) break; - if(mesh.GetDimension()==3 && mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex()) + if (mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex()) continue; multithread.percent = 100.0 * (ei+1) / ne; @@ -678,18 +661,14 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, const Element & elemi = mesh[ei]; if (elemi.IsDeleted()) continue; - - // (*testout) << "check element " << elemi << endl; - int mattyp = elemi.GetIndex(); static const int tetedges[6][2] = { { 0, 1 }, { 0, 2 }, { 0, 3 }, { 1, 2 }, { 1, 3 }, { 2, 3 } }; - pi1 = elemi[tetedges[j][0]]; - pi2 = elemi[tetedges[j][1]]; - + PointIndex pi1 = elemi[tetedges[j][0]]; + PointIndex pi2 = elemi[tetedges[j][1]]; if (pi2 < pi1) Swap (pi1, pi2); @@ -718,21 +697,22 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, if (has1 && has2) { // only once - for (int l = 0; l < hasbothpoints.Size(); l++) - if (hasbothpoints[l] == elnr) - has1 = 0; - + if (hasbothpoints.Contains (elnr)) + has1 = false; + if (has1) - hasbothpoints.Append (elnr); + { + hasbothpoints.Append (elnr); + } } } - bool puretet = 1; - for (int k = 0; k < hasbothpoints.Size(); k++) - if (mesh[hasbothpoints[k]].GetType () != TET) - puretet = 0; - if (!puretet) - continue; + bool puretet = true; + for (ElementIndex ei : hasbothpoints) + if (mesh[ei].GetType () != TET) + puretet = false; + + if (!puretet) continue; int nsuround = hasbothpoints.Size(); @@ -763,10 +743,10 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, for (int k = 0; k < 3; k++) // JS, 201212 { const Element & elemk = mesh[hasbothpoints[k]]; - bool has1 = 0; + bool has1 = false; for (int l = 0; l < 4; l++) if (elemk[l] == pi4) - has1 = 1; + has1 = true; if (has1) { for (int l = 0; l < 4; l++) @@ -775,10 +755,9 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, } } - if(pi5 == 0) + if (!pi5.IsValid()) throw NgException("Illegal state observed in SwapImprove"); - el32[0] = pi1; el32[1] = pi2; @@ -882,7 +861,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, mesh[hasbothpoints[0]] = el21; mesh[hasbothpoints[1]] = el22; for (int l = 0; l < 4; l++) - mesh[hasbothpoints[2]][l] = 0; + mesh[hasbothpoints[2]][l].Invalidate(); mesh[hasbothpoints[2]].Delete(); for (int k = 0; k < 2; k++) @@ -912,14 +891,11 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, el1[3] = pi4; } - pi5 = 0; + pi5.Invalidate(); for (int k = 0; k < 4; k++) { const Element & elem = mesh[hasbothpoints[k]]; - bool has1 = 0; - for (int l = 0; l < 4; l++) - if (elem[l] == pi4) - has1 = 1; + bool has1 = elem.PNums().Contains(pi4); if (has1) { for (int l = 0; l < 4; l++) @@ -928,14 +904,11 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, } } - pi6 = 0; + pi6.Invalidate(); for (int k = 0; k < 4; k++) { const Element & elem = mesh[hasbothpoints[k]]; - bool has1 = 0; - for (int l = 0; l < 4; l++) - if (elem[l] == pi3) - has1 = 1; + bool has1 = elem.PNums().Contains(pi3); if (has1) { for (int l = 0; l < 4; l++) @@ -1202,7 +1175,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, Element hel(TET); ArrayMem suroundpts(nsuround); - ArrayMem tetused(nsuround); + ArrayMem tetused(nsuround); Element & elem = mesh[hasbothpoints[0]]; @@ -1228,33 +1201,32 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, // suroundpts.SetSize (nsuround); + suroundpts = -17; suroundpts[0] = pi3; suroundpts[1] = pi4; - tetused = 0; - tetused[0] = 1; + tetused = false; + tetused[0] = true; for (int l = 2; l < nsuround; l++) { - int oldpi = suroundpts[l-1]; - int newpi = 0; + PointIndex oldpi = suroundpts[l-1]; + PointIndex newpi; + newpi.Invalidate(); - - for (int k = 0; k < nsuround && !newpi; k++) + for (int k = 0; k < nsuround && !newpi.IsValid(); k++) if (!tetused[k]) { const Element & nel = mesh[hasbothpoints[k]]; - - for (int k2 = 0; k2 < 4 && !newpi; k2++) + for (int k2 = 0; k2 < 4 && !newpi.IsValid(); k2++) if (nel[k2] == oldpi) { newpi = nel[0] + nel[1] + nel[2] + nel[3] - pi1 - pi2 - oldpi; - tetused[k] = 1; + tetused[k] = true; suroundpts[l] = newpi; - } } } @@ -1406,8 +1378,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, */ rel.Delete(); for (int k1 = 0; k1 < 4; k1++) - rel[k1] = 0; - + rel[k1].Invalidate(); } } } diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 7475e512..d76248a4 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -5367,30 +5367,9 @@ namespace netgen int Mesh :: MarkIllegalElements () { int cnt = 0; - int i; - - for (i = 1; i <= GetNE(); i++) - { - LegalTet (VolumeElement(i)); - - /* - Element & el = VolumeElement(i); - int leg1 = LegalTet (el); - el.flags.illegal_valid = 0; - int leg2 = LegalTet (el); - - if (leg1 != leg2) - { - cerr << "legal differs!!" << endl; - (*testout) << "legal differs" << endl; - (*testout) << "elnr = " << i << ", el = " << el - << " leg1 = " << leg1 << ", leg2 = " << leg2 << endl; - } - - // el.flags.illegal = !LegalTet (el); - */ - cnt += VolumeElement(i).Illegal(); - } + for (auto & el : VolumeElements()) + if (!LegalTet (el)) + cnt++; return cnt; } diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index d09b7c98..83565105 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -24,8 +24,8 @@ namespace netgen { public: typedef ::netgen::T_POINTS T_POINTS; - typedef Array T_VOLELEMENTS; - typedef Array T_SURFELEMENTS; + typedef Array T_VOLELEMENTS; + typedef Array T_SURFELEMENTS; private: /// point coordinates @@ -192,32 +192,34 @@ namespace netgen void ClearSurfaceElements(); /// - DLL_HEADER void ClearVolumeElements() + DLL_HEADER void ClearVolumeElements() { volelements.SetSize(0); timestamp = NextTimeStamp(); } /// - DLL_HEADER void ClearSegments() + DLL_HEADER void ClearSegments() { segments.SetSize(0); timestamp = NextTimeStamp(); } - + /// bool TestOk () const; void SetAllocSize(int nnodes, int nsegs, int nsel, int nel); - + DLL_HEADER PointIndex AddPoint (const Point3d & p, int layer = 1); DLL_HEADER PointIndex AddPoint (const Point3d & p, int layer, POINTTYPE type); int GetNP () const { return points.Size(); } + // [[deprecated("Use Point(PointIndex) instead of int !")]] MeshPoint & Point(int i) { return points.Elem(i); } MeshPoint & Point(PointIndex pi) { return points[pi]; } + // [[deprecated("Use Point(PointIndex) instead of int !")]] const MeshPoint & Point(int i) const { return points.Get(i); } const MeshPoint & Point(PointIndex pi) const { return points[pi]; } @@ -231,16 +233,20 @@ namespace netgen DLL_HEADER SegmentIndex AddSegment (const Segment & s); void DeleteSegment (int segnr) { - segments.Elem(segnr)[0] = PointIndex::BASE-1; - segments.Elem(segnr)[1] = PointIndex::BASE-1; + segments.Elem(segnr)[0].Invalidate(); + segments.Elem(segnr)[1].Invalidate(); } + /* void FullDeleteSegment (int segnr) // von wem ist das ??? { segments.Delete(segnr-PointIndex::BASE); } + */ int GetNSeg () const { return segments.Size(); } + [[deprecated("Use LineSegment(SegmentIndex) instead of int !")]] Segment & LineSegment(int i) { return segments.Elem(i); } + [[deprecated("Use LineSegment(SegmentIndex) instead of int !")]] const Segment & LineSegment(int i) const { return segments.Get(i); } Segment & LineSegment(SegmentIndex si) { return segments[si]; } @@ -254,22 +260,29 @@ namespace netgen Array pointelements; // only via python interface DLL_HEADER SurfaceElementIndex AddSurfaceElement (const Element2d & el); + + [[deprecated("Use DeleteSurfaceElement(SurfaceElementIndex) instead of int !")]] void DeleteSurfaceElement (int eli) { surfelements.Elem(eli).Delete(); - surfelements.Elem(eli).PNum(1) = -1; - surfelements.Elem(eli).PNum(2) = -1; - surfelements.Elem(eli).PNum(3) = -1; + surfelements.Elem(eli).PNum(1).Invalidate(); + surfelements.Elem(eli).PNum(2).Invalidate(); + surfelements.Elem(eli).PNum(3).Invalidate(); timestamp = NextTimeStamp(); } void DeleteSurfaceElement (SurfaceElementIndex eli) { - DeleteSurfaceElement (int(eli)+1); + for (auto & p : surfelements[eli].PNums()) p.Invalidate(); + surfelements[eli].Delete(); + timestamp = NextTimeStamp(); } int GetNSE () const { return surfelements.Size(); } + + [[deprecated("Use SurfaceElement(SurfaceElementIndex) instead of int !")]] Element2d & SurfaceElement(int i) { return surfelements.Elem(i); } + [[deprecated("Use SurfaceElement(SurfaceElementIndex) instead of int !")]] const Element2d & SurfaceElement(int i) const { return surfelements.Get(i); } Element2d & SurfaceElement(SurfaceElementIndex i) { return surfelements[i]; } const Element2d & SurfaceElement(SurfaceElementIndex i) const { return surfelements[i]; } @@ -290,7 +303,9 @@ namespace netgen int GetNE () const { return volelements.Size(); } + // [[deprecated("Use VolumeElement(ElementIndex) instead of int !")]] Element & VolumeElement(int i) { return volelements.Elem(i); } + // [[deprecated("Use VolumeElement(ElementIndex) instead of int !")]] const Element & VolumeElement(int i) const { return volelements.Get(i); } Element & VolumeElement(ElementIndex i) { return volelements[i]; } const Element & VolumeElement(ElementIndex i) const { return volelements[i]; } @@ -305,9 +320,8 @@ namespace netgen ELEMENTTYPE ElementType (ElementIndex i) const { return (volelements[i].flags.fixed) ? FIXEDELEMENT : FREEELEMENT; } - const T_VOLELEMENTS & VolumeElements() const { return volelements; } - T_VOLELEMENTS & VolumeElements() { return volelements; } - + const auto & VolumeElements() const { return volelements; } + auto & VolumeElements() { return volelements; } /// DLL_HEADER double ElementError (int eli, const MeshingParameters & mp) const; @@ -317,17 +331,13 @@ namespace netgen /// void ClearLockedPoints (); - const Array & LockedPoints() const - { return lockedpoints; } + const auto & LockedPoints() const { return lockedpoints; } /// Returns number of domains int GetNDomains() const; - /// - int GetDimension() const - { return dimension; } - void SetDimension(int dim) - { dimension = dim; } + int GetDimension() const { return dimension; } + void SetDimension (int dim) { dimension = dim; } /// sets internal tables void CalcSurfacesOfNode (); @@ -573,10 +583,10 @@ namespace netgen /// int AddFaceDescriptor(const FaceDescriptor& fd) - { return facedecoding.Append(fd); } + { facedecoding.Append(fd); return facedecoding.Size(); } int AddEdgeDescriptor(const EdgeDescriptor & fd) - { return edgedecoding.Append(fd) - 1; } + { edgedecoding.Append(fd); return edgedecoding.Size() - 1; } /// DLL_HEADER void SetMaterial (int domnr, const string & mat); diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 062c3e99..9af43005 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -708,7 +708,7 @@ namespace netgen break; PrintMessage (5, nillegal, " illegal tets"); - optmesh.SplitImprove (mesh3d, OPT_LEGAL); + optmesh.SplitImprove (mesh3d, OPT_LEGAL); mesh3d.MarkIllegalElements(); // test optmesh.SwapImprove (mesh3d, OPT_LEGAL); diff --git a/libsrc/meshing/meshing3.cpp b/libsrc/meshing/meshing3.cpp index eb756460..5f12c32d 100644 --- a/libsrc/meshing/meshing3.cpp +++ b/libsrc/meshing/meshing3.cpp @@ -176,14 +176,14 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) NgProfiler::RegionTimer reg (meshing3_timer); - Array locpoints; // local points - Array locfaces; // local faces - Array pindex; // mapping from local to front point numbering - Array allowpoint; // point is allowd ? - Array findex; // mapping from local to front face numbering - //INDEX_2_HASHTABLE connectedpairs(100); // connecgted pairs for prism meshing + Array locpoints; // local points + Array locfaces; // local faces + Array pindex; // mapping from local to front point numbering + Array allowpoint; // point is allowd ? + Array findex; // mapping from local to front face numbering + //INDEX_2_HASHTABLE connectedpairs(100); // connecgted pairs for prism meshing - Array plainpoints; // points in reference coordinates + Array plainpoints; // points in reference coordinates Array delpoints, delfaces; // points and lines to be deleted Array locelements; // new generated elements @@ -206,9 +206,9 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) // for star-shaped domain meshing - Array grouppoints; + Array grouppoints; Array groupfaces; - Array grouppindex; + Array grouppindex; Array groupfindex; @@ -269,9 +269,9 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) } const MiniElement2d & bel = adfront->GetFace (baseelem); - const Point3d & p1 = adfront->GetPoint (bel.PNum(1)); - const Point3d & p2 = adfront->GetPoint (bel.PNum(2)); - const Point3d & p3 = adfront->GetPoint (bel.PNum(3)); + const Point3d & p1 = adfront->GetPoint (bel[0]); + const Point3d & p2 = adfront->GetPoint (bel[1]); + const Point3d & p3 = adfront->GetPoint (bel[2]); // (*testout) << endl << "base = " << bel << endl; @@ -302,9 +302,6 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) // (*testout) << "locfaces = " << endl << locfaces << endl; - int pi1 = pindex.Get(locfaces[0].PNum(1)); - int pi2 = pindex.Get(locfaces[0].PNum(2)); - int pi3 = pindex.Get(locfaces[0].PNum(3)); //loktestmode = 1; testmode = loktestmode; //changed @@ -316,6 +313,9 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) if (loktestmode) { (*testout) << "baseel = " << baseelem << ", ind = " << findex.Get(1) << endl; + int pi1 = pindex[locfaces[0].PNum(1)]; + int pi2 = pindex[locfaces[0].PNum(2)]; + int pi3 = pindex[locfaces[0].PNum(3)]; (*testout) << "pi = " << pi1 << ", " << pi2 << ", " << pi3 << endl; } @@ -369,7 +369,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) grouppindex, groupfindex); bool onlytri = 1; - for (i = 0; i < groupfaces.Size(); i++) + for (auto i : groupfaces.Range()) if (groupfaces[i].GetNP() != 3) onlytri = 0; @@ -441,31 +441,28 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) { // set transformatino to reference coordinates - if (locfaces.Get(1).GetNP() == 3) + if (locfaces[0].GetNP() == 3) { - trans.Set (locpoints.Get(locfaces.Get(1).PNumMod(1+rotind)), - locpoints.Get(locfaces.Get(1).PNumMod(2+rotind)), - locpoints.Get(locfaces.Get(1).PNumMod(3+rotind)), hshould); + trans.Set (locpoints[locfaces[0].PNumMod(1+rotind)], + locpoints[locfaces[0].PNumMod(2+rotind)], + locpoints[locfaces[0].PNumMod(3+rotind)], hshould); } else { - trans.Set (locpoints.Get(locfaces.Get(1).PNumMod(1+rotind)), - locpoints.Get(locfaces.Get(1).PNumMod(2+rotind)), - locpoints.Get(locfaces.Get(1).PNumMod(4+rotind)), hshould); + trans.Set (locpoints[locfaces[0].PNumMod(1+rotind)], + locpoints[locfaces[0].PNumMod(2+rotind)], + locpoints[locfaces[0].PNumMod(4+rotind)], hshould); } - trans.ToPlain (locpoints, plainpoints); + // trans.ToPlain (locpoints, plainpoints); - - for (i = 1; i <= allowpoint.Size(); i++) - { - if (plainpoints.Get(i).Z() > 0) - { - //if(loktestmode) - // (*testout) << "plainpoints.Get(i).Z() = " << plainpoints.Get(i).Z() << " > 0" << endl; - allowpoint.Elem(i) = 0; - } - } + plainpoints.SetSize (locpoints.Size()); + for (auto i : locpoints.Range()) + trans.ToPlain (locpoints[i], plainpoints[i]); + + for (auto i : allowpoint.Range()) + if (plainpoints[i].Z() > 0) + allowpoint[i] = false; stat.cnttrials++; @@ -508,7 +505,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) if (found) stat.cntsucc++; locpoints.SetSize (plainpoints.Size()); - for (i = oldnp+1; i <= plainpoints.Size(); i++) + for (int i = oldnp+1; i <= plainpoints.Size(); i++) trans.FromPlain (plainpoints.Elem(i), locpoints.Elem(i)); @@ -516,13 +513,13 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) // avoid meshing from large to small mesh-size if (uselocalh && found && stat.qualclass <= 3) { - for (i = 1; i <= locelements.Size(); i++) + for (int i = 1; i <= locelements.Size(); i++) { - Point3d pmin = locpoints.Get(locelements.Get(i).PNum(1)); + Point3d pmin = locpoints[locelements.Get(i).PNum(1)]; Point3d pmax = pmin; for (j = 2; j <= 4; j++) { - const Point3d & hp = locpoints.Get(locelements.Get(i).PNum(j)); + const Point3d & hp = locpoints[locelements.Get(i).PNum(j)]; pmin.SetToMin (hp); pmax.SetToMax (hp); } @@ -533,10 +530,10 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) } if (found) { - for (i = 1; i <= locelements.Size(); i++) - for (j = 1; j <= 4; j++) + for (int i = 1; i <= locelements.Size(); i++) + for (int j = 1; j <= 4; j++) { - const Point3d & hp = locpoints.Get(locelements.Get(i).PNum(j)); + const Point3d & hp = locpoints[locelements.Get(i).PNum(j)]; if (Dist (hp, pmid) > hinner) found = 0; } @@ -655,25 +652,25 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) pindex.SetSize(locpoints.Size()); - for (i = oldnp+1; i <= locpoints.Size(); i++) + for (int i = oldnp+1; i <= locpoints.Size(); i++) { PointIndex globind = mesh.AddPoint (locpoints.Get(i)); pindex.Elem(i) = adfront -> AddPoint (locpoints.Get(i), globind); } - for (i = 1; i <= locelements.Size(); i++) + for (int i = 1; i <= locelements.Size(); i++) { Point3d * hp1, * hp2, * hp3, * hp4; - hp1 = &locpoints.Elem(locelements.Get(i).PNum(1)); - hp2 = &locpoints.Elem(locelements.Get(i).PNum(2)); - hp3 = &locpoints.Elem(locelements.Get(i).PNum(3)); - hp4 = &locpoints.Elem(locelements.Get(i).PNum(4)); + hp1 = &locpoints[locelements.Get(i).PNum(1)]; + hp2 = &locpoints[locelements.Get(i).PNum(2)]; + hp3 = &locpoints[locelements.Get(i).PNum(3)]; + hp4 = &locpoints[locelements.Get(i).PNum(4)]; tetvol += (1.0 / 6.0) * ( Cross ( *hp2 - *hp1, *hp3 - *hp1) * (*hp4 - *hp1) ); for (j = 1; j <= locelements.Get(i).NP(); j++) locelements.Elem(i).PNum(j) = - adfront -> GetGlobalIndex (pindex.Get(locelements.Get(i).PNum(j))); + adfront -> GetGlobalIndex (pindex[locelements.Get(i).PNum(j)]); mesh.AddVolumeElement (locelements.Get(i)); stat.cntelem++; @@ -683,7 +680,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) { for (j = 1; j <= locfaces.Get(i).GetNP(); j++) locfaces.Elem(i).PNum(j) = - pindex.Get(locfaces.Get(i).PNum(j)); + pindex[locfaces.Get(i).PNum(j)]; // (*testout) << "add face " << locfaces.Get(i) << endl; adfront->AddFace (locfaces.Get(i)); } diff --git a/libsrc/meshing/meshing3.hpp b/libsrc/meshing/meshing3.hpp index 84e9ca6b..2acd209d 100644 --- a/libsrc/meshing/meshing3.hpp +++ b/libsrc/meshing/meshing3.hpp @@ -42,7 +42,8 @@ public: MESHING3_RESULT GenerateMesh (Mesh & mesh, const MeshingParameters & mp); /// - int ApplyRules (Array & lpoints, Array & allowpoint, + int ApplyRules (Array & lpoints, + Array & allowpoint, Array & lfaces, INDEX lfacesplit, INDEX_2_HASHTABLE & connectedpairs, Array & elements, diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index b4476e00..8de4947b 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -538,6 +538,9 @@ namespace netgen shape(3) = (1-p(0))* p(1) ; break; } + + default: + throw NgException ("illegal element type in GetShapeNew"); } } @@ -562,6 +565,8 @@ namespace netgen shape(3) = (1-p(0))* p(1) ; break; } + default: + throw NgException ("illegal element type in GetShapeNew"); } } @@ -641,6 +646,8 @@ namespace netgen dshape(3,1) = (1-p(0)); break; } + default: + throw NgException ("illegal element type in GetDShapeNew"); } } @@ -1739,12 +1746,14 @@ namespace netgen { 0, 0, 1, 1 }, { 0, 0, 0, 1 }, }; - + double * pp = NULL; switch (typ) { case TET: pp = &eltetqp[0][0]; break; case TET10: pp = &eltet10qp[ip-1][0]; break; + default: + throw NgException ("illegal element shape in GetIntegrationPoint"); } p(0) = pp[0]; diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index dfdf889a..a38bcad0 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -125,6 +125,8 @@ namespace netgen PointIndex operator-- (int) { PointIndex hi(*this); i--; return hi; } PointIndex operator++ () { i++; return *this; } PointIndex operator-- () { i--; return *this; } + void Invalidate() { i = PointIndex::BASE-1; } + bool IsValid() const { return i != PointIndex::BASE-1; } #ifdef BASE0 enum { BASE = 0 }; #else @@ -134,7 +136,7 @@ namespace netgen inline istream & operator>> (istream & ist, PointIndex & pi) { - int i; ist >> i; pi = i; return ist; + int i; ist >> i; pi = PointIndex(i); return ist; } inline ostream & operator<< (ostream & ost, const PointIndex & pi) @@ -154,8 +156,10 @@ namespace netgen ElementIndex & operator= (const ElementIndex & ai) { i = ai.i; return *this; } ElementIndex & operator= (int ai) { i = ai; return *this; } operator int () const { return i; } - ElementIndex & operator++ (int) { i++; return *this; } - ElementIndex & operator-- (int) { i--; return *this; } + ElementIndex operator++ (int) { return ElementIndex(i++); } + ElementIndex operator-- (int) { return ElementIndex(i--); } + ElementIndex & operator++ () { ++i; return *this; } + ElementIndex & operator-- () { --i; return *this; } }; inline istream & operator>> (istream & ist, ElementIndex & pi) @@ -179,9 +183,11 @@ namespace netgen { i = ai.i; return *this; } SurfaceElementIndex & operator= (int ai) { i = ai; return *this; } operator int () const { return i; } - SurfaceElementIndex & operator++ (int) { i++; return *this; } + SurfaceElementIndex operator++ (int) { SurfaceElementIndex hi(*this); i++; return hi; } + SurfaceElementIndex operator-- (int) { SurfaceElementIndex hi(*this); i--; return hi; } + SurfaceElementIndex & operator++ () { ++i; return *this; } + SurfaceElementIndex & operator-- () { --i; return *this; } SurfaceElementIndex & operator+= (int inc) { i+=inc; return *this; } - SurfaceElementIndex & operator-- (int) { i--; return *this; } }; inline istream & operator>> (istream & ist, SurfaceElementIndex & pi) @@ -344,6 +350,7 @@ namespace netgen default: PrintSysError ("Element2d::SetType, illegal type ", int(typ)); } + is_curved = (np >= 4); } /// int GetNP() const { return np; } @@ -386,6 +393,8 @@ namespace netgen FlatArray PNums () const { return FlatArray (np, &pnum[0]); } + FlatArray PNums () + { return FlatArray (np, &pnum[0]); } /// PointIndex & PNum (int i) { return pnum[i-1]; } @@ -470,8 +479,13 @@ namespace netgen int pi, Vec2d & dir, double & dd) const; - - void Delete () { deleted = 1; pnum[0] = pnum[1] = pnum[2] = pnum[3] = PointIndex::BASE-1; } + + void Delete () + { + deleted = 1; + for (PointIndex & p : pnum) p.Invalidate(); + } + bool IsDeleted () const { #ifdef DEBUG diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 8ac969c1..10a19466 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -304,8 +304,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) - ExportArray(m); - ExportArray(m); + ExportArray(m); + ExportArray(m); ExportArray(m); ExportArray(m); ExportArray(m); @@ -384,11 +384,11 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def_property("dim", &Mesh::GetDimension, &Mesh::SetDimension) .def("Elements3D", - static_cast&(Mesh::*)()> (&Mesh::VolumeElements), + static_cast&(Mesh::*)()> (&Mesh::VolumeElements), py::return_value_policy::reference) .def("Elements2D", - static_cast&(Mesh::*)()> (&Mesh::SurfaceElements), + static_cast&(Mesh::*)()> (&Mesh::SurfaceElements), py::return_value_policy::reference) .def("Elements1D", diff --git a/libsrc/meshing/refine.cpp b/libsrc/meshing/refine.cpp index 686b7190..4e8a5ae7 100644 --- a/libsrc/meshing/refine.cpp +++ b/libsrc/meshing/refine.cpp @@ -383,11 +383,11 @@ namespace netgen swap (pnums.Elem(3), pnums.Elem(4)); for (int j = 0; j < 6; j++) - { - INDEX_2 i2; - i2.I1() = pnums.Get(betw[j][0]); - i2.I2() = pnums.Get(betw[j][1]); - i2.Sort(); + { + PointIndex pi1 = pnums.Get(betw[j][0]); + PointIndex pi2 = pnums.Get(betw[j][1]); + INDEX_2 i2 (pi1, pi2); + i2.Sort(); /* if (between.Used(i2)) @@ -405,8 +405,8 @@ namespace netgen if (!pointset[pinew]) { pointset[pinew] = true; - mesh.Point(pinew) = Center(mesh.Point(i2.I1()), - mesh.Point(i2.I2())); + mesh.Point(pinew) = Center(mesh.Point(pi1), + mesh.Point(pi2)); } } diff --git a/libsrc/meshing/ruler2.cpp b/libsrc/meshing/ruler2.cpp index f1c84ee4..a6d5410f 100644 --- a/libsrc/meshing/ruler2.cpp +++ b/libsrc/meshing/ruler2.cpp @@ -593,8 +593,9 @@ namespace netgen Point2d np = rule->GetPoint(i); np.X() += newu (2 * (i-oldnp) - 2); np.Y() += newu (2 * (i-oldnp) - 1); - - pmap.Elem(i) = lpoints.Append (np); + + lpoints.Append (np); + pmap.Elem(i) = lpoints.Size(); } } diff --git a/libsrc/meshing/ruler3.cpp b/libsrc/meshing/ruler3.cpp index c83b1db5..e33142bf 100644 --- a/libsrc/meshing/ruler3.cpp +++ b/libsrc/meshing/ruler3.cpp @@ -8,29 +8,29 @@ extern double minother; extern double minwithoutother; -static double CalcElementBadness (const Array & points, - const Element & elem) + static double CalcElementBadness (const Array & points, + const Element & elem) { double vol, l, l4, l5, l6; if (elem.GetNP() != 4) { if (elem.GetNP() == 5) { - double z = points.Get(elem.PNum(5)).Z(); + double z = points[elem.PNum(5)].Z(); if (z > -1e-8) return 1e8; return (-1 / z) - z; // - 2; } return 0; } - Vec3d v1 = points.Get(elem.PNum(2)) - points.Get(elem.PNum(1)); - Vec3d v2 = points.Get(elem.PNum(3)) - points.Get(elem.PNum(1)); - Vec3d v3 = points.Get(elem.PNum(4)) - points.Get(elem.PNum(1)); + Vec3d v1 = points[elem.PNum(2)] - points[elem.PNum(1)]; + Vec3d v2 = points[elem.PNum(3)] - points[elem.PNum(1)]; + Vec3d v3 = points[elem.PNum(4)] - points[elem.PNum(1)]; vol = - (Cross (v1, v2) * v3); - l4 = Dist (points.Get(elem.PNum(2)), points.Get(elem.PNum(3))); - l5 = Dist (points.Get(elem.PNum(2)), points.Get(elem.PNum(4))); - l6 = Dist (points.Get(elem.PNum(3)), points.Get(elem.PNum(4))); + l4 = Dist (points[elem.PNum(2)], points[elem.PNum(3)]); + l5 = Dist (points[elem.PNum(2)], points[elem.PNum(4)]); + l6 = Dist (points[elem.PNum(3)], points[elem.PNum(4)]); l = v1.Length() + v2.Length() + v3.Length() + l4 + l5 + l6; @@ -49,8 +49,8 @@ static double CalcElementBadness (const Array & points, int Meshing3 :: ApplyRules ( - Array & lpoints, // in: local points, out: old+new local points - Array & allowpoint, // in: 2 .. it is allowed to use pointi, 1..will be allowed later, 0..no means + Array & lpoints, // in: local points, out: old+new local points + Array & allowpoint, // in: 2 .. it is allowed to use pointi, 1..will be allowed later, 0..no means Array & lfaces, // in: local faces, out: old+new local faces INDEX lfacesplit, // for local faces in outer radius INDEX_2_HASHTABLE & connectedpairs, // connected pairs for prism-meshing @@ -65,25 +65,23 @@ int Meshing3 :: ApplyRules { NgProfiler::RegionTimer regtot(97); - int i, j, k, ri, nfok, npok, incnpok, refpi, locpi, locfi, locfr; - float hf, err, minerr, teterr, minteterr; + float err, minerr, teterr, minteterr; char ok, found, hc; - vnetrule * rule; + // vnetrule * rule; Vector oldu, newu, newu1, newu2, allp; Vec3d ui; Point3d np; - int oldnp, noldlp, noldlf; const MiniElement2d * locface = NULL; int loktestmode; - Array pused; // point is already mapped - Array fused; // face is already mapped - Array pmap; // map of reference point to local point - Array pfixed; // point mapped by face-map - Array fmapi; // face in reference is mapped to face nr ... - Array fmapr; // face in reference is rotated to map - Array transfreezone; // transformed free-zone + Array pused; // point is already mapped, number of uses + Array fused; // face is already mapped + Array pmap; // map of reference point to local point + Array pfixed; // point mapped by face-map + Array fmapi; // face in reference is mapped to face nr ... + Array fmapr; // face in reference is rotated to map + Array transfreezone; // transformed free-zone INDEX_2_CLOSED_HASHTABLE ledges(100); // edges in local environment Array tempnewpoints; @@ -108,9 +106,10 @@ int Meshing3 :: ApplyRules fnearness.SetSize (lfacesplit); pnearness = INT_MAX/10; - for (j = 0; j < lfaces[0].GetNP(); j++) - pnearness[lfaces[0][j]] = 0; + for (PointIndex pi : lfaces[0].PNums()) + pnearness[pi] = 0; + NgProfiler::RegionTimer reg2(98); NgProfiler::StartTimer (90); @@ -118,24 +117,24 @@ int Meshing3 :: ApplyRules for (int loop = 0; loop < 2; loop++) { - for (i = 0; i < lfacesplit; i++) + for (int i = 0; i < lfacesplit; i++) { const MiniElement2d & hface = lfaces[i]; int minn = INT_MAX-1; - for (j = 0; j < hface.GetNP(); j++) + for (PointIndex pi : hface.PNums()) { - int hi = pnearness[hface[j]]; + int hi = pnearness[pi]; if (hi < minn) minn = hi; } if (minn < INT_MAX/10) - for (j = 0; j < hface.GetNP(); j++) - if (pnearness[hface[j]] > minn+1) - pnearness[hface[j]] = minn+1; + for (PointIndex pi : hface.PNums()) + if (pnearness[pi] > minn+1) + pnearness[pi] = minn+1; } - for (i = 1; i <= connectedpairs.GetNBags(); i++) - for (j = 1; j <= connectedpairs.GetBagSize(i); j++) + for (int i = 1; i <= connectedpairs.GetNBags(); i++) + for (int j = 1; j <= connectedpairs.GetBagSize(i); j++) { INDEX_2 edge; int val; @@ -147,14 +146,13 @@ int Meshing3 :: ApplyRules if (pnearness[edge.I2()] > pnearness[edge.I1()] + 1) pnearness[edge.I2()] = pnearness[edge.I1()] + 1; } - } - for (i = 0; i < fnearness.Size(); i++) + for (int i : fnearness.Range()) { int sum = 0; - for (j = 0; j < lfaces[i].GetNP(); j++) - sum += pnearness[lfaces[i][j]]; + for (PointIndex pi : lfaces[i].PNums()) + sum += pnearness[pi]; fnearness[i] = sum; } @@ -165,34 +163,35 @@ int Meshing3 :: ApplyRules // find bounding boxes of faces triboxes.SetSize (lfaces.Size()); - for (i = 0; i < lfaces.Size(); i++) + // for (int i = 0; i < lfaces.Size(); i++) + for (auto i : lfaces.Range()) { const MiniElement2d & face = lfaces[i]; - triboxes[i].SetPoint (lpoints.Get(face[0])); - for (j = 1; j < face.GetNP(); j++) - triboxes[i].AddPoint (lpoints.Get(face[j])); + triboxes[i].SetPoint (lpoints[face[0]]); + for (int j = 1; j < face.GetNP(); j++) + triboxes[i].AddPoint (lpoints[face[j]]); } NgProfiler::StopTimer (91); NgProfiler::StartTimer (92); - bool useedges = 0; - for (ri = 0; ri < rules.Size(); ri++) - if (rules[ri]->GetNEd()) useedges = 1; + bool useedges = false; + for (int ri = 0; ri < rules.Size(); ri++) + if (rules[ri]->GetNEd()) useedges = true; if (useedges) { ledges.SetSize (5 * lfacesplit); - for (j = 0; j < lfacesplit; j++) + for (int j = 0; j < lfacesplit; j++) // if (fnearness[j] <= 5) { const MiniElement2d & face = lfaces[j]; int newp, oldp; newp = face[face.GetNP()-1]; - for (k = 0; k < face.GetNP(); k++) + for (int k = 0; k < face.GetNP(); k++) { oldp = newp; newp = face[k]; @@ -223,8 +222,8 @@ int Meshing3 :: ApplyRules // check each rule: - for (ri = 1; ri <= rules.Size(); ri++) - { + for (int ri = 1; ri <= rules.Size(); ri++) + { int base = (lfaces[0].GetNP() == 3) ? 100 : 200; NgProfiler::RegionTimer regx1(base); NgProfiler::RegionTimer regx(base+ri); @@ -232,7 +231,7 @@ int Meshing3 :: ApplyRules // sprintf (problems.Elem(ri), ""); *problems.Elem(ri) = '\0'; - rule = rules.Get(ri); + vnetrule * rule = rules.Get(ri); if (rule->GetNP(1) != lfaces[0].GetNP()) continue; @@ -260,21 +259,21 @@ int Meshing3 :: ApplyRules fused = 0; pused = 0; - pmap = 0; + for (auto & p : pmap) p.Invalidate(); fmapi = 0; - for (i = 1; i <= fmapr.Size(); i++) - fmapr.Set(i, rule->GetNP(i)); + + for (int i : fmapr.Range()) + fmapr[i] = rule->GetNP(i+1); fused[0] = 1; fmapi[0] = 1; fmapr[0] = rotind1; - - for (j = 1; j <= lfaces.Get(1).GetNP(); j++) + for (int j = 1; j <= lfaces[0].GetNP(); j++) { - locpi = lfaces[0].PNumMod (j+rotind1); + PointIndex locpi = lfaces[0].PNumMod (j+rotind1); pmap.Set (rule->GetPointNr (1, j), locpi); - pused.Elem(locpi)++; + pused[locpi]++; } /* @@ -282,7 +281,7 @@ int Meshing3 :: ApplyRules nfok .. first nfok-1 faces are mapped properly */ - nfok = 2; + int nfok = 2; NgProfiler::RegionTimer regfa(300); NgProfiler::RegionTimer regx2(base+50+ri); while (nfok >= 2) @@ -293,8 +292,8 @@ int Meshing3 :: ApplyRules // not all faces mapped ok = 0; - locfi = fmapi.Get(nfok); - locfr = fmapr.Get(nfok); + int locfi = fmapi.Get(nfok); + int locfr = fmapr.Get(nfok); int actfnp = rule->GetNP(nfok); @@ -326,28 +325,27 @@ int Meshing3 :: ApplyRules // reference point already mapped differently ? - for (j = 1; j <= actfnp && ok; j++) + for (int j = 1; j <= actfnp && ok; j++) { - locpi = pmap.Get(rule->GetPointNr (nfok, j)); - - if (locpi && locpi != locface->PNumMod(j+locfr)) + PointIndex locpi = pmap.Get(rule->GetPointNr (nfok, j)); + if (locpi.IsValid() && locpi != locface->PNumMod(j+locfr)) ok = 0; } // local point already used or point outside tolerance ? - for (j = 1; j <= actfnp && ok; j++) + for (int j = 1; j <= actfnp && ok; j++) { - refpi = rule->GetPointNr (nfok, j); + int refpi = rule->GetPointNr (nfok, j); - if (pmap.Get(refpi) == 0) + if (!pmap.Get(refpi).IsValid()) { - locpi = locface->PNumMod (j + locfr); + PointIndex locpi = locface->PNumMod (j + locfr); - if (pused.Get(locpi)) + if (pused[locpi]) ok = 0; else { - const Point3d & lp = lpoints.Get(locpi); + const Point3d & lp = lpoints[locpi]; const Point3d & rp = rule->GetPoint(refpi); if ( Dist2 (lp, rp) * rule->PointDistFactor(refpi) > minerr) @@ -370,16 +368,16 @@ int Meshing3 :: ApplyRules fmapr.Set (nfok, locfr); fused.Set (locfi, 1); - for (j = 1; j <= rule->GetNP (nfok); j++) + for (int j = 1; j <= rule->GetNP (nfok); j++) { - locpi = locface->PNumMod(j+locfr); + PointIndex locpi = locface->PNumMod(j+locfr); if (rule->GetPointNr (nfok, j) <= 3 && pmap.Get(rule->GetPointNr(nfok, j)) != locpi) (*testout) << "change face1 point, mark1" << endl; pmap.Set(rule->GetPointNr (nfok, j), locpi); - pused.Elem(locpi)++; + pused[locpi]++; } nfok++; @@ -392,14 +390,15 @@ int Meshing3 :: ApplyRules nfok--; fused.Set (fmapi.Get(nfok), 0); - for (j = 1; j <= rule->GetNP (nfok); j++) + for (int j = 1; j <= rule->GetNP (nfok); j++) { - refpi = rule->GetPointNr (nfok, j); - pused.Elem(pmap.Get(refpi))--; - - if (pused.Get(pmap.Get(refpi)) == 0) + int refpi = rule->GetPointNr (nfok, j); + pused[pmap.Get(refpi)]--; + + if (pused[pmap.Get(refpi)] == 0) { - pmap.Set(refpi, 0); + // pmap.Set(refpi, 0); + pmap.Elem(refpi).Invalidate(); } } } @@ -418,14 +417,18 @@ int Meshing3 :: ApplyRules (*testout) << "Faces Ok" << endl; sprintf (problems.Elem(ri), "Faces Ok"); } - - npok = 1; - incnpok = 1; + + int npok = 1; + int incnpok = 1; pfixed.SetSize (pmap.Size()); - for (i = 1; i <= pmap.Size(); i++) + /* + for (int i = 1; i <= pmap.Size(); i++) pfixed.Set(i, (pmap.Get(i) != 0) ); - + */ + for (int i : pmap.Range()) + pfixed[i] = pmap[i].IsValid(); + while (npok >= 1) { @@ -444,31 +447,31 @@ int Meshing3 :: ApplyRules else { - locpi = pmap.Elem(npok); + PointIndex locpi = pmap.Elem(npok); ok = 0; - if (locpi) - pused.Elem(locpi)--; + if (locpi.IsValid()) + pused[locpi]--; - while (!ok && locpi < lpoints.Size()) + while (!ok && locpi < lpoints.Size()-1+PointIndex::BASE) { ok = 1; locpi++; - if (pused.Get(locpi) || - pnearness.Get(locpi) > rule->GetPNearness(npok)) + if (pused[locpi] || + pnearness[locpi] > rule->GetPNearness(npok)) { ok = 0; } - else if (allowpoint.Get(locpi) != 2) + else if (allowpoint[locpi] != 2) { ok = 0; - if (allowpoint.Get(locpi) == 1) + if (allowpoint[locpi] == 1) impossible = 0; } else { - const Point3d & lp = lpoints.Get(locpi); + const Point3d & lp = lpoints[locpi]; const Point3d & rp = rule->GetPoint(npok); if ( Dist2 (lp, rp) * rule->PointDistFactor(npok) > minerr) @@ -487,7 +490,7 @@ int Meshing3 :: ApplyRules if (npok <= 3) (*testout) << "set face1 point, mark3" << endl; - pused.Elem(locpi)++; + pused[locpi]++; npok++; incnpok = 1; } @@ -495,7 +498,8 @@ int Meshing3 :: ApplyRules else { - pmap.Set (npok, 0); + // pmap.Set (npok, 0); + pmap.Elem(npok).Invalidate(); if (npok <= 3) (*testout) << "set face1 point, mark4" << endl; @@ -516,8 +520,8 @@ int Meshing3 :: ApplyRules if (loktestmode) { (*testout) << "Mapping found!!: Rule " << rule->Name() << endl; - for (i = 1; i <= pmap.Size(); i++) - (*testout) << pmap.Get(i) << " "; + for (auto pi : pmap) + (*testout) << pi << " "; (*testout) << endl; sprintf (problems.Elem(ri), "mapping found"); (*testout) << rule->GetNP(1) << " = " << lfaces.Get(1).GetNP() << endl; @@ -527,7 +531,7 @@ int Meshing3 :: ApplyRules // check mapedges: - for (i = 1; i <= rule->GetNEd(); i++) + for (int i = 1; i <= rule->GetNEd(); i++) { INDEX_2 in2(pmap.Get(rule->GetEdge(i).i1), pmap.Get(rule->GetEdge(i).i2)); @@ -537,12 +541,12 @@ int Meshing3 :: ApplyRules // check prism edges: - for (i = 1; i <= rule->GetNE(); i++) + for (int i = 1; i <= rule->GetNE(); i++) { const Element & el = rule->GetElement (i); if (el.GetType() == PRISM) { - for (j = 1; j <= 3; j++) + for (int j = 1; j <= 3; j++) { INDEX_2 in2(pmap.Get(el.PNum(j)), pmap.Get(el.PNum(j+3))); @@ -554,7 +558,7 @@ int Meshing3 :: ApplyRules { if (loktestmode) (*testout) << "map pyramid, rule = " << rule->Name() << endl; - for (j = 1; j <= 2; j++) + for (int j = 1; j <= 2; j++) { INDEX_2 in2; if (j == 1) @@ -581,7 +585,7 @@ int Meshing3 :: ApplyRules - for (i = rule->GetNOldF() + 1; i <= rule->GetNF(); i++) + for (int i = rule->GetNOldF() + 1; i <= rule->GetNF(); i++) fmapi.Set(i, 0); @@ -599,9 +603,9 @@ int Meshing3 :: ApplyRules newu.SetSize (3 * (rule->GetNP() - rule->GetNOldP())); allp.SetSize (3 * rule->GetNP()); - for (i = 1; i <= rule->GetNOldP(); i++) + for (int i = 1; i <= rule->GetNOldP(); i++) { - const Point3d & lp = lpoints.Get(pmap.Get(i)); + const Point3d & lp = lpoints[pmap.Get(i)]; const Point3d & rp = rule->GetPoint(i); oldu (3*i-3) = lp.X()-rp.X(); oldu (3*i-2) = lp.Y()-rp.Y(); @@ -620,7 +624,7 @@ int Meshing3 :: ApplyRules // int idiff = 3 * (rule->GetNP()-rule->GetNOldP()); int idiff = 3 * rule->GetNOldP(); - for (i = rule->GetNOldP()+1; i <= rule->GetNP(); i++) + for (int i = rule->GetNOldP()+1; i <= rule->GetNP(); i++) { const Point3d & rp = rule->GetPoint(i); allp (3*i-3) = rp.X() + newu(3*i-3 - idiff); @@ -644,14 +648,14 @@ int Meshing3 :: ApplyRules { const Array & fz = rule->GetTransFreeZone(); (*testout) << "Freezone: " << endl; - for (i = 1; i <= fz.Size(); i++) + for (int i = 1; i <= fz.Size(); i++) (*testout) << fz.Get(i) << endl; } // check freezone: - for (i = 1; i <= lpoints.Size(); i++) + for (int i = 1; i <= lpoints.Size(); i++) { if ( !pused.Get(i) ) { @@ -675,7 +679,7 @@ int Meshing3 :: ApplyRules } } - for (i = 1; i <= lfaces.Size() && ok; i++) + for (int i = 1; i <= lfaces.Size() && ok; i++) { static Array lpi(4); @@ -692,7 +696,7 @@ int Meshing3 :: ApplyRules for (li = 1; li <= lfacei.GetNP(); li++) { int lpii = 0; - int pi = lfacei.PNum(li); + PointIndex pi = lfacei.PNum(li); for (lj = 1; lj <= rule->GetNOldP(); lj++) if (pmap.Get(lj) == pi) lpii = lj; @@ -704,19 +708,19 @@ int Meshing3 :: ApplyRules { triin = rule->IsTriangleInFreeZone ( - lpoints.Get(lfacei.PNum(1)), - lpoints.Get(lfacei.PNum(2)), - lpoints.Get(lfacei.PNum(3)), lpi, 1 + lpoints[lfacei.PNum(1)], + lpoints[lfacei.PNum(2)], + lpoints[lfacei.PNum(3)], lpi, 1 ); } else { triin = rule->IsQuadInFreeZone ( - lpoints.Get(lfacei.PNum(1)), - lpoints.Get(lfacei.PNum(2)), - lpoints.Get(lfacei.PNum(3)), - lpoints.Get(lfacei.PNum(4)), + lpoints[lfacei.PNum(1)], + lpoints[lfacei.PNum(2)], + lpoints[lfacei.PNum(3)], + lpoints[lfacei.PNum(4)], lpi, 1 ); } @@ -741,7 +745,7 @@ int Meshing3 :: ApplyRules << lfaces.Get(i).PNum(3) << " - " << lfaces.Get(i).PNum(4) << endl; for (int lj = 1; lj <= lfaces.Get(i).GetNP(); lj++) - (*testout) << lpoints.Get(lfaces.Get(i).PNum(lj)) << " "; + (*testout) << lpoints[lfaces.Get(i).PNum(lj)] << " "; (*testout) << endl; @@ -759,9 +763,9 @@ int Meshing3 :: ApplyRules << lfacei.PNum(2) << " - " << lfacei.PNum(3) << ", or " - << lpoints.Get(lfacei.PNum(1)) << " - " - << lpoints.Get(lfacei.PNum(2)) << " - " - << lpoints.Get(lfacei.PNum(3)) + << lpoints[lfacei.PNum(1)] << " - " + << lpoints[lfacei.PNum(2)] << " - " + << lpoints[lfacei.PNum(3)] << endl; (*testout) << "lpi = " << lpi.Get(1) << ", " << lpi.Get(2) << ", " << lpi.Get(3) << endl; @@ -773,10 +777,10 @@ int Meshing3 :: ApplyRules << lfacei.PNum(3) << " - " << lfacei.PNum(4) << ", or " - << lpoints.Get(lfacei.PNum(1)) << " - " - << lpoints.Get(lfacei.PNum(2)) << " - " - << lpoints.Get(lfacei.PNum(3)) << " - " - << lpoints.Get(lfacei.PNum(4)) + << lpoints[lfacei.PNum(1)] << " - " + << lpoints[lfacei.PNum(2)] << " - " + << lpoints[lfacei.PNum(3)] << " - " + << lpoints[lfacei.PNum(4)] << endl; sprintf (problems.Elem(ri), "triangle (%d, %d, %d) in Freezone", @@ -786,13 +790,13 @@ int Meshing3 :: ApplyRules } hc = 0; - for (k = rule->GetNOldF() + 1; k <= rule->GetNF(); k++) + for (int k = rule->GetNOldF() + 1; k <= rule->GetNF(); k++) { if (rule->GetPointNr(k, 1) <= rule->GetNOldP() && rule->GetPointNr(k, 2) <= rule->GetNOldP() && rule->GetPointNr(k, 3) <= rule->GetNOldP()) { - for (j = 1; j <= 3; j++) + for (int j = 1; j <= 3; j++) if (lfaces.Get(i).PNumMod(j ) == pmap.Get(rule->GetPointNr(k, 1)) && lfaces.Get(i).PNumMod(j+1) == pmap.Get(rule->GetPointNr(k, 3)) && lfaces.Get(i).PNumMod(j+2) == pmap.Get(rule->GetPointNr(k, 2))) @@ -839,9 +843,9 @@ int Meshing3 :: ApplyRules if (ok) { err = 0; - for (i = 1; i <= rule->GetNOldP(); i++) + for (int i = 1; i <= rule->GetNOldP(); i++) { - hf = rule->CalcPointDist (i, lpoints.Get(pmap.Get(i))); + double hf = rule->CalcPointDist (i, lpoints[pmap.Get(i)]); if (hf > err) err = hf; } @@ -856,40 +860,37 @@ int Meshing3 :: ApplyRules // newu = rule->GetOldUToNewU() * oldu; // set new points: + int oldnp = rule->GetNOldP(); + int noldlp = lpoints.Size(); + int noldlf = lfaces.Size(); - oldnp = rule->GetNOldP(); - noldlp = lpoints.Size(); - noldlf = lfaces.Size(); - - - for (i = oldnp + 1; i <= rule->GetNP(); i++) + for (int i = oldnp + 1; i <= rule->GetNP(); i++) { np = rule->GetPoint(i); np.X() += newu (3 * (i-oldnp) - 3); np.Y() += newu (3 * (i-oldnp) - 2); np.Z() += newu (3 * (i-oldnp) - 1); - - pmap.Elem(i) = lpoints.Append (np); + lpoints.Append (np); + pmap.Elem(i) = lpoints.Size()-1+PointIndex::BASE; } // Set new Faces: - for (i = rule->GetNOldF() + 1; i <= rule->GetNF(); i++) + for (int i = rule->GetNOldF() + 1; i <= rule->GetNF(); i++) if (!fmapi.Get(i)) { MiniElement2d nface(rule->GetNP(i)); - for (j = 1; j <= nface.GetNP(); j++) + for (int j = 1; j <= nface.GetNP(); j++) nface.PNum(j) = pmap.Get(rule->GetPointNr (i, j)); lfaces.Append (nface); } - // Delete old Faces: - for (i = 1; i <= rule->GetNDelF(); i++) + for (int i = 1; i <= rule->GetNDelF(); i++) delfaces.Append (fmapi.Get(rule->GetDelFace(i))); - for (i = rule->GetNOldF()+1; i <= rule->GetNF(); i++) + for (int i = rule->GetNOldF()+1; i <= rule->GetNF(); i++) if (fmapi.Get(i)) { delfaces.Append (fmapi.Get(i)); @@ -898,17 +899,17 @@ int Meshing3 :: ApplyRules // check orientation - for (i = 1; i <= rule->GetNO() && ok; i++) + for (int i = 1; i <= rule->GetNO() && ok; i++) { const fourint * fouri; fouri = &rule->GetOrientation(i); - Vec3d v1 (lpoints.Get(pmap.Get(fouri->i1)), - lpoints.Get(pmap.Get(fouri->i2))); - Vec3d v2 (lpoints.Get(pmap.Get(fouri->i1)), - lpoints.Get(pmap.Get(fouri->i3))); - Vec3d v3 (lpoints.Get(pmap.Get(fouri->i1)), - lpoints.Get(pmap.Get(fouri->i4))); + Vec3d v1 (lpoints[pmap.Get(fouri->i1)], + lpoints[pmap.Get(fouri->i2)]); + Vec3d v2 (lpoints[pmap.Get(fouri->i1)], + lpoints[pmap.Get(fouri->i3)]); + Vec3d v3 (lpoints[pmap.Get(fouri->i1)], + lpoints[pmap.Get(fouri->i4)]); Vec3d n; Cross (v1, v2, n); @@ -927,7 +928,7 @@ int Meshing3 :: ApplyRules // new points in free-zone ? - for (i = rule->GetNOldP() + 1; i <= rule->GetNP() && ok; i++) + for (int i = rule->GetNOldP() + 1; i <= rule->GetNP() && ok; i++) if (!rule->IsInFreeZone (lpoints.Get(pmap.Get(i)))) { if (loktestmode) @@ -942,10 +943,10 @@ int Meshing3 :: ApplyRules // insert new elements - for (i = 1; i <= rule->GetNE(); i++) + for (int i = 1; i <= rule->GetNE(); i++) { elements.Append (rule->GetElement(i)); - for (j = 1; j <= elements.Get(i).NP(); j++) + for (int j = 1; j <= elements.Get(i).NP(); j++) elements.Elem(i).PNum(j) = pmap.Get(elements.Get(i).PNum(j)); } @@ -953,9 +954,9 @@ int Meshing3 :: ApplyRules // Calculate Element badness teterr = 0; - for (i = 1; i <= elements.Size(); i++) + for (int i = 1; i <= elements.Size(); i++) { - hf = CalcElementBadness (lpoints, elements.Get(i)); + double hf = CalcElementBadness (lpoints, elements.Get(i)); if (hf > teterr) teterr = hf; } @@ -978,29 +979,29 @@ int Meshing3 :: ApplyRules double oldlen = 0; double newlen = 0; - for (i = 1; i <= rule->GetNDelF(); i++) + for (int i = 1; i <= rule->GetNDelF(); i++) { const Element2d & face = rule->GetFace (rule->GetDelFace(i)); - for (j = 1; j <= 3; j++) + for (int j = 1; j <= 3; j++) { const Point3d & p1 = - lpoints.Get(pmap.Get(face.PNumMod(j))); + lpoints[pmap.Get(face.PNumMod(j))]; const Point3d & p2 = - lpoints.Get(pmap.Get(face.PNumMod(j+1))); + lpoints[pmap.Get(face.PNumMod(j+1))]; oldlen += Dist(p1, p2); } } - for (i = rule->GetNOldF()+1; i <= rule->GetNF(); i++) + for (int i = rule->GetNOldF()+1; i <= rule->GetNF(); i++) { const Element2d & face = rule->GetFace (i); - for (j = 1; j <= 3; j++) + for (int j = 1; j <= 3; j++) { const Point3d & p1 = - lpoints.Get(pmap.Get(face.PNumMod(j))); + lpoints[pmap.Get(face.PNumMod(j))]; const Point3d & p2 = - lpoints.Get(pmap.Get(face.PNumMod(j+1))); + lpoints[pmap.Get(face.PNumMod(j+1))]; newlen += Dist(p1, p2); } } @@ -1057,7 +1058,7 @@ int Meshing3 :: ApplyRules if (testmode) { - for (i = 1; i <= rule->GetNOldP(); i++) + for (int i = 1; i <= rule->GetNOldP(); i++) { (*testout) << "P" << i << ": Ref: " << rule->GetPoint (i) << " is: " @@ -1066,19 +1067,19 @@ int Meshing3 :: ApplyRules } tempnewpoints.SetSize (0); - for (i = noldlp+1; i <= lpoints.Size(); i++) + for (int i = noldlp+1; i <= lpoints.Size(); i++) tempnewpoints.Append (lpoints.Get(i)); tempnewfaces.SetSize (0); - for (i = noldlf+1; i <= lfaces.Size(); i++) + for (int i = noldlf+1; i <= lfaces.Size(); i++) tempnewfaces.Append (lfaces.Get(i)); tempdelfaces.SetSize (0); - for (i = 1; i <= delfaces.Size(); i++) + for (int i = 1; i <= delfaces.Size(); i++) tempdelfaces.Append (delfaces.Get(i)); tempelements.SetSize (0); - for (i = 1; i <= elements.Size(); i++) + for (int i = 1; i <= elements.Size(); i++) tempelements.Append (elements.Get(i)); } @@ -1096,15 +1097,13 @@ int Meshing3 :: ApplyRules nfok = rule->GetNOldF(); - for (j = 1; j <= rule->GetNP (nfok); j++) + for (int j = 1; j <= rule->GetNP (nfok); j++) { - refpi = rule->GetPointNr (nfok, j); - pused.Elem(pmap.Get(refpi))--; + int refpi = rule->GetPointNr (nfok, j); + pused[pmap.Get(refpi)]--; - if (pused.Get(pmap.Get(refpi)) == 0) - { - pmap.Set(refpi, 0); - } + if (pused[pmap.Get(refpi)] == 0) + pmap.Elem(refpi).Invalidate(); } } @@ -1115,15 +1114,32 @@ int Meshing3 :: ApplyRules if (found) { + /* for (i = 1; i <= tempnewpoints.Size(); i++) lpoints.Append (tempnewpoints.Get(i)); + */ + for (Point3d p : tempnewpoints) + lpoints.Append(p); + /* for (i = 1; i <= tempnewfaces.Size(); i++) if (tempnewfaces.Get(i).PNum(1)) lfaces.Append (tempnewfaces.Get(i)); + */ + for (int i : tempnewfaces.Range()) + if (tempnewfaces[i].PNum(1).IsValid()) + lfaces.Append (tempnewfaces[i]); + /* for (i = 1; i <= tempdelfaces.Size(); i++) delfaces.Append (tempdelfaces.Get(i)); + */ + for (int i : tempdelfaces.Range()) + delfaces.Append (tempdelfaces[i]); + /* for (i = 1; i <= tempelements.Size(); i++) elements.Append (tempelements.Get(i)); + */ + for (int i : tempelements.Range()) + elements.Append (tempelements[i]); } retminerr = minerr; diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp index 8cb318ea..6efb8f9b 100644 --- a/libsrc/meshing/smoothing3.cpp +++ b/libsrc/meshing/smoothing3.cpp @@ -1422,19 +1422,18 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) if(lochfunc) { - for(int i=1; i<=points.Size(); i++) - pointh[i] = GetH(points.Get(i)); + for (PointIndex pi : points.Range()) + pointh[pi] = GetH(points[pi]); } else { pointh = 0; - for(int i=0; i pointh[el.PNum(j)]) - pointh[el.PNum(j)] = h; + for (PointIndex pi : el.PNums()) + if (h > pointh[pi]) + pointh[pi] = h; } } @@ -1455,21 +1454,15 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) const char * savetask = multithread.task; multithread.task = "Smooth Mesh"; - - for (PointIndex pi = points.Begin(); pi < points.End(); pi++) + + for (PointIndex pi : points.Range()) if ( (*this)[pi].Type() == INNERPOINT && perrs[pi] > 0.01 * badmax) { if (multithread.terminate) throw NgException ("Meshing stopped"); multithread.percent = 100.0 * (pi+1-PointIndex::BASE) / points.Size(); - /* - if (points.Size() < 1000) - PrintDot (); - else - if ( (i+1-PointIndex::BASE) % 10 == 0) - PrintDot ('+'); - */ + if ( (pi+1-PointIndex::BASE) % printmod == 0) PrintDot (printdot); double lh = pointh[pi]; @@ -1503,7 +1496,6 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) } PrintDot ('\n'); - delete pf; multithread.task = savetask; diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp index 01346f40..338caedf 100644 --- a/libsrc/meshing/topology.cpp +++ b/libsrc/meshing/topology.cpp @@ -144,7 +144,7 @@ namespace netgen { swap (face.I2(), face.I3()); facedir += 2; } if (face.I1() > face.I2()) { swap (face.I1(), face.I2()); facedir += 4; } - + if (face.I1() != v) continue; func (face, elnr, j, true, facedir); @@ -608,15 +608,15 @@ namespace netgen { case 3: edges[elnr][loc_edge].nr = edgenum; - edges[elnr][loc_edge].orient = edgedir; + // edges[elnr][loc_edge].orient = edgedir; break; case 2: surfedges[elnr][loc_edge].nr = edgenum; - surfedges[elnr][loc_edge].orient = edgedir; + // surfedges[elnr][loc_edge].orient = edgedir; break; case 1: segedges[elnr].nr = edgenum; - segedges[elnr].orient = edgedir; + // segedges[elnr].orient = edgedir; break; } }); @@ -783,12 +783,12 @@ namespace netgen if (volume) { faces[elnr][j].fnr = facenum; - faces[elnr][j].forient = facedir; + // faces[elnr][j].forient = facedir; } else { surffaces[elnr].fnr = facenum; - surffaces[elnr].forient = facedir; + // surffaces[elnr].forient = facedir; } }); } @@ -1346,7 +1346,8 @@ namespace netgen eorient.SetSize (ned); for (int i = 1; i <= ned; i++) // eorient.Elem(i) = (edges.Get(elnr)[i-1] > 0) ? 1 : -1; - eorient.Elem(i) = (edges.Get(elnr)[i-1].orient) ? -1 : 1; + // eorient.Elem(i) = (edges.Get(elnr)[i-1].orient) ? -1 : 1; + eorient.Elem(i) = GetElementEdgeOrientation (elnr, i-1) ? -1 : 1; } void MeshTopology :: GetElementFaceOrientations (int elnr, Array & forient) const @@ -1354,8 +1355,9 @@ namespace netgen int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType()); forient.SetSize (nfa); for (int i = 1; i <= nfa; i++) - forient.Elem(i) = faces.Get(elnr)[i-1].forient; - // forient.Elem(i) = (faces.Get(elnr)[i-1]-1) % 8; + // forient.Elem(i) = faces.Get(elnr)[i-1].forient; + // forient.Elem(i) = (faces.Get(elnr)[i-1]-1) % 8; + forient.Elem(i) = GetElementFaceOrientation(elnr, i-1); } @@ -1377,7 +1379,8 @@ namespace netgen */ if (edges.Get(elnr)[i].nr == -1) return i; eledges[i] = edges.Get(elnr)[i].nr+1; - orient[i] = edges.Get(elnr)[i].orient ? -1 : 1; + // orient[i] = edges.Get(elnr)[i].orient ? -1 : 1; + orient[i] = GetElementEdgeOrientation(elnr, i) ? -1 : 1; } } else @@ -1432,7 +1435,8 @@ namespace netgen */ if (faces.Get(elnr)[i].fnr == -1) return i; elfaces[i] = faces.Get(elnr)[i].fnr+1; - orient[i] = faces.Get(elnr)[i].forient; + // orient[i] = faces.Get(elnr)[i].forient; + orient[i] = GetElementFaceOrientation (elnr, i); } } else @@ -1484,13 +1488,15 @@ namespace netgen eorient.SetSize (ned); for (int i = 0; i < ned; i++) // eorient[i] = (surfedges.Get(elnr)[i] > 0) ? 1 : -1; - eorient[i] = (surfedges.Get(elnr)[i].orient) ? -1 : 1; + // eorient[i] = (surfedges.Get(elnr)[i].orient) ? -1 : 1; + eorient[i] = GetSurfaceElementEdgeOrientation(elnr, i) ? -1 : 1; } int MeshTopology :: GetSurfaceElementFaceOrientation (int elnr) const { // return (surffaces.Get(elnr)-1) % 8; - return surffaces.Get(elnr).forient; + // return surffaces.Get(elnr).forient; + return GetSurfaceElementFaceOrientation2(elnr); } int MeshTopology :: GetSurfaceElementEdges (int elnr, int * eledges, int * orient) const @@ -1509,7 +1515,8 @@ namespace netgen */ if (surfedges.Get(elnr)[i].nr == -1) return i; eledges[i] = surfedges.Get(elnr)[i].nr+1; - orient[i] = (surfedges.Get(elnr)[i].orient) ? -1 : 1; + // orient[i] = (surfedges.Get(elnr)[i].orient) ? -1 : 1; + orient[i] = GetSurfaceElementEdgeOrientation(elnr, i) ? -1 : 1; } } @@ -1536,12 +1543,159 @@ namespace netgen */ eledges[0] = segedges.Get(elnr).nr+1; if (orient) - orient[0] = segedges.Get(elnr).orient ? -1 : 1; + // orient[0] = segedges.Get(elnr).orient ? -1 : 1; + orient[0] = GetSegmentEdgeOrientation(elnr) ? -1 : 1; } return 1; } + int MeshTopology :: GetElementEdgeOrientation (int elnr, int locedgenr) const + { + const Element & el = mesh.VolumeElement (elnr); + const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType()); + + int k = locedgenr; + INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); + int edgedir = (edge.I1() > edge.I2()); + return edgedir; + } + + + int MeshTopology :: GetElementFaceOrientation (int elnr, int locfacenr) const + { + const Element & el = mesh.VolumeElement (elnr); + + const ELEMENT_FACE * elfaces = MeshTopology::GetFaces0 (el.GetType()); + + int j = locfacenr; + if (elfaces[j][3] < 0) + { // triangle + INDEX_4 face(el[elfaces[j][0]], el[elfaces[j][1]], + el[elfaces[j][2]], 0); + + int facedir = 0; + if (face.I1() > face.I2()) + { swap (face.I1(), face.I2()); facedir += 1; } + if (face.I2() > face.I3()) + { swap (face.I2(), face.I3()); facedir += 2; } + if (face.I1() > face.I2()) + { swap (face.I1(), face.I2()); facedir += 4; } + + return facedir; + } + else + { + // quad + int facenum; + INDEX_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()) > + min2 (face4.I4(), face4.I3())) + { // z - flip + facedir += 1; + swap (face4.I1(), face4.I4()); + swap (face4.I2(), face4.I3()); + } + if (min2 (face4.I1(), face4.I4()) > + min2 (face4.I2(), face4.I3())) + { // x - flip + facedir += 2; + swap (face4.I1(), face4.I2()); + swap (face4.I3(), face4.I4()); + } + if (face4.I2() > face4.I4()) + { // diagonal flip + facedir += 4; + swap (face4.I2(), face4.I4()); + } + + return facedir; + } + } + + + + int MeshTopology :: GetSurfaceElementEdgeOrientation (int elnr, int locedgenr) const + { + const Element2d & el = mesh.SurfaceElement (elnr); + const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType()); + + int k = locedgenr; + INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); + int edgedir = (edge.I1() > edge.I2()); + return edgedir; + } + + int MeshTopology :: GetSurfaceElementFaceOrientation2 (int elnr) const + { + const Element2d & el = mesh.SurfaceElement (elnr); + + const ELEMENT_FACE * elfaces = MeshTopology::GetFaces0 (el.GetType()); + + int j = 0; + if (elfaces[j][3] < 0) + { // triangle + INDEX_4 face(el[elfaces[j][0]], el[elfaces[j][1]], + el[elfaces[j][2]], 0); + + int facedir = 0; + if (face.I1() > face.I2()) + { swap (face.I1(), face.I2()); facedir += 1; } + if (face.I2() > face.I3()) + { swap (face.I2(), face.I3()); facedir += 2; } + if (face.I1() > face.I2()) + { swap (face.I1(), face.I2()); facedir += 4; } + + return facedir; + } + else + { + // quad + int facenum; + INDEX_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()) > + min2 (face4.I4(), face4.I3())) + { // z - flip + facedir += 1; + swap (face4.I1(), face4.I4()); + swap (face4.I2(), face4.I3()); + } + if (min2 (face4.I1(), face4.I4()) > + min2 (face4.I2(), face4.I3())) + { // x - flip + facedir += 2; + swap (face4.I1(), face4.I2()); + swap (face4.I3(), face4.I4()); + } + if (face4.I2() > face4.I4()) + { // diagonal flip + facedir += 4; + swap (face4.I2(), face4.I4()); + } + + return facedir; + } + } + + int MeshTopology :: GetSegmentEdgeOrientation (int elnr) const + { + const Segment & el = mesh.LineSegment (elnr); + const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType()); + + int k = 0; + INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); + int edgedir = (edge.I1() > edge.I2()); + return edgedir; + } + + + void MeshTopology :: GetFaceVertices (int fnr, Array & vertices) const { vertices.SetSize(4); diff --git a/libsrc/meshing/topology.hpp b/libsrc/meshing/topology.hpp index b41e0e92..21da2b23 100644 --- a/libsrc/meshing/topology.hpp +++ b/libsrc/meshing/topology.hpp @@ -17,14 +17,14 @@ namespace netgen struct T_EDGE { - int orient:1; - int nr:31; // 0-based + // int orient:1; + int nr; // 0-based }; struct T_FACE { - int forient:3; - int fnr:29; // 0-based + // int forient:3; + int fnr; // 0-based }; @@ -87,7 +87,8 @@ public: void GetSegmentEdge (int segnr, int & enr, int & orient) const { enr = segedges.Get(segnr).nr+1; - orient = segedges.Get(segnr).orient; + // orient = segedges.Get(segnr).orient; + orient = GetSegmentEdgeOrientation(segnr); } void GetElementEdges (int elnr, Array & edges) const; @@ -98,6 +99,13 @@ public: int GetElementEdges (int elnr, int * edges, int * orient) const; int GetElementFaces (int elnr, int * faces, int * orient) const; + int GetElementEdgeOrientation (int elnr, int locedgenr) const; // old style + int GetElementFaceOrientation (int elnr, int locfacenr) const; // old style + int GetSurfaceElementEdgeOrientation (int elnr, int locedgenr) const; // old style + int GetSurfaceElementFaceOrientation2 (int elnr) const; // old style + int GetSegmentEdgeOrientation (int elnr) const; // old style + + void GetFaceVertices (int fnr, Array & vertices) const; void GetFaceVertices (int fnr, int * vertices) const; void GetEdgeVertices (int enr, int & v1, int & v2) const; diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index 9e79c2a1..ca74b5b0 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -441,8 +441,7 @@ namespace netgen // External access to the mesh generation functions within the OCC // subsystem (Not sure if this is the best way to implement this....!!) extern int OCCGenerateMesh (OCCGeometry & occgeometry, shared_ptr & mesh, - MeshingParameters & mparam, - int perfstepsstart, int perfstepsend); + MeshingParameters & mparam); extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh); diff --git a/libsrc/stlgeom/stlgeom.cpp b/libsrc/stlgeom/stlgeom.cpp index 91f3ab22..a89a1757 100644 --- a/libsrc/stlgeom/stlgeom.cpp +++ b/libsrc/stlgeom/stlgeom.cpp @@ -455,7 +455,8 @@ int STLGeometry :: AddEdge(int ap1, int ap2) STLEdge e(ap1,ap2); e.SetLeftTrig(GetLeftTrig(ap1,ap2)); e.SetRightTrig(GetRightTrig(ap1,ap2)); - return edges.Append(e); + edges.Append(e); + return edges.Size(); } void STLGeometry :: STLDoctorConfirmEdge() diff --git a/libsrc/stlgeom/stlgeom.hpp b/libsrc/stlgeom/stlgeom.hpp index 59a3342e..cb51cfae 100644 --- a/libsrc/stlgeom/stlgeom.hpp +++ b/libsrc/stlgeom/stlgeom.hpp @@ -302,11 +302,11 @@ namespace netgen DLL_HEADER int GetNodeOfSelTrig() const; - int AddNormal(const Vec3d& n) {return normals.Append(n);} + int AddNormal(const Vec3d& n) { normals.Append(n); return normals.Size(); } const Vec3d & GetNormal(int nr) const {return normals.Get(nr);} void SetNormal(int nr, const Vec3d& n) {normals.Elem(nr) = n;} - int AddEdge(const STLEdge& v) {return edges.Append(v);} + int AddEdge(const STLEdge& v) { edges.Append(v); return edges.Size(); } int AddEdge(int p1, int p2); STLEdge GetEdge(int nr) {return edges.Get(nr);} @@ -434,7 +434,7 @@ namespace netgen int ProjectOnWholeSurface (Point<3> & p3d) const; int GetNLines() const {return lines.Size();} - int AddLine(STLLine* line) {return lines.Append(line);} + int AddLine(STLLine* line) { lines.Append(line); return lines.Size(); } STLLine* GetLine(int nr) const {return lines.Get(nr);} int GetLineP(int lnr, int pnr) const {return lines.Get(lnr)->PNum(pnr);} int GetLineNP(int nr) const {return lines.Get(nr)->NP();} diff --git a/libsrc/stlgeom/stltool.cpp b/libsrc/stlgeom/stltool.cpp index 855cec42..ee42da1b 100644 --- a/libsrc/stlgeom/stltool.cpp +++ b/libsrc/stlgeom/stltool.cpp @@ -19,7 +19,8 @@ int AddPointIfNotExists(Array& ap, const Point3d& p, double eps) for (int i = 1; i <= ap.Size(); i++) if (Dist2(ap.Get(i),p) <= eps2 ) return i; - return ap.Append(p); + ap.Append(p); + return ap.Size(); } //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ diff --git a/libsrc/stlgeom/stltopology.cpp b/libsrc/stlgeom/stltopology.cpp index f4078ead..53d3d9f0 100644 --- a/libsrc/stlgeom/stltopology.cpp +++ b/libsrc/stlgeom/stltopology.cpp @@ -592,7 +592,8 @@ void STLTopology :: FindNeighbourTrigs() } else { - enr = topedges.Append (STLTopEdge (pi1, pi2, i, 0)); + topedges.Append (STLTopEdge (pi1, pi2, i, 0)); + enr = topedges.Size(); ht_topedges->Set (i2, enr); trig.EdgeNum(j) = enr; } diff --git a/libsrc/stlgeom/stltopology.hpp b/libsrc/stlgeom/stltopology.hpp index 37c13603..7351615d 100644 --- a/libsrc/stlgeom/stltopology.hpp +++ b/libsrc/stlgeom/stltopology.hpp @@ -293,7 +293,7 @@ public: int GetNP() const { return points.Size(); } - int AddPoint(const Point<3> & p) { return points.Append(p); } + int AddPoint(const Point<3> & p) { points.Append(p); return points.Size(); } const Point<3> & GetPoint(int nr) const { return points.Get(nr); } int GetPointNum (const Point<3> & p); void SetPoint(int nr, const Point<3> & p) { points.Elem(nr) = p; } diff --git a/libsrc/visualization/mvdraw.cpp b/libsrc/visualization/mvdraw.cpp index 04d98e6c..f5d00922 100644 --- a/libsrc/visualization/mvdraw.cpp +++ b/libsrc/visualization/mvdraw.cpp @@ -63,8 +63,8 @@ namespace netgen int VisualScene :: selface; int VisualScene :: selelement; - int VisualScene :: selpoint; - int VisualScene :: selpoint2; + PointIndex VisualScene :: selpoint; + PointIndex VisualScene :: selpoint2; int VisualScene :: locpi; int VisualScene :: seledge; diff --git a/libsrc/visualization/mvdraw.hpp b/libsrc/visualization/mvdraw.hpp index a1977572..39efbc5f 100644 --- a/libsrc/visualization/mvdraw.hpp +++ b/libsrc/visualization/mvdraw.hpp @@ -26,8 +26,8 @@ namespace netgen static int DLL_HEADER selface; static int selelement; - static int DLL_HEADER selpoint; - static int selpoint2; + static PointIndex DLL_HEADER selpoint; + static PointIndex selpoint2; static int locpi; static int DLL_HEADER seledge; @@ -235,8 +235,8 @@ namespace netgen const Point3d & center, const double rad, const int displaylist, - int & selelement, int & selface, int & seledge, int & selpoint, - int & selpoint2, int & locpi); + int & selelement, int & selface, int & seledge, PointIndex & selpoint, + PointIndex & selpoint2, int & locpi); } diff --git a/libsrc/visualization/vispar.hpp b/libsrc/visualization/vispar.hpp index 1661ac3b..94026440 100644 --- a/libsrc/visualization/vispar.hpp +++ b/libsrc/visualization/vispar.hpp @@ -76,7 +76,7 @@ public: int drawcurveprojedge; - int centerpoint; + PointIndex centerpoint; int drawelement; // stl: diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp index 2a58e6e6..08c2f6de 100644 --- a/libsrc/visualization/vsmesh.cpp +++ b/libsrc/visualization/vsmesh.cpp @@ -440,12 +440,12 @@ namespace netgen char buf[30]; if (vispar.drawpointnumbers) - for (int i = 1; i <= mesh->GetNP(); i++) + for (PointIndex pi : mesh->Points().Range()) { - const Point3d & p = mesh->Point(i); + const Point3d & p = mesh->Point(pi); glRasterPos3d (p.X(), p.Y(), p.Z()); - sprintf (buf, "%d", i); + sprintf (buf, "%d", int(pi)); // glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf); MyOpenGLText (buf); @@ -663,12 +663,12 @@ namespace netgen - for (int i = 1; i <= mesh->GetNE(); i++) + for (ElementIndex ei : mesh->VolumeElements().Range()) { - if (mesh->VolumeElement(i).flags.badel) + if (mesh->VolumeElement(ei).flags.badel) { // copy to be thread-safe - Element el = mesh->VolumeElement (i); + Element el = mesh->VolumeElement (ei); if ( (el.GetNP() == 4) || (el.GetNP() == 10)) { glBegin (GL_LINES); @@ -747,17 +747,16 @@ namespace netgen } } - - for (int i = 1; i <= mesh->GetNSE(); i++) + for (SurfaceElementIndex sei : mesh->SurfaceElements().Range()) { - Element2d el = mesh->SurfaceElement(i); + Element2d el = mesh->SurfaceElement(sei); // copy to be thread-safe if (!el.BadElement()) continue; - int drawel = 1; + bool drawel = true; for (int j = 1; j <= el.GetNP(); j++) - if (!el.PNum(j)) - drawel = 0; + if (!el.PNum(j).IsValid()) + drawel = false; if (!drawel) continue; @@ -3343,8 +3342,8 @@ namespace netgen const Point3d & center, const double rad, const int displaylist, - int & selelement, int & selface, int & seledge, int & selpoint, - int & selpoint2, int & locpi) + int & selelement, int & selface, int & seledge, PointIndex & selpoint, + PointIndex & selpoint2, int & locpi) { auto mesh = vsmesh.GetMesh(); diff --git a/libsrc/visualization/vssolution.cpp b/libsrc/visualization/vssolution.cpp index a83fd9c0..18b217a4 100644 --- a/libsrc/visualization/vssolution.cpp +++ b/libsrc/visualization/vssolution.cpp @@ -4216,8 +4216,8 @@ namespace netgen Point<3> p2 = locgrid[teti[pi2]]; cppt.lami = p2 + edgelam[ednr] * (p1-p2); - - pnr = pts.Append (cppt)-1; + pts.Append (cppt); + pnr = pts.Size()-1; edges.Set (pair, pnr); }