ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-11-11 16:49:16 +05:00
Code Issues Packages Projects Releases Wiki Activity

T_POINTS are now ngcore::Array

Browse Source
This commit is contained in:
Joachim Schöberl 2019-08-10 00:21:37 +02:00
parent 7b9fbbf0b1
commit d215ac1025
8 changed files with 186 additions and 147 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

82
libsrc/core/array.hpp
View File

@ -59,6 +59,28 @@ namespace ngcore
}
template <typename AO>
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,25 +121,10 @@ namespace ngcore
// NETGEN_INLINE auto & operator[] (size_t i) { return Spec()[i]; }
NETGEN_INLINE auto operator[] (size_t i) const { return Spec()[i]; }
};
template <typename AO>
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; }
// NETGEN_INLINE auto begin() const { return Spec().begin(); }
// NETGEN_INLINE auto end() const { return Spec().end(); }
NETGEN_INLINE auto begin () const { return AOWrapperIterator<BaseArrayObject> (*this, 0); }
NETGEN_INLINE auto end () const { return AOWrapperIterator<BaseArrayObject> (*this, Size()); }
};
@ -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<T> begin() const { return first; }
NETGEN_INLINE ArrayRangeIterator<T> 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<T(int)> & 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<size_t> Range () const
{
return T_Range<size_t> (0, Size());
return T_Range<size_t> (BASE, Size()+BASE);
}
NETGEN_INLINE const CArray<T> Addr (size_t pos) const
{
return CArray<T> (data+pos);
return CArray<T> (data+pos-BASE);
}
// const CArray<T> operator+ (int pos)
@ -606,6 +615,7 @@ namespace ngcore
using FlatArray<T,IndexType>::size;
using FlatArray<T,IndexType>::data;
using FlatArray<T,IndexType>::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<T> 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<T>::operator= (val);
FlatArray<T,IndexType>::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;
}

12
libsrc/meshing/findip.hpp
View File

@ -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;
}
}

28
libsrc/meshing/improve2.cpp
View File

@ -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<SurfaceElementIndex,PointIndex::BASE> elementsonnode(np);
NgArray<SurfaceElementIndex> 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<bool,PointIndex::BASE> fixed(np);
Array<bool,PointIndex> 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<Vec<3>,PointIndex::BASE> normals(np);

37
libsrc/meshing/meshclass.cpp
View File

@ -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,8 +2415,12 @@ 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++)
{
@ -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<BoxTree<3>> (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<BoxTree<3>> (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);
}
}

22
libsrc/meshing/meshtype.cpp
View File

@ -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();

32
libsrc/meshing/meshtype.hpp
View File

@ -369,8 +369,8 @@ namespace netgen
typedef NgArray<MeshPoint, PointIndex::BASE, PointIndex> T_POINTS;
// typedef Array<MeshPoint, PointIndex> T_POINTS;
// typedef NgArray<MeshPoint, PointIndex::BASE, PointIndex> T_POINTS;
typedef Array<MeshPoint, PointIndex> T_POINTS;
@ -500,18 +500,12 @@ namespace netgen
///
const PointIndex & operator[] (int i) const { return pnum[i]; }
NgFlatArray<const PointIndex> PNums () const
{ return NgFlatArray<const PointIndex> (np, &pnum[0]); }
NgFlatArray<PointIndex> PNums ()
{ return NgFlatArray<PointIndex> (np, &pnum[0]); }
auto Vertices() const
{ return NgFlatArray<const PointIndex> (GetNV(), &pnum[0]); }
NgFlatArray<const PointGeomInfo> GeomInfo() const
{ return NgFlatArray<const PointGeomInfo> (np, &geominfo[0]); }
NgFlatArray<PointGeomInfo> GeomInfo()
{ return NgFlatArray<PointGeomInfo> (np, &geominfo[0]); }
auto PNums () const { return FlatArray<const PointIndex> (np, &pnum[0]); }
auto PNums () { return FlatArray<PointIndex> (np, &pnum[0]); }
auto Vertices() const { return FlatArray<const PointIndex> (GetNV(), &pnum[0]); }
auto GeomInfo() const { return FlatArray<const PointGeomInfo> (np, &geominfo[0]); }
auto GeomInfo() { return FlatArray<PointGeomInfo> (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<const PointIndex> PNums () const
{ return NgFlatArray<const PointIndex> (np, &pnum[0]); }
auto PNums () const
{ return FlatArray<const PointIndex> (np, &pnum[0]); }
NgFlatArray<const PointIndex> Vertices() const { return { GetNV(), &pnum[0] }; }
FlatArray<const PointIndex> Vertices() const { return { GetNV(), &pnum[0] }; }
///
PointIndex & PNum (int i) { return pnum[i-1]; }

8
libsrc/meshing/parallelmesh.cpp
View File

@ -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<const PointIndex> 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<const PointIndex> 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<int> pp_data;
MyMPI_Recv(pp_data, 0, MPI_TAG_MESH+1, comm);

108
libsrc/meshing/smoothing3.cpp
View File

@ -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<double, PointIndex::BASE> pointh (points.Size());
NgArray<double, PointIndex::BASE, PointIndex> 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<GetNE(); i++)
for (int i=0; i<GetNE(); i++)
{
const Element & el = VolumeElement(i+1);
double h = pow(el.Volume(points),1./3.);
for(j=1; j<=el.GetNV(); j++)
for(int j=1; j<=el.GetNV(); j++)
if(h > 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<int,PointIndex::BASE>* > * 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<GetNE(); i++)
for(int i=0; i<GetNE(); i++)
{
const Element & el = VolumeElement(i+1);
double h = pow(el.Volume(points),1./3.);
for(j=1; j<=el.GetNV(); j++)
for(int j=1; j<=el.GetNV(); j++)
if(h > 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 && j<used_idmaps->Size(); j++)
for(int j=0; brother == -1 && j<used_idmaps->Size(); 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<locidmaps.Size(); i++)
for(int i=0; i<locidmaps.Size(); i++)
delete locidmaps[i];
multithread.task = savetask;