From d215ac1025641472cc2c21b810e4c551e32470a6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Joachim=20Sch=C3=B6berl?= Date: Sat, 10 Aug 2019 00:21:37 +0200 Subject: [PATCH] T_POINTS are now ngcore::Array --- libsrc/core/array.hpp | 84 +++++++++++++++---------- libsrc/meshing/findip.hpp | 12 ++-- libsrc/meshing/improve2.cpp | 28 +++++---- libsrc/meshing/meshclass.cpp | 39 +++++++++--- libsrc/meshing/meshtype.cpp | 22 +++---- libsrc/meshing/meshtype.hpp | 32 ++++------ libsrc/meshing/parallelmesh.cpp | 8 +-- libsrc/meshing/smoothing3.cpp | 108 ++++++++++++++++---------------- 8 files changed, 186 insertions(+), 147 deletions(-) diff --git a/libsrc/core/array.hpp b/libsrc/core/array.hpp index b47ae182..06a8dc40 100644 --- a/libsrc/core/array.hpp +++ b/libsrc/core/array.hpp @@ -59,6 +59,28 @@ namespace ngcore } + + template + class AOWrapperIterator + { + const AO & ao; + size_t ind; + public: + NETGEN_INLINE AOWrapperIterator (const AO & aao, size_t ai) + : ao(aao), ind(ai) { ; } + NETGEN_INLINE AOWrapperIterator operator++ (int) + { return AOWrapperIterator(ao, ind++); } + NETGEN_INLINE AOWrapperIterator operator++ () + { return AOWrapperIterator(ao, ++ind); } + NETGEN_INLINE auto operator*() const -> decltype(ao[ind]) { return ao[ind]; } + NETGEN_INLINE auto operator*() -> decltype(ao[ind]) { return ao[ind]; } + NETGEN_INLINE bool operator != (AOWrapperIterator d2) { return ind != d2.ind; } + NETGEN_INLINE bool operator == (AOWrapperIterator d2) { return ind == d2.ind; } + }; + + + + /* Some class which can be treated as array */ @@ -99,27 +121,12 @@ namespace ngcore // NETGEN_INLINE auto & operator[] (size_t i) { return Spec()[i]; } NETGEN_INLINE auto operator[] (size_t i) const { return Spec()[i]; } + // NETGEN_INLINE auto begin() const { return Spec().begin(); } + // NETGEN_INLINE auto end() const { return Spec().end(); } + NETGEN_INLINE auto begin () const { return AOWrapperIterator (*this, 0); } + NETGEN_INLINE auto end () const { return AOWrapperIterator (*this, Size()); } }; - - - template - class AOWrapperIterator - { - const AO & ao; - size_t ind; - public: - NETGEN_INLINE AOWrapperIterator (const AO & aao, size_t ai) - : ao(aao), ind(ai) { ; } - NETGEN_INLINE AOWrapperIterator operator++ (int) - { return AOWrapperIterator(ao, ind++); } - NETGEN_INLINE AOWrapperIterator operator++ () - { return AOWrapperIterator(ao, ++ind); } - NETGEN_INLINE auto operator*() const -> decltype(ao[ind]) { return ao[ind]; } - NETGEN_INLINE auto operator*() -> decltype(ao[ind]) { return ao[ind]; } - NETGEN_INLINE bool operator != (AOWrapperIterator d2) { return ind != d2.ind; } - NETGEN_INLINE bool operator == (AOWrapperIterator d2) { return ind == d2.ind; } - }; - + template @@ -262,7 +269,8 @@ namespace ngcore NETGEN_INLINE auto Size() const { return next-first; } NETGEN_INLINE T operator[] (T i) const { return first+i; } NETGEN_INLINE bool Contains (T i) const { return ((i >= first) && (i < next)); } - + NETGEN_INLINE T_Range Modify(int inc_beg, int inc_end) const + { return T_Range(first+inc_beg, next+inc_end); } NETGEN_INLINE ArrayRangeIterator begin() const { return first; } NETGEN_INLINE ArrayRangeIterator end() const { return next; } @@ -437,7 +445,7 @@ namespace ngcore { size_t hsize = size; T * hdata = data; - for (size_t i = 0; i < hsize; i++) hdata[i] = a2[i]; + for (size_t i = 0; i < hsize; i++) hdata[i] = a2.data[i]; return *this; } @@ -446,14 +454,15 @@ namespace ngcore { size_t hsize = size; T * hdata = data; - for (size_t i = 0; i < hsize; i++) hdata[i] = a2[i]; + auto p2 = a2.begin(); + for (size_t i = 0; i < hsize; i++, p2++) hdata[i] = *p2; return *this; } NETGEN_INLINE const FlatArray & operator= (const std::function & func) const { for (size_t i = 0; i < size; i++) - data[i] = func(i); + data[i] = func(i+BASE); return *this; } @@ -481,18 +490,18 @@ namespace ngcore /// Access array. range check by macro NETGEN_CHECK_RANGE NETGEN_INLINE T & operator[] (IndexType i) const { - NETGEN_CHECK_RANGE(i,0,size); + NETGEN_CHECK_RANGE(i,BASE,size+BASE); return data[i-BASE]; } NETGEN_INLINE T_Range Range () const { - return T_Range (0, Size()); + return T_Range (BASE, Size()+BASE); } NETGEN_INLINE const CArray Addr (size_t pos) const { - return CArray (data+pos); + return CArray (data+pos-BASE); } // const CArray operator+ (int pos) @@ -606,6 +615,7 @@ namespace ngcore using FlatArray::size; using FlatArray::data; + using FlatArray::BASE; public: /// Generate array of logical and physical size asize @@ -674,8 +684,14 @@ namespace ngcore { allocsize = size; mem_to_delete = data; + /* for (size_t i = 0; i < size; i++) data[i] = a2[i]; + */ + auto p2 = a2.begin(); + for (size_t i = 0; i < size; i++, p2++) + data[i] = *p2; + } Array (std::initializer_list list) @@ -697,9 +713,9 @@ namespace ngcore allocsize = size; mem_to_delete = data; for(size_t i = 0; i < a2.Size(); i++) - (*this)[i] = a2[i]; + data[i] = a2[i]; for (size_t i = a2.Size(), j=0; i < size; i++,j++) - (*this)[i] = a3[j]; + data[i] = a3[j]; } /// if responsible, deletes memory @@ -842,8 +858,8 @@ namespace ngcore /// Delete element i. Move last element to position i. NETGEN_INLINE void DeleteElement (size_t i) { - NETGEN_CHECK_RANGE(i,0,size); - data[i] = std::move(data[size-1]); + NETGEN_CHECK_RANGE(i,BASE,BASE+size); + data[i-BASE] = std::move(data[size-1]); size--; } @@ -876,7 +892,7 @@ namespace ngcore /// Fill array with val NETGEN_INLINE Array & operator= (const T & val) { - FlatArray::operator= (val); + FlatArray::operator= (val); return *this; } @@ -886,7 +902,7 @@ namespace ngcore SetSize0 (); SetSize (a2.Size()); for (size_t i = 0; i < size; i++) - (*this)[i] = a2[i]; + data[i] = a2.data[i]; return *this; } @@ -906,7 +922,7 @@ namespace ngcore { SetSize (a2.Size()); for (size_t i = 0; i < size; i++) - (*this)[i] = a2[i]; + data[i] = a2[i]; return *this; } diff --git a/libsrc/meshing/findip.hpp b/libsrc/meshing/findip.hpp index 10caa73f..dc3c6744 100644 --- a/libsrc/meshing/findip.hpp +++ b/libsrc/meshing/findip.hpp @@ -90,9 +90,9 @@ inline int FindInnerPoint (POINTArray & points, for (int i = 0; i < nf; i++) { - Point3d p1 = points.Get(faces[i][0]); - a[i] = Cross (points.Get(faces[i][1]) - p1, - points.Get(faces[i][2]) - p1); + Point3d p1 = points[faces[i][0]]; + a[i] = Cross (points[faces[i][1]] - p1, + points[faces[i][2]] - p1); a[i] /= a[i].Length(); c[i] = - (a[i].X() * p1.X() + a[i].Y() * p1.Y() + a[i].Z() * p1.Z()); } @@ -107,7 +107,7 @@ inline int FindInnerPoint (POINTArray & points, center = 0; for (int i = 0; i < faces.Size(); i++) for (int j = 0; j < 3; j++) - center += Vec<3> (points.Get(faces[i][j])); + center += Vec<3> (points[faces[i][j]]); center /= (3*faces.Size()); @@ -120,8 +120,8 @@ inline int FindInnerPoint (POINTArray & points, // (*testout) << "el[" << i << "] = " << el << endl; for (int j = 1; j <= 3; j++) { - double hi = Dist (points.Get(faces[i].PNumMod(j)), - points.Get(faces[i].PNumMod(j+1))); + double hi = Dist (points[faces[i].PNumMod(j)], + points[faces[i].PNumMod(j+1)]); if (hi > hmax) hmax = hi; } } diff --git a/libsrc/meshing/improve2.cpp b/libsrc/meshing/improve2.cpp index 7afc65dc..332ae955 100644 --- a/libsrc/meshing/improve2.cpp +++ b/libsrc/meshing/improve2.cpp @@ -332,7 +332,7 @@ namespace netgen (nvp4 * nv4 > critval); - double horder = Dist (mesh.Point(pi1), mesh.Point(pi2)); + double horder = Dist (mesh[pi1], mesh[pi2]); if ( // nv1 * nv2 >= 0 && nv1.Length() > 1e-3 * horder * horder && @@ -343,8 +343,8 @@ namespace netgen { int e = pdef[pi1] + pdef[pi2] - pdef[pi3] - pdef[pi4]; double d = - Dist2 (mesh.Point(pi1), mesh.Point(pi2)) - - Dist2 (mesh.Point(pi3), mesh.Point(pi4)); + Dist2 (mesh[pi1], mesh[pi2]) - + Dist2 (mesh[pi3], mesh[pi4]); should = e >= t && (e > 2 || d > 0); } @@ -468,14 +468,16 @@ namespace netgen TABLE elementsonnode(np); NgArray hasonepi, hasbothpi; - for (int i = 0; i < seia.Size(); i++) + for (SurfaceElementIndex sei : seia) { - Element2d & el = mesh[seia[i]]; - for (int j = 0; j < el.GetNP(); j++) - elementsonnode.Add (el[j], seia[i]); + // Element2d & el = mesh[sei]; + // for (int j = 0; j < el.GetNP(); j++) + // elementsonnode.Add (el[j], sei); + for (PointIndex pi : mesh[sei].PNums()) + elementsonnode.Add (pi, sei); } - NgArray fixed(np); + Array fixed(np); fixed = false; NgProfiler::StopTimer (timerstart1); @@ -489,9 +491,9 @@ namespace netgen } */ - for (int i = 0; i < seia.Size(); i++) + for (SurfaceElementIndex sei : seia) { - Element2d & sel = mesh[seia[i]]; + Element2d & sel = mesh[sei]; for (int j = 0; j < sel.GetNP(); j++) { PointIndex pi1 = sel.PNumMod(j+2); @@ -505,10 +507,12 @@ namespace netgen } - + /* for(int i = 0; i < mesh.LockedPoints().Size(); i++) fixed[mesh.LockedPoints()[i]] = true; - + */ + for (PointIndex pi : mesh.LockedPoints()) + fixed[pi] = true; NgArray,PointIndex::BASE> normals(np); diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 0333ccc3..5ae25122 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -1824,7 +1824,7 @@ namespace netgen } } - + // BitArray base is PointIndex::BASE ... void Mesh :: FixPoints (const BitArray & fixpoints) { if (fixpoints.Size() != GetNP()) @@ -1833,11 +1833,16 @@ namespace netgen return; } int np = GetNP(); + /* for (int i = 1; i <= np; i++) if (fixpoints.Test(i)) { points.Elem(i).SetType (FIXEDPOINT); } + */ + for (PointIndex pi : points.Range()) + if (fixpoints.Test(pi)) + points[pi].SetType(FIXEDPOINT); } @@ -2410,9 +2415,13 @@ namespace netgen ptyps.Elem(seg[1]) = EDGEPOINT; } */ + /* for (int i = 1; i <= points.Size(); i++) points.Elem(i).SetType(SURFACEPOINT); - + */ + for (auto & p : points) + p.SetType (SURFACEPOINT); + for (int i = 1; i <= GetNSeg(); i++) { const Segment & seg = LineSegment (i); @@ -3229,8 +3238,8 @@ namespace netgen double Mesh :: ElementError (int eli, const MeshingParameters & mp) const { const Element & el = volelements[eli-1]; - return CalcTetBadness (points.Get(el[0]), points.Get(el[1]), - points.Get(el[2]), points.Get(el[3]), -1, mp); + return CalcTetBadness (points[el[0]], points[el[1]], + points[el[2]], points[el[3]], -1, mp); } void Mesh :: AddLockedPoint (PointIndex pi) @@ -4341,14 +4350,21 @@ namespace netgen { Box<3> box (Box<3>::EMPTY_BOX); for (SurfaceElementIndex sei = 0; sei < ne; sei++) - box.Add (points[surfelements[sei].PNums()]); + // box.Add (points[surfelements[sei].PNums()]); + for (auto pi : surfelements[sei].PNums()) + box.Add (points[pi]); box.Increase (1.01 * box.Diam()); elementsearchtree = make_unique> (box); for (SurfaceElementIndex sei = 0; sei < ne; sei++) { - box.Set (points[surfelements[sei].PNums()]); + // box.Set (points[surfelements[sei].PNums()]); + + Box<3> box (Box<3>::EMPTY_BOX); + for (auto pi : surfelements[sei].PNums()) + box.Add (points[pi]); + elementsearchtree -> Insert (box, sei+1); } } @@ -4356,14 +4372,21 @@ namespace netgen { Box<3> box (Box<3>::EMPTY_BOX); for (ElementIndex ei = 0; ei < ne; ei++) - box.Add (points[volelements[ei].PNums()]); + // box.Add (points[volelements[ei].PNums()]); + for (auto pi : volelements[ei].PNums()) + box.Add (points[pi]); box.Increase (1.01 * box.Diam()); elementsearchtree = make_unique> (box); for (ElementIndex ei = 0; ei < ne; ei++) { - box.Set (points[volelements[ei].PNums()]); + // box.Set (points[volelements[ei].PNums()]); + + Box<3> box (Box<3>::EMPTY_BOX); + for (auto pi : volelements[ei].PNums()) + box.Add (points[pi]); + elementsearchtree -> Insert (box, ei+1); } } diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index 95127c7c..b32068a0 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -295,9 +295,9 @@ namespace netgen void Element2d :: GetBox (const T_POINTS & points, Box3d & box) const { - box.SetPoint (points.Get(pnum[0])); + box.SetPoint (points[pnum[0]]); for (unsigned i = 1; i < np; i++) - box.AddPoint (points.Get(pnum[i])); + box.AddPoint (points[pnum[i]]); } bool Element2d :: operator==(const Element2d & el2) const @@ -879,7 +879,7 @@ namespace netgen for (i = 1; i <= GetNP(); i++) { - Point3d p = points.Get(PNum(i)); + Point3d p = points[PNum(i)]; pmat.Elem(1, i) = p.X() * t1(0) + p.Y() * t1(1) + p.Z() * t1(2); pmat.Elem(2, i) = p.X() * t2(0) + p.Y() * t2(1) + p.Z() * t2(2); } @@ -1172,17 +1172,17 @@ namespace netgen void Element :: GetBox (const T_POINTS & points, Box3d & box) const { - box.SetPoint (points.Get(PNum(1))); - box.AddPoint (points.Get(PNum(2))); - box.AddPoint (points.Get(PNum(3))); - box.AddPoint (points.Get(PNum(4))); + box.SetPoint (points[PNum(1)]); + box.AddPoint (points[PNum(2)]); + box.AddPoint (points[PNum(3)]); + box.AddPoint (points[PNum(4)]); } double Element :: Volume (const T_POINTS & points) const { - Vec<3> v1 = points.Get(PNum(2)) - points.Get(PNum(1)); - Vec<3> v2 = points.Get(PNum(3)) - points.Get(PNum(1)); - Vec<3> v3 = points.Get(PNum(4)) - points.Get(PNum(1)); + Vec<3> v1 = points[PNum(2)] - points[PNum(1)]; + Vec<3> v2 = points[PNum(3)] - points[PNum(1)]; + Vec<3> v3 = points[PNum(4)] - points[PNum(1)]; return -(Cross (v1, v2) * v3) / 6; } @@ -2179,7 +2179,7 @@ namespace netgen int np = GetNP(); for (int i = 1; i <= np; i++) { - const Point3d & p = points.Get(PNum(i)); + const Point3d & p = points[PNum(i)]; pmat.Elem(1, i) = p.X(); pmat.Elem(2, i) = p.Y(); pmat.Elem(3, i) = p.Z(); diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 75b2d34c..2476a421 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -369,8 +369,8 @@ namespace netgen - typedef NgArray T_POINTS; - // typedef Array T_POINTS; + // typedef NgArray T_POINTS; + typedef Array T_POINTS; @@ -500,18 +500,12 @@ namespace netgen /// const PointIndex & operator[] (int i) const { return pnum[i]; } - NgFlatArray PNums () const - { return NgFlatArray (np, &pnum[0]); } - NgFlatArray PNums () - { return NgFlatArray (np, &pnum[0]); } - auto Vertices() const - { return NgFlatArray (GetNV(), &pnum[0]); } + auto PNums () const { return FlatArray (np, &pnum[0]); } + auto PNums () { return FlatArray (np, &pnum[0]); } + auto Vertices() const { return FlatArray (GetNV(), &pnum[0]); } - NgFlatArray GeomInfo() const - { return NgFlatArray (np, &geominfo[0]); } - NgFlatArray GeomInfo() - { return NgFlatArray (np, &geominfo[0]); } - + auto GeomInfo() const { return FlatArray (np, &geominfo[0]); } + auto GeomInfo() { return FlatArray (np, &geominfo[0]); } /// PointIndex & PNum (int i) { return pnum[i-1]; } @@ -629,17 +623,17 @@ namespace netgen // Philippose - 08 August 2010 // Access functions for the new property: visible - void Visible(bool vis = 1) + void Visible(bool vis = true) { visible = vis; } bool IsVisible () const { return visible; } - void SetRefinementFlag (bool rflag = 1) + void SetRefinementFlag (bool rflag = true) { refflag = rflag; } bool TestRefinementFlag () const { return refflag; } - void SetStrongRefinementFlag (bool rflag = 1) + void SetStrongRefinementFlag (bool rflag = true) { strongrefflag = rflag; } bool TestStrongRefinementFlag () const { return strongrefflag; } @@ -798,10 +792,10 @@ namespace netgen /// const PointIndex & operator[] (int i) const { return pnum[i]; } - NgFlatArray PNums () const - { return NgFlatArray (np, &pnum[0]); } + auto PNums () const + { return FlatArray (np, &pnum[0]); } - NgFlatArray Vertices() const { return { GetNV(), &pnum[0] }; } + FlatArray Vertices() const { return { GetNV(), &pnum[0] }; } /// PointIndex & PNum (int i) { return pnum[i-1]; } diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp index 8867dbc9..c4dd2d96 100644 --- a/libsrc/meshing/parallelmesh.cpp +++ b/libsrc/meshing/parallelmesh.cpp @@ -314,7 +314,7 @@ namespace netgen MPI_Type_commit (&newtype); MPI_Request request; - MPI_Isend( &points[0], 1, newtype, dest, MPI_TAG_MESH+1, comm, &request); + MPI_Isend( &points[PointIndex(0)], 1, newtype, dest, MPI_TAG_MESH+1, comm, &request); sendrequests.Append (request); } @@ -474,7 +474,7 @@ namespace netgen { if(ided_sel[sei]!=-1) continue; const Element2d & sel = (*this)[sei]; - NgFlatArray points = sel.PNums(); + auto points = sel.PNums(); auto ided1 = per_verts[points[0]]; os1.SetSize(0); for (int j = 0; j < ided1.Size(); j++) @@ -498,7 +498,7 @@ namespace netgen throw NgException("SurfaceElement identified with more than one other??"); } const Element2d & sel2 = (*this)[sei]; - NgFlatArray points2 = sel2.PNums(); + auto points2 = sel2.PNums(); has_ided_sels = true; ided_sel[sei] = os1[0]; ided_sel[os1[0]] = sei; @@ -842,7 +842,7 @@ namespace netgen MPI_Datatype mptype = MeshPoint::MyGetMPIType(); MPI_Status status; - MPI_Recv( &points[1], numvert, mptype, 0, MPI_TAG_MESH+1, comm, &status); + MPI_Recv( &points[PointIndex(PointIndex::BASE)], numvert, mptype, 0, MPI_TAG_MESH+1, comm, &status); NgArray pp_data; MyMPI_Recv(pp_data, 0, MPI_TAG_MESH+1, comm); diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp index b4cbfc89..cb34b10c 100644 --- a/libsrc/meshing/smoothing3.cpp +++ b/libsrc/meshing/smoothing3.cpp @@ -1023,12 +1023,12 @@ double JacobianPointFunction :: Func (const Vector & v) const int j; double badness = 0; - Point<3> hp = points.Elem(actpind); + Point<3> hp = points[actpind]; - points.Elem(actpind) = hp + Vec<3> (v(0), v(1), v(2)); + points[actpind] = hp + Vec<3> (v(0), v(1), v(2)); if(onplane) - points.Elem(actpind) -= (v(0)*nv(0)+v(1)*nv(1)+v(2)*nv(2)) * nv; + points[actpind] -= (v(0)*nv(0)+v(1)*nv(1)+v(2)*nv(2)) * nv; for (j = 1; j <= elementsonpoint.EntrySize(actpind); j++) @@ -1037,7 +1037,7 @@ double JacobianPointFunction :: Func (const Vector & v) const badness += elements[eli-1].CalcJacobianBadness (points); } - points.Elem(actpind) = hp; + points[actpind] = hp; return badness; } @@ -1053,11 +1053,11 @@ FuncGrad (const Vector & x, Vector & g) const int lpi; double badness = 0;//, hbad; - Point<3> hp = points.Elem(actpind); - points.Elem(actpind) = hp + Vec<3> (x(0), x(1), x(2)); + Point<3> hp = points[actpind]; + points[actpind] = hp + Vec<3> (x(0), x(1), x(2)); if(onplane) - points.Elem(actpind) -= (x(0)*nv(0)+x(1)*nv(1)+x(2)*nv(2)) * nv; + points[actpind] -= (x(0)*nv(0)+x(1)*nv(1)+x(2)*nv(2)) * nv; Vec<3> hderiv; //Vec3d vdir; @@ -1108,7 +1108,7 @@ FuncGrad (const Vector & x, Vector & g) const //(*testout) << "g = " << g << endl; - points.Elem(actpind) = hp; + points[actpind] = hp; return badness; } @@ -1121,11 +1121,11 @@ FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const int lpi; double badness = 0; - Point<3> hp = points.Elem(actpind); - points.Elem(actpind) = Point<3> (hp + Vec3d (x(0), x(1), x(2))); + Point<3> hp = points[actpind]; + points[actpind] = Point<3> (hp + Vec3d (x(0), x(1), x(2))); if(onplane) - points.Elem(actpind) -= (Vec3d (x(0), x(1), x(2))*nv) * nv; + points[actpind] -= (Vec3d (x(0), x(1), x(2))*nv) * nv; double hderiv; deriv = 0; @@ -1153,7 +1153,7 @@ FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const deriv += hderiv; } - points.Elem(actpind) = hp; + points[actpind] = hp; return badness; @@ -1476,7 +1476,7 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp, OPTIMIZEGOAL goal, const BitArray * usepoint) { - int i, j; + // int i, j; (*testout) << "Improve Mesh Jacobian" << "\n"; PrintMessage (3, "ImproveMesh Jacobian"); @@ -1499,30 +1499,31 @@ void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp, BitArray badnodes(np); badnodes.Clear(); - for (i = 1; i <= ne; i++) + for (int i = 1; i <= ne; i++) { const Element & el = VolumeElement(i); double bad = el.CalcJacobianBadness (Points()); if (bad > 1) - for (j = 1; j <= el.GetNP(); j++) + for (int j = 1; j <= el.GetNP(); j++) badnodes.Set (el.PNum(j)); } - NgArray pointh (points.Size()); + NgArray pointh (points.Size()); if(lochfunc) { - for(i = 1; i<=points.Size(); i++) - pointh[i] = GetH(points.Get(i)); + // for(i = 1; i<=points.Size(); i++) + for (PointIndex pi : points.Range()) + pointh[pi] = GetH(points[pi]); } else { pointh = 0; - for(i=0; i pointh[el.PNum(j)]) pointh[el.PNum(j)] = h; } @@ -1539,12 +1540,12 @@ void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp, if ((*this)[pi].Type() != INNERPOINT) continue; - if(usepoint && !usepoint->Test(i)) + if(usepoint && !usepoint->Test(pi)) continue; //(*testout) << "improvejac, p = " << i << endl; - if (goal == OPT_WORSTCASE && !badnodes.Test(i)) + if (goal == OPT_WORSTCASE && !badnodes.Test(pi)) continue; // (*testout) << "smooth p " << i << endl; @@ -1555,15 +1556,15 @@ void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp, if (multithread.terminate) throw NgException ("Meshing stopped"); - multithread.percent = 100.0 * i / points.Size(); + multithread.percent = 100.0 * pi / points.Size(); if (points.Size() < 1000) PrintDot (); else - if (i % 10 == 0) + if (pi % 10 == 0) PrintDot ('+'); - double lh = pointh[i]; + double lh = pointh[pi]; par.typx = lh; pf.SetPointIndex (pi); @@ -1576,9 +1577,9 @@ void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp, //*testout << "start BFGS, Jacobian" << endl; BFGS (x, pf, par); //*testout << "end BFGS, Jacobian" << endl; - points.Elem(i)(0) += x(0); - points.Elem(i)(1) += x(1); - points.Elem(i)(2) += x(2); + points[pi](0) += x(0); + points[pi](1) += x(1); + points[pi](2) += x(2); } else { @@ -1601,7 +1602,7 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp, OPTIMIZEGOAL goal, const NgArray< NgArray* > * idmaps) { - int i, j; + // int i, j; (*testout) << "Improve Mesh Jacobian" << "\n"; PrintMessage (3, "ImproveMesh Jacobian"); @@ -1625,7 +1626,7 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp, { used_idmaps = &locidmaps; - for(i=1; i<=GetIdentifications().GetMaxNr(); i++) + for(int i=1; i<=GetIdentifications().GetMaxNr(); i++) { if(GetIdentifications().GetType(i) == Identifications::PERIODIC) { @@ -1655,12 +1656,12 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp, BitArray badnodes(np); badnodes.Clear(); - for (i = 1; i <= ne; i++) + for (int i = 1; i <= ne; i++) { const Element & el = VolumeElement(i); double bad = el.CalcJacobianBadness (Points()); if (bad > 1) - for (j = 1; j <= el.GetNP(); j++) + for (int j = 1; j <= el.GetNP(); j++) badnodes.Set (el.PNum(j)); } @@ -1668,17 +1669,18 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp, if(lochfunc) { - for(i=1; i<=points.Size(); i++) - pointh[i] = GetH(points.Get(i)); + // for(i=1; i<=points.Size(); i++) + for (PointIndex pi : points.Range()) + pointh[pi] = GetH(points[pi]); } else { pointh = 0; - for(i=0; i pointh[el.PNum(j)]) pointh[el.PNum(j)] = h; } @@ -1690,11 +1692,11 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp, // for (PointIndex pi = points.Begin(); pi <= points.End(); pi++) for (PointIndex pi : points.Range()) - if ( usepoint.Test(i) ) + if ( usepoint.Test(pi) ) { //(*testout) << "improvejac, p = " << i << endl; - if (goal == OPT_WORSTCASE && !badnodes.Test(i)) + if (goal == OPT_WORSTCASE && !badnodes.Test(pi)) continue; // (*testout) << "smooth p " << i << endl; @@ -1705,15 +1707,15 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp, if (multithread.terminate) throw NgException ("Meshing stopped"); - multithread.percent = 100.0 * i / points.Size(); + multithread.percent = 100.0 * pi / points.Size(); if (points.Size() < 1000) PrintDot (); else - if (i % 10 == 0) + if (pi % 10 == 0) PrintDot ('+'); - double lh = pointh[i];//GetH(points.Get(i)); + double lh = pointh[pi];//GetH(points.Get(i)); par.typx = lh; pf.SetPointIndex (pi); @@ -1721,29 +1723,29 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp, PointIndex brother (-1); if(usesum) { - for(j=0; brother == -1 && jSize(); j++) + for(int j=0; brother == -1 && jSize(); j++) { - if(i < (*used_idmaps)[j]->Size() + PointIndex::BASE) + if(pi < (*used_idmaps)[j]->Size() + PointIndex::BASE) { - brother = (*(*used_idmaps)[j])[i]; - if(brother == i || brother == 0) + brother = (*(*used_idmaps)[j])[pi]; + if(brother == pi || brother == 0) brother = -1; } } - if(brother >= i) + if(brother >= pi) { pf2ptr->SetPointIndex(brother); pf2ptr->SetNV(*nv[brother-1]); } } - if(usesum && brother < i) + if(usesum && brother < pi) continue; //pf.UnSetNV(); x = 0; //(*testout) << "before " << pf.Func(x); - pf.SetNV(*nv[i-1]); + pf.SetNV(*nv[pi-1]); x = 0; int pok = (brother == -1) ? (pf.Func (x) < 1e10) : (pf_sum.Func (x) < 1e10); @@ -1757,12 +1759,12 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp, BFGS (x, pf_sum, par); - for(j=0; j<3; j++) - points.Elem(i)(j) += x(j);// - scal*nv[i-1].X(j); + for(int j=0; j<3; j++) + points[pi](j) += x(j);// - scal*nv[i-1].X(j); if(brother != -1) - for(j=0; j<3; j++) - points.Elem(brother)(j) += x(j);// - scal*nv[brother-1].X(j); + for(int j=0; j<3; j++) + points[brother](j) += x(j);// - scal*nv[brother-1].X(j); } @@ -1780,7 +1782,7 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp, PrintDot ('\n'); delete pf2ptr; - for(i=0; i