ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-01-27 13:20:34 +05:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' into stl_to_python

Browse Source
This commit is contained in:
Christopher Lackner 2016-12-12 10:39:51 +01:00
commit 6e60338987
46 changed files with 834 additions and 613 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
View File

@ -172,7 +172,6 @@ if (USE_GUI)
get_filename_component(MY_LIB_DIR ${TK_LIBRARY} DIRECTORY) 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}/tcl8.6" DESTINATION lib COMPONENT netgen )
install( DIRECTORY "${MY_LIB_DIR}/tk8.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 ) install( DIRECTORY "${MY_LIB_DIR}/../bin" DESTINATION . COMPONENT netgen )
add_definitions(-DTOGL_WGL) add_definitions(-DTOGL_WGL)
else(WIN32) else(WIN32)

3
libsrc/csg/csgeom.cpp
View File

@ -648,7 +648,8 @@ namespace netgen
int CSGeometry :: SetTopLevelObject (Solid * sol, Surface * surf) 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 :: TopLevelObject * CSGeometry ::

12
libsrc/csg/edgeflw.cpp
View File

@ -617,7 +617,9 @@ namespace netgen
(*testout) << "Point on edge" << endl (*testout) << "Point on edge" << endl
<< "seg = " << i2 << ", p = " << pi << endl << "seg = " << i2 << ", p = " << pi << endl
<< "pos = " << p << ", projected = " << hp << 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 (seg.edgenr >= 1 && seg.edgenr <= cntedge)
{ {
if (osedges.Get(seg.edgenr) != -1) if (osedges.Get(seg.edgenr) != -1)
@ -1156,7 +1159,8 @@ namespace netgen
//seg.surfnr2 = s2_rep; //seg.surfnr2 = s2_rep;
seg.surfnr1 = s1; seg.surfnr1 = s1;
seg.surfnr2 = s2; seg.surfnr2 = s2;
hi = refedges.Append (seg); refedges.Append (seg);
hi = refedges.Size();
refedgesinv.Append (edgeinv); refedgesinv.Append (edgeinv);
edges_priority.Append((pre_ok[k-1]) ? 1 : 0); edges_priority.Append((pre_ok[k-1]) ? 1 : 0);
} }

6
libsrc/csg/genmesh.cpp
View File

@ -759,8 +759,9 @@ namespace netgen
<< "time = " << GetTime() << " sec" << endl << "time = " << GetTime() << " sec" << endl
<< "points: " << mesh->GetNP() << endl; << "points: " << mesh->GetNP() << endl;
#endif #endif
if (mparam.uselocalh && 0) /*
if (mparam.uselocalh)
{ {
mesh->CalcLocalH(mparam.grading); mesh->CalcLocalH(mparam.grading);
mesh->DeleteMesh(); mesh->DeleteMesh();
@ -773,6 +774,7 @@ namespace netgen
MeshSurface (geom, *mesh, mparam); MeshSurface (geom, *mesh, mparam);
if (multithread.terminate) return TCL_OK; if (multithread.terminate) return TCL_OK;
} }
*/
#ifdef LOG_STREAM #ifdef LOG_STREAM
(*logout) << "Surfaces remeshed" << endl (*logout) << "Surfaces remeshed" << endl

3
libsrc/csg/polyhedra.cpp
View File

@ -499,7 +499,8 @@ int Polyhedra :: AddPoint (const Point<3> & p)
else else
poly_bbox.Add(p); 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) int Polyhedra :: AddFace (int pi1, int pi2, int pi3, int inputnum)

3
libsrc/csg/specpoin.cpp
View File

@ -2027,7 +2027,8 @@ namespace netgen
if (spi == -1) if (spi == -1)
{ {
spi = specpoints.Append (SpecialPoint()) - 1; specpoints.Append (SpecialPoint());
spi = specpoints.Size()-1;
specpoint2point.Append (i); specpoint2point.Append (i);
specpoints.Last().unconditional = 0; specpoints.Last().unconditional = 0;
searchtree.Insert (apoints[i], spi); searchtree.Insert (apoints[i], spi);

40
libsrc/csg/vscsg.cpp
View File

@ -130,18 +130,18 @@ namespace netgen
int hasp = 0; int hasp = 0;
for (int i = 0; i < geometry->GetNTopLevelObjects(); i++) for (int i = 0; i < geometry->GetNTopLevelObjects(); i++)
{ {
const TriangleApproximation & ta = const TriangleApproximation * ta =
*geometry->GetTriApprox(i); geometry->GetTriApprox(i);
if (!&ta) continue; if (!ta) continue;
for (int j = 0; j < ta.GetNP(); j++) for (int j = 0; j < ta->GetNP(); j++)
{ {
if (hasp) if (hasp)
box.Add (ta.GetPoint(j)); box.Add (ta->GetPoint(j));
else else
{ {
hasp = 1; hasp = 1;
box.Set (ta.GetPoint(j)); box.Set (ta->GetPoint(j));
} }
} }
} }
@ -167,18 +167,18 @@ namespace netgen
trilists.Append (glGenLists (1)); trilists.Append (glGenLists (1));
glNewList (trilists.Last(), GL_COMPILE); glNewList (trilists.Last(), GL_COMPILE);
glEnable (GL_NORMALIZE); glEnable (GL_NORMALIZE);
const TriangleApproximation & ta = const TriangleApproximation * ta =
*geometry->GetTriApprox(i); geometry->GetTriApprox(i);
if (&ta) if (ta)
{ {
glEnableClientState(GL_VERTEX_ARRAY); 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); 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++) for (int j = 0; j < ta->GetNT(); j++)
glDrawElements(GL_TRIANGLES, 3, GL_UNSIGNED_INT, & (ta.GetTriangle(j)[0])); glDrawElements(GL_TRIANGLES, 3, GL_UNSIGNED_INT, & (ta->GetTriangle(j)[0]));
glDisableClientState(GL_VERTEX_ARRAY); glDisableClientState(GL_VERTEX_ARRAY);
glDisableClientState(GL_NORMAL_ARRAY); glDisableClientState(GL_NORMAL_ARRAY);
@ -382,13 +382,20 @@ namespace netgen
glDisable (GL_COLOR_MATERIAL); glDisable (GL_COLOR_MATERIAL);
glDisable (GL_LIGHTING); glDisable (GL_LIGHTING);
glDisable (GL_CLIP_PLANE0); glDisable (GL_CLIP_PLANE0);
/*
for (int i = 1; i <= mesh -> GetNP(); i++) for (int i = 1; i <= mesh -> GetNP(); i++)
{ {
const Point3d & p = mesh -> Point(i); const Point3d & p = mesh -> Point(i);
glRasterPos3d (p.X(), p.Y(), p.Z()); glRasterPos3d (p.X(), p.Y(), p.Z());
glBitmap (7, 7, 3, 3, 0, 0, &knoedel[0]); 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) if (vispar.drawedpointnrs)
@ -403,12 +410,13 @@ namespace netgen
// glListBase (fontbase); // glListBase (fontbase);
char buf[20]; 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); const Point3d & p = mesh->Point(i);
glRasterPos3d (p.X(), p.Y(), p.Z()); glRasterPos3d (p.X(), p.Y(), p.Z());
sprintf (buf, "%d", i); sprintf (buf, "%d", int(i));
// glCallLists (GLsizei(strlen (buf)), GL_UNSIGNED_BYTE, buf); // glCallLists (GLsizei(strlen (buf)), GL_UNSIGNED_BYTE, buf);
MyOpenGLText (buf); MyOpenGLText (buf);
} }

4
libsrc/general/array.hpp
View File

@ -293,13 +293,13 @@ namespace netgen
/// Add element at end of array. reallocation if necessary. /// Add element at end of array. reallocation if necessary.
int Append (const T & el) void Append (const T & el)
{ {
if (size == allocsize) if (size == allocsize)
ReSize (size+1); ReSize (size+1);
data[size] = el; data[size] = el;
size++; size++;
return size; // return size;
} }
template <typename T2, int B2> template <typename T2, int B2>

3
libsrc/include/mydefs.hpp
View File

@ -38,7 +38,8 @@
// #define BASE0
// #define DEBUG
#define noDEMOVERSION #define noDEMOVERSION

33
libsrc/include/nginterface_v2.hpp
View File

@ -14,15 +14,20 @@
namespace netgen namespace netgen
{ {
static constexpr int POINTINDEX_BASE = 1;
struct T_EDGE2 struct T_EDGE2
{ {
int orient:1; // int orient:1;
int nr:31; // 0-based // int nr:31; // 0-based
int nr; // 0-based
}; };
struct T_FACE2 struct T_FACE2
{ {
int orient:3; // int orient:3;
int nr:29; // 0-based // int nr:29; // 0-based
int nr; // 0-based
}; };
class Ng_Element class Ng_Element
@ -35,7 +40,7 @@ namespace netgen
const int * ptr; const int * ptr;
int Size() const { return num; } 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; const int * ptr;
int Size() const { return num; } 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 class Ng_Edges
@ -71,6 +76,17 @@ namespace netgen
int operator[] (int i) const { return ptr[i].nr; } 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: public:
NG_ELEMENT_TYPE type; NG_ELEMENT_TYPE type;
int index; // material / boundary condition int index; // material / boundary condition
@ -81,6 +97,7 @@ namespace netgen
Ng_Vertices vertices; Ng_Vertices vertices;
Ng_Edges edges; Ng_Edges edges;
Ng_Faces faces; Ng_Faces faces;
Ng_Facets facets;
bool is_curved; bool is_curved;
}; };
@ -131,7 +148,7 @@ namespace netgen
const int * ptr; const int * ptr;
int Size() const { return 2; } 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; const int * ptr;
int Size() const { return nv; } 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 class Ng_Edges

29
libsrc/include/nginterface_v2_impl.hpp
View File

@ -1,6 +1,6 @@
NGX_INLINE DLL_HEADER Ng_Point Ngx_Mesh :: GetPoint (int nr) const 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.num = 0;
ret.faces.ptr = NULL; 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 = mesh->GetCurvedElements().IsSegmentCurved(nr);
ret.is_curved = el.IsCurved(); 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.num = MeshTopology::GetNFaces (el.GetType());
ret.faces.ptr = (T_FACE2*)mesh->GetTopology().GetSurfaceElementFacesPtr (nr); 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(); ret.is_curved = el.IsCurved();
return ret; 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.num = MeshTopology::GetNFaces (el.GetType());
ret.faces.ptr = (T_FACE2*)mesh->GetTopology().GetElementFacesPtr (nr); 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(); ret.is_curved = el.IsCurved();
return ret; return ret;
} }

2
libsrc/interface/writejcm.cpp
View File

@ -33,7 +33,7 @@ void WriteJCMFormat (const Mesh & mesh,
int np = mesh.GetNP(); int np = mesh.GetNP();
// Identic points // Identic points
Array<int,1> identmap1, identmap2, identmap3; Array<int,PointIndex::BASE> identmap1, identmap2, identmap3;
mesh.GetIdentifications().GetMap(1, identmap1); mesh.GetIdentifications().GetMap(1, identmap1);
mesh.GetIdentifications().GetMap(2, identmap2); mesh.GetIdentifications().GetMap(2, identmap2);
mesh.GetIdentifications().GetMap(3, identmap3); mesh.GetIdentifications().GetMap(3, identmap3);

6
libsrc/interface/wuchemnitz.cpp
View File

@ -202,7 +202,8 @@ namespace netgen
fa.p1 = i3.I1(); fa.p1 = i3.I1();
fa.p2 = i3.I2(); fa.p2 = i3.I2();
fa.p3 = i3.I3(); fa.p3 = i3.I3();
facei = faces.Append (fa); faces.Append (fa);
facei = faces.Size();
faceindex.Set (i3, facei); faceindex.Set (i3, facei);
} }
@ -241,7 +242,8 @@ namespace netgen
EDGE ed; EDGE ed;
ed.p1 = i2.I1(); ed.p1 = i2.I1();
ed.p2 = i2.I2(); ed.p2 = i2.I2();
edgei = edges.Append (ed); edges.Append (ed);
edgei = edges.Size();
edgeindex.Set (i2, edgei); edgeindex.Set (i2, edgei);
} }

12
libsrc/meshing/adfront2.cpp
View File

@ -88,7 +88,8 @@ namespace netgen
} }
else else
{ {
pi = points.Append (FrontPoint2 (p, globind, mgi, pointonsurface)) - 1; points.Append (FrontPoint2 (p, globind, mgi, pointonsurface));
pi = points.Size()-1;
} }
if (mgi) if (mgi)
@ -128,7 +129,8 @@ namespace netgen
} }
else 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); pindex.Append (pi);
invpindex[pi] = pindex.Size(); invpindex[pi] = pindex.Size();
loclines[i][j] = locpoints.Append (points[pi].P()); locpoints.Append (points[pi].P());
loclines[i][j] = locpoints.Size();
} }
else else
loclines[i][j] = invpindex[pi]; loclines[i][j] = invpindex[pi];
@ -349,7 +352,8 @@ namespace netgen
// Dist2 (points.Get(i).P(), p0) <= xh2 && // Dist2 (points.Get(i).P(), p0) <= xh2 &&
invpindex[i] <= 0) invpindex[i] <= 0)
{ {
invpindex[i] = locpoints.Append (points[i].P()); locpoints.Append (points[i].P());
invpindex[i] = locpoints.Size();
pindex.Append(i); pindex.Append(i);
} }
} }

2
libsrc/meshing/adfront2.hpp
View File

@ -35,7 +35,7 @@
/// ///
FrontPoint2 () FrontPoint2 ()
{ {
globalindex = -1; globalindex.Invalidate(); // = -1;
nlinetopoint = 0; nlinetopoint = 0;
frontnr = INT_MAX-10; // attention: overflow on calculating INT_MAX + 1 frontnr = INT_MAX-10; // attention: overflow on calculating INT_MAX + 1
mgi = NULL; mgi = NULL;

28
libsrc/meshing/adfront3.cpp
View File

@ -9,7 +9,7 @@ namespace netgen
FrontPoint3 :: FrontPoint3 () FrontPoint3 :: FrontPoint3 ()
{ {
globalindex = -1; globalindex.Invalidate(); // = -1;
nfacetopoint = 0; nfacetopoint = 0;
frontnr = 1000; frontnr = 1000;
cluster = 0; cluster = 0;
@ -159,8 +159,9 @@ INDEX AdFront3 :: AddFace (const MiniElement2d & aface)
for (i = 1; i <= aface.GetNP(); i++) for (i = 1; i <= aface.GetNP(); i++)
points[aface.PNum(i)].DecFrontNr (minfn+1); 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; faces.Elem(nfn).cluster = cluster;
if (hashon && hashcreated) if (hashon && hashcreated)
@ -484,9 +485,9 @@ int AdFront3 :: SelectBaseElement ()
int AdFront3 :: GetLocals (int fstind, int AdFront3 :: GetLocals (int fstind,
Array<Point3d > & locpoints, Array<Point3d, PointIndex::BASE> & locpoints,
Array<MiniElement2d> & locfaces, // local index Array<MiniElement2d> & locfaces, // local index
Array<PointIndex> & pindex, Array<PointIndex, PointIndex::BASE> & pindex,
Array<INDEX> & findex, Array<INDEX> & findex,
INDEX_2_HASHTABLE<int> & getconnectedpairs, INDEX_2_HASHTABLE<int> & getconnectedpairs,
float xh, float xh,
@ -600,12 +601,13 @@ int AdFront3 :: GetLocals (int fstind,
if (invpindex[pi] == -1) if (invpindex[pi] == -1)
{ {
pindex.Append (pi); pindex.Append (pi);
invpindex[pi] = pindex.Size(); // -1+PointIndex::BASE; locpoints.Append (points[pi].P());
locfaces.Elem(i).PNum(j) = locpoints.Append (points[pi].P()); invpindex[pi] = pindex.Size()-1+PointIndex::BASE;
} }
else // locfaces.Elem(i).PNum(j) = locpoints.Append (points[pi].P());
locfaces.Elem(i).PNum(j) = invpindex[pi]; // }
// else
locfaces.Elem(i).PNum(j) = invpindex[pi];
} }
} }
@ -655,9 +657,9 @@ int AdFront3 :: GetLocals (int fstind,
// returns all points connected with fi // returns all points connected with fi
void AdFront3 :: GetGroup (int fi, void AdFront3 :: GetGroup (int fi,
Array<MeshPoint> & grouppoints, Array<MeshPoint, PointIndex::BASE> & grouppoints,
Array<MiniElement2d> & groupelements, Array<MiniElement2d> & groupelements,
Array<PointIndex> & pindex, Array<PointIndex, PointIndex::BASE> & pindex,
Array<INDEX> & findex) Array<INDEX> & findex)
{ {
// static Array<char> pingroup; // static Array<char> pingroup;

68
libsrc/meshing/adfront3.hpp
View File

@ -17,15 +17,15 @@
class FrontPoint3 class FrontPoint3
{ {
/// coordinates /// coordinates
Point<3> p; Point<3> p;
/// global node index /// global node index
PointIndex globalindex; PointIndex globalindex;
/// number of faces connected to point /// number of faces connected to point
int nfacetopoint; int nfacetopoint;
/// distance to original boundary /// distance to original boundary
int frontnr; int frontnr;
/// ///
int cluster; int cluster;
public: public:
/// ///
FrontPoint3 (); FrontPoint3 ();
@ -43,7 +43,7 @@ public:
void AddFace () void AddFace ()
{ nfacetopoint++; } { nfacetopoint++; }
/// /// if last face is removed, then point is invalidated
void RemoveFace() void RemoveFace()
{ {
nfacetopoint--; nfacetopoint--;
@ -51,7 +51,7 @@ public:
} }
/// ///
int Valid () const bool Valid () const
{ return nfacetopoint >= 0; } { return nfacetopoint >= 0; }
/// ///
@ -74,7 +74,7 @@ class MiniElement2d
{ {
protected: protected:
int np; int np;
PointIndex pnum[4]; PointIndex pnum[4]; // can be global or local nums
bool deleted; bool deleted;
public: public:
MiniElement2d () MiniElement2d ()
@ -89,8 +89,8 @@ public:
const PointIndex PNum (int i) const { return pnum[i-1]; } const PointIndex PNum (int i) const { return pnum[i-1]; }
PointIndex & PNum (int i) { return pnum[i-1]; } PointIndex & PNum (int i) { return pnum[i-1]; }
const PointIndex PNumMod (int i) const { return pnum[(i-1)%np]; } const PointIndex PNumMod (int i) const { return pnum[(i-1)%np]; }
auto PNums() const { return FlatArray<const PointIndex> (np, &pnum[0]); }
void Delete () { deleted = 1; pnum[0] = pnum[1] = pnum[2] = pnum[3] = PointIndex::BASE-1; } void Delete () { deleted = true; for (PointIndex & p : pnum) p.Invalidate(); }
bool IsDeleted () const { return deleted; } bool IsDeleted () const { return deleted; }
}; };
@ -120,7 +120,7 @@ private:
int hashvalue; int hashvalue;
/// ///
int cluster; int cluster;
public: public:
/// ///
FrontFace (); FrontFace ();
@ -178,43 +178,43 @@ class AdFront3
/// ///
Array<FrontPoint3, PointIndex::BASE, PointIndex> points; Array<FrontPoint3, PointIndex::BASE, PointIndex> points;
/// ///
Array<FrontFace> faces; Array<FrontFace> faces;
/// ///
Array<PointIndex> delpointl; Array<PointIndex> delpointl;
/// which points are connected to pi ? /// which points are connected to pi ?
TABLE<int, PointIndex::BASE> * connectedpairs; TABLE<int, PointIndex::BASE> * connectedpairs;
/// number of total front faces; /// number of total front faces;
int nff; int nff;
/// number of quads in front /// 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 /// counter for rebuilding internal tables
int rebuildcounter; int rebuildcounter;
/// last base element /// last base element
int lasti; int lasti;
/// minimal selection-value of baseelements /// minimal selection-value of baseelements
int minval; int minval;
Array<int, PointIndex::BASE, PointIndex> invpindex; Array<PointIndex, PointIndex::BASE, PointIndex> invpindex;
Array<char> pingroup; Array<char> pingroup;
/// ///
class Box3dTree * facetree; class Box3dTree * facetree;
public: public:
/// ///
AdFront3 (); AdFront3 ();
/// ///
@ -260,9 +260,9 @@ public:
/// ///
int GetLocals (int baseelement, int GetLocals (int baseelement,
Array<Point3d > & locpoints, Array<Point3d, PointIndex::BASE> & locpoints,
Array<MiniElement2d> & locfaces, // local index Array<MiniElement2d> & locfaces, // local index
Array<PointIndex> & pindex, Array<PointIndex, PointIndex::BASE> & pindex,
Array<INDEX> & findex, Array<INDEX> & findex,
INDEX_2_HASHTABLE<int> & connectedpairs, INDEX_2_HASHTABLE<int> & connectedpairs,
float xh, float xh,
@ -271,9 +271,9 @@ public:
/// ///
void GetGroup (int fi, void GetGroup (int fi,
Array<MeshPoint> & grouppoints, Array<MeshPoint, PointIndex::BASE> & grouppoints,
Array<MiniElement2d> & groupelements, Array<MiniElement2d> & groupelements,
Array<PointIndex> & pindex, Array<PointIndex, PointIndex::BASE> & pindex,
Array<INDEX> & findex); Array<INDEX> & findex);
/// ///

31
libsrc/meshing/clusters.cpp
View File

@ -96,7 +96,7 @@ namespace netgen
nnums.SetSize(elnv+elned+elnfa+1); nnums.SetSize(elnv+elned+elnfa+1);
for (int j = 1; j <= elnv; j++) 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++) for (int j = 1; j <= elned; j++)
nnums.Elem(elnv+j) = nv+ednums.Elem(j); nnums.Elem(elnv+j) = nv+ednums.Elem(j);
for (int j = 1; j <= elnfa; j++) for (int j = 1; j <= elnfa; j++)
@ -108,7 +108,6 @@ namespace netgen
} }
}); });
NgProfiler::StopTimer(timer1); NgProfiler::StopTimer(timer1);
NgProfiler::StartTimer(timer2); NgProfiler::StartTimer(timer2);
@ -125,7 +124,7 @@ namespace netgen
nnums.SetSize(elnv+elned+1); nnums.SetSize(elnv+elned+1);
for (int j = 1; j <= elnv; j++) 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++) for (int j = 1; j <= elned; j++)
nnums.Elem(elnv+j) = nv+ednums.Elem(j); nnums.Elem(elnv+j) = nv+ednums.Elem(j);
nnums.Elem(elnv+elned+1) = fanum; nnums.Elem(elnv+elned+1) = fanum;
@ -197,7 +196,7 @@ namespace netgen
nnums.SetSize(elnv+elned+elnfa+1); nnums.SetSize(elnv+elned+elnfa+1);
for (int j = 1; j <= elnv; j++) 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++) for (int j = 1; j <= elned; j++)
nnums.Elem(elnv+j) = nv+ednums.Elem(j); nnums.Elem(elnv+j) = nv+ednums.Elem(j);
for (int j = 1; j <= elnfa; j++) for (int j = 1; j <= elnfa; j++)
@ -220,23 +219,23 @@ namespace netgen
break; break;
case TET: case TET:
case TET10: case TET10:
if (cluster_reps.Get(el.PNum(1)) == if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) ==
cluster_reps.Get(el.PNum(2))) cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE))
clustertab = tet_cluster12; clustertab = tet_cluster12;
else if (cluster_reps.Get(el.PNum(1)) == else if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) ==
cluster_reps.Get(el.PNum(3))) cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE))
clustertab = tet_cluster13; clustertab = tet_cluster13;
else if (cluster_reps.Get(el.PNum(1)) == else if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) ==
cluster_reps.Get(el.PNum(4))) cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE))
clustertab = tet_cluster14; clustertab = tet_cluster14;
else if (cluster_reps.Get(el.PNum(2)) == else if (cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE) ==
cluster_reps.Get(el.PNum(3))) cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE))
clustertab = tet_cluster23; clustertab = tet_cluster23;
else if (cluster_reps.Get(el.PNum(2)) == else if (cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE) ==
cluster_reps.Get(el.PNum(4))) cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE))
clustertab = tet_cluster24; clustertab = tet_cluster24;
else if (cluster_reps.Get(el.PNum(3)) == else if (cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE) ==
cluster_reps.Get(el.PNum(4))) cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE))
clustertab = tet_cluster34; clustertab = tet_cluster34;
else else

2
libsrc/meshing/curvedelems.cpp
View File

@ -3873,7 +3873,7 @@ namespace netgen
Mat<DIM_SPACE,2,T> _dxdxi; Mat<DIM_SPACE,2,T> _dxdxi;
if (!EvaluateMapping (info, _xi, _x, _dxdxi)) if (!EvaluateMapping (info, _xi, _x, _dxdxi))
{ ok = false; break; } { ok = false; break; }
// cout << "x = " << _x << ", dxdxi = " << _dxdxi << endl; // *testout << "x = " << _x << ", dxdxi = " << _dxdxi << endl;
if (x) if (x)
for (int j = 0; j < DIM_SPACE; j++) for (int j = 0; j < DIM_SPACE; j++)
x[i*sx+j] = _x[j]; x[i*sx+j] = _x[j];

22
libsrc/meshing/delaunay.cpp
View File

@ -415,9 +415,9 @@ namespace netgen
INDEX_3 i3 = tempels.Get(helind).GetFace (k-1); INDEX_3 i3 = tempels.Get(helind).GetFace (k-1);
const Point3d & p1 = mesh.Point ( i3.I1()); const Point3d & p1 = mesh.Point ( PointIndex (i3.I1()) );
const Point3d & p2 = mesh.Point ( i3.I2()); const Point3d & p2 = mesh.Point ( PointIndex (i3.I2()) );
const Point3d & p3 = mesh.Point ( i3.I3()); const Point3d & p3 = mesh.Point ( PointIndex (i3.I3()) );
Vec3d v1(p1, p2); Vec3d v1(p1, p2);
@ -612,20 +612,20 @@ namespace netgen
for (int i = 1; i <= adfront->GetNF(); i++) for (int i = 1; i <= adfront->GetNF(); i++)
{ {
const MiniElement2d & face = adfront->GetFace(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])); pmin.SetToMin (mesh.Point (pi));
pmax.SetToMax (mesh.Point (face[j])); 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])); pmin.SetToMin (mesh.Point (pi));
pmax.SetToMax (mesh.Point (mesh.LockedPoints()[i])); pmax.SetToMax (mesh.Point (pi));
} }
Vec3d vdiag(pmin, pmax); Vec3d vdiag(pmin, pmax);
// double r1 = vdiag.Length(); // double r1 = vdiag.Length();

235
libsrc/meshing/improve3.cpp
View File

@ -274,16 +274,8 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
void MeshOptimize3d :: SplitImprove (Mesh & mesh, void MeshOptimize3d :: SplitImprove (Mesh & mesh,
OPTIMIZEGOAL goal) OPTIMIZEGOAL goal)
{ {
int j, k, l;
Point3d p1, p2, pnew;
ElementIndex ei;
SurfaceElementIndex sei;
PointIndex pi1, pi2;
double bad1, bad2, badmax, badlimit; double bad1, bad2, badmax, badlimit;
int cnt = 0; int cnt = 0;
int np = mesh.GetNP(); int np = mesh.GetNP();
int ne = mesh.GetNE(); int ne = mesh.GetNE();
@ -291,7 +283,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
TABLE<ElementIndex,PointIndex::BASE> elementsonnode(np); TABLE<ElementIndex,PointIndex::BASE> elementsonnode(np);
Array<ElementIndex> hasbothpoints; Array<ElementIndex> hasbothpoints;
BitArray origpoint(np), boundp(np); BitArray origpoint(np+1), boundp(np+1); // big enough for 0 and 1-based
origpoint.Set(); origpoint.Set();
Array<double> elerrs(ne); Array<double> elerrs(ne);
@ -301,7 +293,6 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
const char * savetask = multithread.task; const char * savetask = multithread.task;
multithread.task = "Split Improve"; multithread.task = "Split Improve";
PrintMessage (3, "SplitImprove"); PrintMessage (3, "SplitImprove");
(*testout) << "start SplitImprove" << "\n"; (*testout) << "start SplitImprove" << "\n";
@ -311,10 +302,11 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
bad1 = 0; bad1 = 0;
badmax = 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; continue;
elerrs[ei] = CalcBad (mesh.Points(), mesh[ei], 0); elerrs[ei] = CalcBad (mesh.Points(), mesh[ei], 0);
bad1 += elerrs[ei]; bad1 += elerrs[ei];
if (elerrs[ei] > badmax) badmax = elerrs[ei]; if (elerrs[ei] > badmax) badmax = elerrs[ei];
@ -325,41 +317,40 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
boundp.Clear(); boundp.Clear();
for (sei = 0; sei < mesh.GetNSE(); sei++) for (auto & el : mesh.SurfaceElements())
for (j = 0; j < 3; j++) for (PointIndex pi : el.PNums())
boundp.Set (mesh[sei][j]); boundp.Set (pi);
if (goal == OPT_QUALITY) if (goal == OPT_QUALITY)
{ {
bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements()); bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
(*testout) << "Total badness = " << bad1 << endl; (*testout) << "Total badness = " << bad1 << endl;
} }
for (ei = 0; ei < ne; ei++) for (ElementIndex ei : mesh.VolumeElements().Range())
for (j = 0; j < mesh[ei].GetNP(); j++) for (PointIndex pi : mesh[ei].PNums())
elementsonnode.Add (mesh[ei][j], ei); elementsonnode.Add (pi, ei);
mesh.MarkIllegalElements(); mesh.MarkIllegalElements();
if (goal == OPT_QUALITY || goal == OPT_LEGAL) if (goal == OPT_QUALITY || goal == OPT_LEGAL)
{ {
int cntill = 0; int cntill = 0;
for (ei = 0; ei < ne; ei++) for (ElementIndex ei = 0; ei < ne; ei++)
{ {
// if (!LegalTet (volelements.Get(i))) // if (!LegalTet (volelements.Get(i)))
if (mesh[ei].flags.illegal) if (mesh[ei].flags.illegal)
{ {
cntill++; cntill++;
illegaltet.Set (ei+1); illegaltet.Set (ei);
} }
} }
// (*mycout) << cntill << " illegal tets" << endl;
} }
for (ElementIndex ei : mesh.VolumeElements().Range())
for (ei = 0; ei < ne; ei++)
{ {
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; continue;
if (multithread.terminate) if (multithread.terminate)
break; break;
@ -368,82 +359,74 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
bool ltestmode = 0; bool ltestmode = 0;
if (elerrs[ei] < badlimit && !illegaltet.Test(ei)) continue;
if (elerrs[ei] < badlimit && !illegaltet.Test(ei+1)) continue;
if ((goal == OPT_LEGAL) && if ((goal == OPT_LEGAL) &&
!illegaltet.Test(ei+1) && !illegaltet.Test(ei) &&
CalcBad (mesh.Points(), mesh[ei], 0) < 1e3) CalcBad (mesh.Points(), elem, 0) < 1e3)
continue; continue;
Element & elem = mesh[ei];
if (ltestmode) if (ltestmode)
{ {
(*testout) << "test el " << ei << endl; (*testout) << "test el " << ei << endl;
for (j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
(*testout) << elem[j] << " "; (*testout) << elem[j] << " ";
(*testout) << endl; (*testout) << endl;
} }
for (int j = 0; j < 6; j++)
for (j = 0; j < 6; j++)
{ {
static const int tetedges[6][2] = static const int tetedges[6][2] =
{ { 0, 1 }, { 0, 2 }, { 0, 3 }, { { 0, 1 }, { 0, 2 }, { 0, 3 },
{ 1, 2 }, { 1, 3 }, { 2, 3 } }; { 1, 2 }, { 1, 3 }, { 2, 3 } };
pi1 = elem[tetedges[j][0]]; PointIndex pi1 = elem[tetedges[j][0]];
pi2 = elem[tetedges[j][1]]; PointIndex pi2 = elem[tetedges[j][1]];
if (pi2 < pi1) Swap (pi1, pi2); 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)) if (!origpoint.Test(pi1) || !origpoint.Test(pi2))
continue; continue;
INDEX_2 i2(pi1, pi2); INDEX_2 i2(pi1, pi2);
i2.Sort(); i2.Sort();
if (mesh.BoundaryEdge (pi1, pi2)) continue; 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); edgetested.Set (i2, 1);
hasbothpoints.SetSize (0); 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); 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 (has1 && has2)
if (!hasbothpoints.Contains (elnr)) if (!hasbothpoints.Contains (ei))
hasbothpoints.Append (elnr); hasbothpoints.Append (ei);
} }
bad1 = 0; 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; bool puretet = 1;
for (k = 0; k < hasbothpoints.Size(); k++) for (ElementIndex ei : hasbothpoints)
if (mesh[hasbothpoints[k]].GetType() != TET) if (mesh[ei].GetType() != TET)
puretet = 0; puretet = 0;
if (!puretet) continue; if (!puretet) continue;
p1 = mesh[pi1]; Point3d p1 = mesh[pi1];
p2 = mesh[pi2]; Point3d p2 = mesh[pi2];
/* /*
pnew = Center (p1, p2); pnew = Center (p1, p2);
@ -463,13 +446,12 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
points.Elem(pi2) = p2; points.Elem(pi2) = p2;
*/ */
locfaces.SetSize (0); locfaces.SetSize (0);
for (k = 0; k < hasbothpoints.Size(); k++) for (ElementIndex ei : hasbothpoints)
{ {
const Element & el = mesh[hasbothpoints[k]]; const Element & el = mesh[ei];
for (l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
if (el[l] == pi1 || el[l] == pi2) if (el[l] == pi1 || el[l] == pi2)
{ {
INDEX_3 i3; INDEX_3 i3;
@ -486,7 +468,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
par.maxit_linsearch = 50; par.maxit_linsearch = 50;
par.maxit_bfgs = 20; par.maxit_bfgs = 20;
pnew = Center (p1, p2); Point3d pnew = Center (p1, p2);
Vector px(3); Vector px(3);
px(0) = pnew.X(); px(0) = pnew.X();
px(1) = pnew.Y(); px(1) = pnew.Y();
@ -501,11 +483,10 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
pnew.Y() = px(1); pnew.Y() = px(1);
pnew.Z() = px(2); pnew.Z() = px(2);
PointIndex hpinew = mesh.AddPoint (pnew);
int hpinew = mesh.AddPoint (pnew);
// ptyps.Append (INNERPOINT); // ptyps.Append (INNERPOINT);
for (k = 0; k < hasbothpoints.Size(); k++) for (int k = 0; k < hasbothpoints.Size(); k++)
{ {
Element & oldel = mesh[hasbothpoints[k]]; Element & oldel = mesh[hasbothpoints[k]];
Element newel1 = oldel; Element newel1 = oldel;
@ -515,7 +496,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
newel1.flags.illegal_valid = 0; newel1.flags.illegal_valid = 0;
newel2.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 (newel1[l] == pi2) newel1[l] = hpinew;
if (newel2[l] == pi1) newel2[l] = hpinew; if (newel2[l] == pi1) newel2[l] = hpinew;
@ -536,15 +517,15 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
PointIndex pinew = mesh.AddPoint (pnew); 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; Element newel = oldel;
newel.flags.illegal_valid = 0; newel.flags.illegal_valid = 0;
oldel.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]); origpoint.Clear (oldel[l]);
@ -554,7 +535,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
mesh.AddVolumeElement (newel); mesh.AddVolumeElement (newel);
} }
j = 10; j = 10; // end j-loop
} }
} }
} }
@ -572,11 +553,9 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
int cntill = 0; int cntill = 0;
ne = mesh.GetNE(); ne = mesh.GetNE();
for (ei = 0; ei < ne; ei++) for (ElementIndex ei = 0; ei < ne; ei++)
{ if (!mesh.LegalTet (mesh[ei]))
if (!mesh.LegalTet (mesh[ei])) cntill++;
cntill++;
}
// cout << cntill << " illegal tets" << endl; // cout << cntill << " illegal tets" << endl;
} }
@ -590,7 +569,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
const BitArray * working_elements) 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; int cnt = 0;
Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET); 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 // find elements on node
for (ElementIndex ei = 0; ei < ne; ei++) 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++) for (int j = 0; j < mesh[ei].GetNP(); j++)
elementsonnode.Add (mesh[ei][j], ei); elementsonnode.Add (mesh[ei][j], ei);
*/
// INDEX_2_HASHTABLE<int> edgeused(2 * ne + 5); // INDEX_2_HASHTABLE<int> edgeused(2 * ne + 5);
@ -648,7 +631,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
if (multithread.terminate) if (multithread.terminate)
break; 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; continue;
multithread.percent = 100.0 * (ei+1) / ne; multithread.percent = 100.0 * (ei+1) / ne;
@ -678,18 +661,14 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
const Element & elemi = mesh[ei]; const Element & elemi = mesh[ei];
if (elemi.IsDeleted()) continue; if (elemi.IsDeleted()) continue;
// (*testout) << "check element " << elemi << endl;
int mattyp = elemi.GetIndex(); int mattyp = elemi.GetIndex();
static const int tetedges[6][2] = static const int tetedges[6][2] =
{ { 0, 1 }, { 0, 2 }, { 0, 3 }, { { 0, 1 }, { 0, 2 }, { 0, 3 },
{ 1, 2 }, { 1, 3 }, { 2, 3 } }; { 1, 2 }, { 1, 3 }, { 2, 3 } };
pi1 = elemi[tetedges[j][0]]; PointIndex pi1 = elemi[tetedges[j][0]];
pi2 = elemi[tetedges[j][1]]; PointIndex pi2 = elemi[tetedges[j][1]];
if (pi2 < pi1) Swap (pi1, pi2); if (pi2 < pi1) Swap (pi1, pi2);
@ -718,21 +697,22 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
if (has1 && has2) if (has1 && has2)
{ // only once { // only once
for (int l = 0; l < hasbothpoints.Size(); l++) if (hasbothpoints.Contains (elnr))
if (hasbothpoints[l] == elnr) has1 = false;
has1 = 0;
if (has1) if (has1)
hasbothpoints.Append (elnr); {
hasbothpoints.Append (elnr);
}
} }
} }
bool puretet = 1; bool puretet = true;
for (int k = 0; k < hasbothpoints.Size(); k++) for (ElementIndex ei : hasbothpoints)
if (mesh[hasbothpoints[k]].GetType () != TET) if (mesh[ei].GetType () != TET)
puretet = 0; puretet = false;
if (!puretet)
continue; if (!puretet) continue;
int nsuround = hasbothpoints.Size(); int nsuround = hasbothpoints.Size();
@ -763,10 +743,10 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
for (int k = 0; k < 3; k++) // JS, 201212 for (int k = 0; k < 3; k++) // JS, 201212
{ {
const Element & elemk = mesh[hasbothpoints[k]]; const Element & elemk = mesh[hasbothpoints[k]];
bool has1 = 0; bool has1 = false;
for (int l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
if (elemk[l] == pi4) if (elemk[l] == pi4)
has1 = 1; has1 = true;
if (has1) if (has1)
{ {
for (int l = 0; l < 4; l++) 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"); throw NgException("Illegal state observed in SwapImprove");
el32[0] = pi1; el32[0] = pi1;
el32[1] = pi2; el32[1] = pi2;
@ -882,7 +861,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
mesh[hasbothpoints[0]] = el21; mesh[hasbothpoints[0]] = el21;
mesh[hasbothpoints[1]] = el22; mesh[hasbothpoints[1]] = el22;
for (int l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
mesh[hasbothpoints[2]][l] = 0; mesh[hasbothpoints[2]][l].Invalidate();
mesh[hasbothpoints[2]].Delete(); mesh[hasbothpoints[2]].Delete();
for (int k = 0; k < 2; k++) for (int k = 0; k < 2; k++)
@ -912,14 +891,11 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
el1[3] = pi4; el1[3] = pi4;
} }
pi5 = 0; pi5.Invalidate();
for (int k = 0; k < 4; k++) for (int k = 0; k < 4; k++)
{ {
const Element & elem = mesh[hasbothpoints[k]]; const Element & elem = mesh[hasbothpoints[k]];
bool has1 = 0; bool has1 = elem.PNums().Contains(pi4);
for (int l = 0; l < 4; l++)
if (elem[l] == pi4)
has1 = 1;
if (has1) if (has1)
{ {
for (int l = 0; l < 4; l++) 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++) for (int k = 0; k < 4; k++)
{ {
const Element & elem = mesh[hasbothpoints[k]]; const Element & elem = mesh[hasbothpoints[k]];
bool has1 = 0; bool has1 = elem.PNums().Contains(pi3);
for (int l = 0; l < 4; l++)
if (elem[l] == pi3)
has1 = 1;
if (has1) if (has1)
{ {
for (int l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
@ -1202,7 +1175,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
Element hel(TET); Element hel(TET);
ArrayMem<PointIndex, 50> suroundpts(nsuround); ArrayMem<PointIndex, 50> suroundpts(nsuround);
ArrayMem<char, 50> tetused(nsuround); ArrayMem<bool, 50> tetused(nsuround);
Element & elem = mesh[hasbothpoints[0]]; Element & elem = mesh[hasbothpoints[0]];
@ -1228,33 +1201,32 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
// suroundpts.SetSize (nsuround); // suroundpts.SetSize (nsuround);
suroundpts = -17;
suroundpts[0] = pi3; suroundpts[0] = pi3;
suroundpts[1] = pi4; suroundpts[1] = pi4;
tetused = 0; tetused = false;
tetused[0] = 1; tetused[0] = true;
for (int l = 2; l < nsuround; l++) for (int l = 2; l < nsuround; l++)
{ {
int oldpi = suroundpts[l-1]; PointIndex oldpi = suroundpts[l-1];
int newpi = 0; PointIndex newpi;
newpi.Invalidate();
for (int k = 0; k < nsuround && !newpi.IsValid(); k++)
for (int k = 0; k < nsuround && !newpi; k++)
if (!tetused[k]) if (!tetused[k])
{ {
const Element & nel = mesh[hasbothpoints[k]]; const Element & nel = mesh[hasbothpoints[k]];
for (int k2 = 0; k2 < 4 && !newpi.IsValid(); k2++)
for (int k2 = 0; k2 < 4 && !newpi; k2++)
if (nel[k2] == oldpi) if (nel[k2] == oldpi)
{ {
newpi = newpi =
nel[0] + nel[1] + nel[2] + nel[3] nel[0] + nel[1] + nel[2] + nel[3]
- pi1 - pi2 - oldpi; - pi1 - pi2 - oldpi;
tetused[k] = 1; tetused[k] = true;
suroundpts[l] = newpi; suroundpts[l] = newpi;
} }
} }
} }
@ -1406,8 +1378,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
*/ */
rel.Delete(); rel.Delete();
for (int k1 = 0; k1 < 4; k1++) for (int k1 = 0; k1 < 4; k1++)
rel[k1] = 0; rel[k1].Invalidate();
} }
} }
} }

27
libsrc/meshing/meshclass.cpp
View File

@ -5367,30 +5367,9 @@ namespace netgen
int Mesh :: MarkIllegalElements () int Mesh :: MarkIllegalElements ()
{ {
int cnt = 0; int cnt = 0;
int i; for (auto & el : VolumeElements())
if (!LegalTet (el))
for (i = 1; i <= GetNE(); i++) cnt++;
{
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();
}
return cnt; return cnt;
} }

58
libsrc/meshing/meshclass.hpp
View File

@ -24,8 +24,8 @@ namespace netgen
{ {
public: public:
typedef ::netgen::T_POINTS T_POINTS; typedef ::netgen::T_POINTS T_POINTS;
typedef Array<Element> T_VOLELEMENTS; typedef Array<Element, 0, ElementIndex> T_VOLELEMENTS;
typedef Array<Element2d> T_SURFELEMENTS; typedef Array<Element2d, 0, SurfaceElementIndex> T_SURFELEMENTS;
private: private:
/// point coordinates /// point coordinates
@ -192,32 +192,34 @@ namespace netgen
void ClearSurfaceElements(); void ClearSurfaceElements();
/// ///
DLL_HEADER void ClearVolumeElements() DLL_HEADER void ClearVolumeElements()
{ {
volelements.SetSize(0); volelements.SetSize(0);
timestamp = NextTimeStamp(); timestamp = NextTimeStamp();
} }
/// ///
DLL_HEADER void ClearSegments() DLL_HEADER void ClearSegments()
{ {
segments.SetSize(0); segments.SetSize(0);
timestamp = NextTimeStamp(); timestamp = NextTimeStamp();
} }
/// ///
bool TestOk () const; bool TestOk () const;
void SetAllocSize(int nnodes, int nsegs, int nsel, int nel); 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 = 1);
DLL_HEADER PointIndex AddPoint (const Point3d & p, int layer, POINTTYPE type); DLL_HEADER PointIndex AddPoint (const Point3d & p, int layer, POINTTYPE type);
int GetNP () const { return points.Size(); } int GetNP () const { return points.Size(); }
// [[deprecated("Use Point(PointIndex) instead of int !")]]
MeshPoint & Point(int i) { return points.Elem(i); } MeshPoint & Point(int i) { return points.Elem(i); }
MeshPoint & Point(PointIndex pi) { return points[pi]; } 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(int i) const { return points.Get(i); }
const MeshPoint & Point(PointIndex pi) const { return points[pi]; } const MeshPoint & Point(PointIndex pi) const { return points[pi]; }
@ -231,16 +233,20 @@ namespace netgen
DLL_HEADER SegmentIndex AddSegment (const Segment & s); DLL_HEADER SegmentIndex AddSegment (const Segment & s);
void DeleteSegment (int segnr) void DeleteSegment (int segnr)
{ {
segments.Elem(segnr)[0] = PointIndex::BASE-1; segments.Elem(segnr)[0].Invalidate();
segments.Elem(segnr)[1] = PointIndex::BASE-1; segments.Elem(segnr)[1].Invalidate();
} }
/*
void FullDeleteSegment (int segnr) // von wem ist das ??? void FullDeleteSegment (int segnr) // von wem ist das ???
{ {
segments.Delete(segnr-PointIndex::BASE); segments.Delete(segnr-PointIndex::BASE);
} }
*/
int GetNSeg () const { return segments.Size(); } int GetNSeg () const { return segments.Size(); }
[[deprecated("Use LineSegment(SegmentIndex) instead of int !")]]
Segment & LineSegment(int i) { return segments.Elem(i); } 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); } const Segment & LineSegment(int i) const { return segments.Get(i); }
Segment & LineSegment(SegmentIndex si) { return segments[si]; } Segment & LineSegment(SegmentIndex si) { return segments[si]; }
@ -254,22 +260,29 @@ namespace netgen
Array<Element0d> pointelements; // only via python interface Array<Element0d> pointelements; // only via python interface
DLL_HEADER SurfaceElementIndex AddSurfaceElement (const Element2d & el); DLL_HEADER SurfaceElementIndex AddSurfaceElement (const Element2d & el);
[[deprecated("Use DeleteSurfaceElement(SurfaceElementIndex) instead of int !")]]
void DeleteSurfaceElement (int eli) void DeleteSurfaceElement (int eli)
{ {
surfelements.Elem(eli).Delete(); surfelements.Elem(eli).Delete();
surfelements.Elem(eli).PNum(1) = -1; surfelements.Elem(eli).PNum(1).Invalidate();
surfelements.Elem(eli).PNum(2) = -1; surfelements.Elem(eli).PNum(2).Invalidate();
surfelements.Elem(eli).PNum(3) = -1; surfelements.Elem(eli).PNum(3).Invalidate();
timestamp = NextTimeStamp(); timestamp = NextTimeStamp();
} }
void DeleteSurfaceElement (SurfaceElementIndex eli) 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(); } int GetNSE () const { return surfelements.Size(); }
[[deprecated("Use SurfaceElement(SurfaceElementIndex) instead of int !")]]
Element2d & SurfaceElement(int i) { return surfelements.Elem(i); } 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); } const Element2d & SurfaceElement(int i) const { return surfelements.Get(i); }
Element2d & SurfaceElement(SurfaceElementIndex i) { return surfelements[i]; } Element2d & SurfaceElement(SurfaceElementIndex i) { return surfelements[i]; }
const Element2d & SurfaceElement(SurfaceElementIndex i) const { return surfelements[i]; } const Element2d & SurfaceElement(SurfaceElementIndex i) const { return surfelements[i]; }
@ -290,7 +303,9 @@ namespace netgen
int GetNE () const { return volelements.Size(); } int GetNE () const { return volelements.Size(); }
// [[deprecated("Use VolumeElement(ElementIndex) instead of int !")]]
Element & VolumeElement(int i) { return volelements.Elem(i); } 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); } const Element & VolumeElement(int i) const { return volelements.Get(i); }
Element & VolumeElement(ElementIndex i) { return volelements[i]; } Element & VolumeElement(ElementIndex i) { return volelements[i]; }
const Element & VolumeElement(ElementIndex i) const { return volelements[i]; } const Element & VolumeElement(ElementIndex i) const { return volelements[i]; }
@ -305,9 +320,8 @@ namespace netgen
ELEMENTTYPE ElementType (ElementIndex i) const ELEMENTTYPE ElementType (ElementIndex i) const
{ return (volelements[i].flags.fixed) ? FIXEDELEMENT : FREEELEMENT; } { return (volelements[i].flags.fixed) ? FIXEDELEMENT : FREEELEMENT; }
const T_VOLELEMENTS & VolumeElements() const { return volelements; } const auto & VolumeElements() const { return volelements; }
T_VOLELEMENTS & VolumeElements() { return volelements; } auto & VolumeElements() { return volelements; }
/// ///
DLL_HEADER double ElementError (int eli, const MeshingParameters & mp) const; DLL_HEADER double ElementError (int eli, const MeshingParameters & mp) const;
@ -317,17 +331,13 @@ namespace netgen
/// ///
void ClearLockedPoints (); void ClearLockedPoints ();
const Array<PointIndex> & LockedPoints() const const auto & LockedPoints() const { return lockedpoints; }
{ return lockedpoints; }
/// Returns number of domains /// Returns number of domains
int GetNDomains() const; int GetNDomains() const;
/// ///
int GetDimension() const int GetDimension() const { return dimension; }
{ return dimension; } void SetDimension (int dim) { dimension = dim; }
void SetDimension(int dim)
{ dimension = dim; }
/// sets internal tables /// sets internal tables
void CalcSurfacesOfNode (); void CalcSurfacesOfNode ();
@ -573,10 +583,10 @@ namespace netgen
/// ///
int AddFaceDescriptor(const FaceDescriptor& fd) int AddFaceDescriptor(const FaceDescriptor& fd)
{ return facedecoding.Append(fd); } { facedecoding.Append(fd); return facedecoding.Size(); }
int AddEdgeDescriptor(const EdgeDescriptor & fd) 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); DLL_HEADER void SetMaterial (int domnr, const string & mat);

2
libsrc/meshing/meshfunc.cpp
View File

@ -708,7 +708,7 @@ namespace netgen
break; break;
PrintMessage (5, nillegal, " illegal tets"); PrintMessage (5, nillegal, " illegal tets");
optmesh.SplitImprove (mesh3d, OPT_LEGAL); optmesh.SplitImprove (mesh3d, OPT_LEGAL);
mesh3d.MarkIllegalElements(); // test mesh3d.MarkIllegalElements(); // test
optmesh.SwapImprove (mesh3d, OPT_LEGAL); optmesh.SwapImprove (mesh3d, OPT_LEGAL);

95
libsrc/meshing/meshing3.cpp
View File

@ -176,14 +176,14 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
NgProfiler::RegionTimer reg (meshing3_timer); NgProfiler::RegionTimer reg (meshing3_timer);
Array<Point3d > locpoints; // local points Array<Point3d, PointIndex::BASE> locpoints; // local points
Array<MiniElement2d> locfaces; // local faces Array<MiniElement2d> locfaces; // local faces
Array<PointIndex> pindex; // mapping from local to front point numbering Array<PointIndex, PointIndex::BASE> pindex; // mapping from local to front point numbering
Array<int> allowpoint; // point is allowd ? Array<int, PointIndex::BASE> allowpoint; // point is allowd ?
Array<INDEX> findex; // mapping from local to front face numbering Array<INDEX> findex; // mapping from local to front face numbering
//INDEX_2_HASHTABLE<int> connectedpairs(100); // connecgted pairs for prism meshing //INDEX_2_HASHTABLE<int> connectedpairs(100); // connecgted pairs for prism meshing
Array<Point3d > plainpoints; // points in reference coordinates Array<Point3d, PointIndex::BASE> plainpoints; // points in reference coordinates
Array<int> delpoints, delfaces; // points and lines to be deleted Array<int> delpoints, delfaces; // points and lines to be deleted
Array<Element> locelements; // new generated elements Array<Element> locelements; // new generated elements
@ -206,9 +206,9 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
// for star-shaped domain meshing // for star-shaped domain meshing
Array<MeshPoint> grouppoints; Array<MeshPoint, PointIndex::BASE> grouppoints;
Array<MiniElement2d> groupfaces; Array<MiniElement2d> groupfaces;
Array<PointIndex> grouppindex; Array<PointIndex, PointIndex::BASE> grouppindex;
Array<INDEX> groupfindex; Array<INDEX> groupfindex;
@ -269,9 +269,9 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
} }
const MiniElement2d & bel = adfront->GetFace (baseelem); const MiniElement2d & bel = adfront->GetFace (baseelem);
const Point3d & p1 = adfront->GetPoint (bel.PNum(1)); const Point3d & p1 = adfront->GetPoint (bel[0]);
const Point3d & p2 = adfront->GetPoint (bel.PNum(2)); const Point3d & p2 = adfront->GetPoint (bel[1]);
const Point3d & p3 = adfront->GetPoint (bel.PNum(3)); const Point3d & p3 = adfront->GetPoint (bel[2]);
// (*testout) << endl << "base = " << bel << endl; // (*testout) << endl << "base = " << bel << endl;
@ -302,9 +302,6 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
// (*testout) << "locfaces = " << endl << locfaces << endl; // (*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; //loktestmode = 1;
testmode = loktestmode; //changed testmode = loktestmode; //changed
@ -316,6 +313,9 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
if (loktestmode) if (loktestmode)
{ {
(*testout) << "baseel = " << baseelem << ", ind = " << findex.Get(1) << endl; (*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; (*testout) << "pi = " << pi1 << ", " << pi2 << ", " << pi3 << endl;
} }
@ -369,7 +369,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
grouppindex, groupfindex); grouppindex, groupfindex);
bool onlytri = 1; bool onlytri = 1;
for (i = 0; i < groupfaces.Size(); i++) for (auto i : groupfaces.Range())
if (groupfaces[i].GetNP() != 3) if (groupfaces[i].GetNP() != 3)
onlytri = 0; onlytri = 0;
@ -441,31 +441,28 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
{ {
// set transformatino to reference coordinates // 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)), trans.Set (locpoints[locfaces[0].PNumMod(1+rotind)],
locpoints.Get(locfaces.Get(1).PNumMod(2+rotind)), locpoints[locfaces[0].PNumMod(2+rotind)],
locpoints.Get(locfaces.Get(1).PNumMod(3+rotind)), hshould); locpoints[locfaces[0].PNumMod(3+rotind)], hshould);
} }
else else
{ {
trans.Set (locpoints.Get(locfaces.Get(1).PNumMod(1+rotind)), trans.Set (locpoints[locfaces[0].PNumMod(1+rotind)],
locpoints.Get(locfaces.Get(1).PNumMod(2+rotind)), locpoints[locfaces[0].PNumMod(2+rotind)],
locpoints.Get(locfaces.Get(1).PNumMod(4+rotind)), hshould); locpoints[locfaces[0].PNumMod(4+rotind)], hshould);
} }
trans.ToPlain (locpoints, plainpoints); // trans.ToPlain (locpoints, plainpoints);
plainpoints.SetSize (locpoints.Size());
for (i = 1; i <= allowpoint.Size(); i++) for (auto i : locpoints.Range())
{ trans.ToPlain (locpoints[i], plainpoints[i]);
if (plainpoints.Get(i).Z() > 0)
{ for (auto i : allowpoint.Range())
//if(loktestmode) if (plainpoints[i].Z() > 0)
// (*testout) << "plainpoints.Get(i).Z() = " << plainpoints.Get(i).Z() << " > 0" << endl; allowpoint[i] = false;
allowpoint.Elem(i) = 0;
}
}
stat.cnttrials++; stat.cnttrials++;
@ -508,7 +505,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
if (found) stat.cntsucc++; if (found) stat.cntsucc++;
locpoints.SetSize (plainpoints.Size()); 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)); 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 // avoid meshing from large to small mesh-size
if (uselocalh && found && stat.qualclass <= 3) 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; Point3d pmax = pmin;
for (j = 2; j <= 4; j++) 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); pmin.SetToMin (hp);
pmax.SetToMax (hp); pmax.SetToMax (hp);
} }
@ -533,10 +530,10 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
} }
if (found) if (found)
{ {
for (i = 1; i <= locelements.Size(); i++) for (int i = 1; i <= locelements.Size(); i++)
for (j = 1; j <= 4; j++) 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) if (Dist (hp, pmid) > hinner)
found = 0; found = 0;
} }
@ -655,25 +652,25 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
pindex.SetSize(locpoints.Size()); 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)); PointIndex globind = mesh.AddPoint (locpoints.Get(i));
pindex.Elem(i) = adfront -> AddPoint (locpoints.Get(i), globind); 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; Point3d * hp1, * hp2, * hp3, * hp4;
hp1 = &locpoints.Elem(locelements.Get(i).PNum(1)); hp1 = &locpoints[locelements.Get(i).PNum(1)];
hp2 = &locpoints.Elem(locelements.Get(i).PNum(2)); hp2 = &locpoints[locelements.Get(i).PNum(2)];
hp3 = &locpoints.Elem(locelements.Get(i).PNum(3)); hp3 = &locpoints[locelements.Get(i).PNum(3)];
hp4 = &locpoints.Elem(locelements.Get(i).PNum(4)); hp4 = &locpoints[locelements.Get(i).PNum(4)];
tetvol += (1.0 / 6.0) * ( Cross ( *hp2 - *hp1, *hp3 - *hp1) * (*hp4 - *hp1) ); tetvol += (1.0 / 6.0) * ( Cross ( *hp2 - *hp1, *hp3 - *hp1) * (*hp4 - *hp1) );
for (j = 1; j <= locelements.Get(i).NP(); j++) for (j = 1; j <= locelements.Get(i).NP(); j++)
locelements.Elem(i).PNum(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)); mesh.AddVolumeElement (locelements.Get(i));
stat.cntelem++; stat.cntelem++;
@ -683,7 +680,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
{ {
for (j = 1; j <= locfaces.Get(i).GetNP(); j++) for (j = 1; j <= locfaces.Get(i).GetNP(); j++)
locfaces.Elem(i).PNum(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; // (*testout) << "add face " << locfaces.Get(i) << endl;
adfront->AddFace (locfaces.Get(i)); adfront->AddFace (locfaces.Get(i));
} }

3
libsrc/meshing/meshing3.hpp
View File

@ -42,7 +42,8 @@ public:
MESHING3_RESULT GenerateMesh (Mesh & mesh, const MeshingParameters & mp); MESHING3_RESULT GenerateMesh (Mesh & mesh, const MeshingParameters & mp);
/// ///
int ApplyRules (Array<Point3d> & lpoints, Array<int> & allowpoint, int ApplyRules (Array<Point3d, PointIndex::BASE> & lpoints,
Array<int, PointIndex::BASE> & allowpoint,
Array<MiniElement2d> & lfaces, INDEX lfacesplit, Array<MiniElement2d> & lfaces, INDEX lfacesplit,
INDEX_2_HASHTABLE<int> & connectedpairs, INDEX_2_HASHTABLE<int> & connectedpairs,
Array<Element> & elements, Array<Element> & elements,

11
libsrc/meshing/meshtype.cpp
View File

@ -538,6 +538,9 @@ namespace netgen
shape(3) = (1-p(0))* p(1) ; shape(3) = (1-p(0))* p(1) ;
break; break;
} }
default:
throw NgException ("illegal element type in GetShapeNew");
} }
} }
@ -562,6 +565,8 @@ namespace netgen
shape(3) = (1-p(0))* p(1) ; shape(3) = (1-p(0))* p(1) ;
break; break;
} }
default:
throw NgException ("illegal element type in GetShapeNew");
} }
} }
@ -641,6 +646,8 @@ namespace netgen
dshape(3,1) = (1-p(0)); dshape(3,1) = (1-p(0));
break; break;
} }
default:
throw NgException ("illegal element type in GetDShapeNew");
} }
} }
@ -1739,12 +1746,14 @@ namespace netgen
{ 0, 0, 1, 1 }, { 0, 0, 1, 1 },
{ 0, 0, 0, 1 }, { 0, 0, 0, 1 },
}; };
double * pp = NULL; double * pp = NULL;
switch (typ) switch (typ)
{ {
case TET: pp = &eltetqp[0][0]; break; case TET: pp = &eltetqp[0][0]; break;
case TET10: pp = &eltet10qp[ip-1][0]; break; case TET10: pp = &eltet10qp[ip-1][0]; break;
default:
throw NgException ("illegal element shape in GetIntegrationPoint");
} }
p(0) = pp[0]; p(0) = pp[0];

28
libsrc/meshing/meshtype.hpp
View File

@ -125,6 +125,8 @@ namespace netgen
PointIndex operator-- (int) { PointIndex hi(*this); i--; return hi; } PointIndex operator-- (int) { PointIndex hi(*this); i--; return hi; }
PointIndex operator++ () { i++; return *this; } PointIndex operator++ () { i++; return *this; }
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 #ifdef BASE0
enum { BASE = 0 }; enum { BASE = 0 };
#else #else
@ -134,7 +136,7 @@ namespace netgen
inline istream & operator>> (istream & ist, PointIndex & pi) 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) 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= (const ElementIndex & ai) { i = ai.i; return *this; }
ElementIndex & operator= (int ai) { i = ai; return *this; } ElementIndex & operator= (int ai) { i = ai; return *this; }
operator int () const { return i; } operator int () const { return i; }
ElementIndex & operator++ (int) { i++; return *this; } ElementIndex operator++ (int) { return ElementIndex(i++); }
ElementIndex & operator-- (int) { i--; return *this; } ElementIndex operator-- (int) { return ElementIndex(i--); }
ElementIndex & operator++ () { ++i; return *this; }
ElementIndex & operator-- () { --i; return *this; }
}; };
inline istream & operator>> (istream & ist, ElementIndex & pi) inline istream & operator>> (istream & ist, ElementIndex & pi)
@ -179,9 +183,11 @@ namespace netgen
{ i = ai.i; return *this; } { i = ai.i; return *this; }
SurfaceElementIndex & operator= (int ai) { i = ai; return *this; } SurfaceElementIndex & operator= (int ai) { i = ai; return *this; }
operator int () const { return i; } 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 inc) { i+=inc; return *this; }
SurfaceElementIndex & operator-- (int) { i--; return *this; }
}; };
inline istream & operator>> (istream & ist, SurfaceElementIndex & pi) inline istream & operator>> (istream & ist, SurfaceElementIndex & pi)
@ -344,6 +350,7 @@ namespace netgen
default: default:
PrintSysError ("Element2d::SetType, illegal type ", int(typ)); PrintSysError ("Element2d::SetType, illegal type ", int(typ));
} }
is_curved = (np >= 4);
} }
/// ///
int GetNP() const { return np; } int GetNP() const { return np; }
@ -386,6 +393,8 @@ namespace netgen
FlatArray<const PointIndex> PNums () const FlatArray<const PointIndex> PNums () const
{ return FlatArray<const PointIndex> (np, &pnum[0]); } { return FlatArray<const PointIndex> (np, &pnum[0]); }
FlatArray<PointIndex> PNums ()
{ return FlatArray<PointIndex> (np, &pnum[0]); }
/// ///
PointIndex & PNum (int i) { return pnum[i-1]; } PointIndex & PNum (int i) { return pnum[i-1]; }
@ -470,8 +479,13 @@ namespace netgen
int pi, Vec2d & dir, double & dd) const; 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 bool IsDeleted () const
{ {
#ifdef DEBUG #ifdef DEBUG

8
libsrc/meshing/python_mesh.cpp
View File

@ -304,8 +304,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
ExportArray<Element>(m); ExportArray<Element,0,ElementIndex>(m);
ExportArray<Element2d>(m); ExportArray<Element2d,0,SurfaceElementIndex>(m);
ExportArray<Segment>(m); ExportArray<Segment>(m);
ExportArray<Element0d>(m); ExportArray<Element0d>(m);
ExportArray<MeshPoint,PointIndex::BASE,PointIndex>(m); ExportArray<MeshPoint,PointIndex::BASE,PointIndex>(m);
@ -384,11 +384,11 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
.def_property("dim", &Mesh::GetDimension, &Mesh::SetDimension) .def_property("dim", &Mesh::GetDimension, &Mesh::SetDimension)
.def("Elements3D", .def("Elements3D",
static_cast<Array<Element>&(Mesh::*)()> (&Mesh::VolumeElements), static_cast<Array<Element,0,ElementIndex>&(Mesh::*)()> (&Mesh::VolumeElements),
py::return_value_policy::reference) py::return_value_policy::reference)
.def("Elements2D", .def("Elements2D",
static_cast<Array<Element2d>&(Mesh::*)()> (&Mesh::SurfaceElements), static_cast<Array<Element2d,0,SurfaceElementIndex>&(Mesh::*)()> (&Mesh::SurfaceElements),
py::return_value_policy::reference) py::return_value_policy::reference)
.def("Elements1D", .def("Elements1D",

14
libsrc/meshing/refine.cpp
View File

@ -383,11 +383,11 @@ namespace netgen
swap (pnums.Elem(3), pnums.Elem(4)); swap (pnums.Elem(3), pnums.Elem(4));
for (int j = 0; j < 6; j++) for (int j = 0; j < 6; j++)
{ {
INDEX_2 i2; PointIndex pi1 = pnums.Get(betw[j][0]);
i2.I1() = pnums.Get(betw[j][0]); PointIndex pi2 = pnums.Get(betw[j][1]);
i2.I2() = pnums.Get(betw[j][1]); INDEX_2 i2 (pi1, pi2);
i2.Sort(); i2.Sort();
/* /*
if (between.Used(i2)) if (between.Used(i2))
@ -405,8 +405,8 @@ namespace netgen
if (!pointset[pinew]) if (!pointset[pinew])
{ {
pointset[pinew] = true; pointset[pinew] = true;
mesh.Point(pinew) = Center(mesh.Point(i2.I1()), mesh.Point(pinew) = Center(mesh.Point(pi1),
mesh.Point(i2.I2())); mesh.Point(pi2));
} }
} }

5
libsrc/meshing/ruler2.cpp
View File

@ -593,8 +593,9 @@ namespace netgen
Point2d np = rule->GetPoint(i); Point2d np = rule->GetPoint(i);
np.X() += newu (2 * (i-oldnp) - 2); np.X() += newu (2 * (i-oldnp) - 2);
np.Y() += newu (2 * (i-oldnp) - 1); np.Y() += newu (2 * (i-oldnp) - 1);
pmap.Elem(i) = lpoints.Append (np); lpoints.Append (np);
pmap.Elem(i) = lpoints.Size();
} }
} }

358
libsrc/meshing/ruler3.cpp
View File

@ -8,29 +8,29 @@ extern double minother;
extern double minwithoutother; extern double minwithoutother;
static double CalcElementBadness (const Array<Point3d> & points, static double CalcElementBadness (const Array<Point3d, PointIndex::BASE> & points,
const Element & elem) const Element & elem)
{ {
double vol, l, l4, l5, l6; double vol, l, l4, l5, l6;
if (elem.GetNP() != 4) if (elem.GetNP() != 4)
{ {
if (elem.GetNP() == 5) 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; if (z > -1e-8) return 1e8;
return (-1 / z) - z; // - 2; return (-1 / z) - z; // - 2;
} }
return 0; return 0;
} }
Vec3d v1 = points.Get(elem.PNum(2)) - points.Get(elem.PNum(1)); Vec3d v1 = points[elem.PNum(2)] - points[elem.PNum(1)];
Vec3d v2 = points.Get(elem.PNum(3)) - points.Get(elem.PNum(1)); Vec3d v2 = points[elem.PNum(3)] - points[elem.PNum(1)];
Vec3d v3 = points.Get(elem.PNum(4)) - points.Get(elem.PNum(1)); Vec3d v3 = points[elem.PNum(4)] - points[elem.PNum(1)];
vol = - (Cross (v1, v2) * v3); vol = - (Cross (v1, v2) * v3);
l4 = Dist (points.Get(elem.PNum(2)), points.Get(elem.PNum(3))); l4 = Dist (points[elem.PNum(2)], points[elem.PNum(3)]);
l5 = Dist (points.Get(elem.PNum(2)), points.Get(elem.PNum(4))); l5 = Dist (points[elem.PNum(2)], points[elem.PNum(4)]);
l6 = Dist (points.Get(elem.PNum(3)), points.Get(elem.PNum(4))); l6 = Dist (points[elem.PNum(3)], points[elem.PNum(4)]);
l = v1.Length() + v2.Length() + v3.Length() + l4 + l5 + l6; l = v1.Length() + v2.Length() + v3.Length() + l4 + l5 + l6;
@ -49,8 +49,8 @@ static double CalcElementBadness (const Array<Point3d> & points,
int Meshing3 :: ApplyRules int Meshing3 :: ApplyRules
( (
Array<Point3d> & lpoints, // in: local points, out: old+new local points Array<Point3d, PointIndex::BASE> & lpoints, // in: local points, out: old+new local points
Array<int> & allowpoint, // in: 2 .. it is allowed to use pointi, 1..will be allowed later, 0..no means Array<int, PointIndex::BASE> & allowpoint, // in: 2 .. it is allowed to use pointi, 1..will be allowed later, 0..no means
Array<MiniElement2d> & lfaces, // in: local faces, out: old+new local faces Array<MiniElement2d> & lfaces, // in: local faces, out: old+new local faces
INDEX lfacesplit, // for local faces in outer radius INDEX lfacesplit, // for local faces in outer radius
INDEX_2_HASHTABLE<int> & connectedpairs, // connected pairs for prism-meshing INDEX_2_HASHTABLE<int> & connectedpairs, // connected pairs for prism-meshing
@ -65,25 +65,23 @@ int Meshing3 :: ApplyRules
{ {
NgProfiler::RegionTimer regtot(97); NgProfiler::RegionTimer regtot(97);
int i, j, k, ri, nfok, npok, incnpok, refpi, locpi, locfi, locfr; float err, minerr, teterr, minteterr;
float hf, err, minerr, teterr, minteterr;
char ok, found, hc; char ok, found, hc;
vnetrule * rule; // vnetrule * rule;
Vector oldu, newu, newu1, newu2, allp; Vector oldu, newu, newu1, newu2, allp;
Vec3d ui; Vec3d ui;
Point3d np; Point3d np;
int oldnp, noldlp, noldlf;
const MiniElement2d * locface = NULL; const MiniElement2d * locface = NULL;
int loktestmode; int loktestmode;
Array<int> pused; // point is already mapped Array<int, PointIndex::BASE> pused; // point is already mapped, number of uses
Array<char> fused; // face is already mapped Array<char> fused; // face is already mapped
Array<int> pmap; // map of reference point to local point Array<PointIndex> pmap; // map of reference point to local point
Array<char> pfixed; // point mapped by face-map Array<bool, PointIndex::BASE> pfixed; // point mapped by face-map
Array<int> fmapi; // face in reference is mapped to face nr ... Array<int> fmapi; // face in reference is mapped to face nr ...
Array<int> fmapr; // face in reference is rotated to map Array<int> fmapr; // face in reference is rotated to map
Array<Point3d> transfreezone; // transformed free-zone Array<Point3d> transfreezone; // transformed free-zone
INDEX_2_CLOSED_HASHTABLE<int> ledges(100); // edges in local environment INDEX_2_CLOSED_HASHTABLE<int> ledges(100); // edges in local environment
Array<Point3d> tempnewpoints; Array<Point3d> tempnewpoints;
@ -108,9 +106,10 @@ int Meshing3 :: ApplyRules
fnearness.SetSize (lfacesplit); fnearness.SetSize (lfacesplit);
pnearness = INT_MAX/10; 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::RegionTimer reg2(98);
NgProfiler::StartTimer (90); NgProfiler::StartTimer (90);
@ -118,24 +117,24 @@ int Meshing3 :: ApplyRules
for (int loop = 0; loop < 2; loop++) 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]; const MiniElement2d & hface = lfaces[i];
int minn = INT_MAX-1; 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 (hi < minn) minn = hi;
} }
if (minn < INT_MAX/10) if (minn < INT_MAX/10)
for (j = 0; j < hface.GetNP(); j++) for (PointIndex pi : hface.PNums())
if (pnearness[hface[j]] > minn+1) if (pnearness[pi] > minn+1)
pnearness[hface[j]] = minn+1; pnearness[pi] = minn+1;
} }
for (i = 1; i <= connectedpairs.GetNBags(); i++) for (int i = 1; i <= connectedpairs.GetNBags(); i++)
for (j = 1; j <= connectedpairs.GetBagSize(i); j++) for (int j = 1; j <= connectedpairs.GetBagSize(i); j++)
{ {
INDEX_2 edge; INDEX_2 edge;
int val; int val;
@ -147,14 +146,13 @@ int Meshing3 :: ApplyRules
if (pnearness[edge.I2()] > pnearness[edge.I1()] + 1) if (pnearness[edge.I2()] > pnearness[edge.I1()] + 1)
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; int sum = 0;
for (j = 0; j < lfaces[i].GetNP(); j++) for (PointIndex pi : lfaces[i].PNums())
sum += pnearness[lfaces[i][j]]; sum += pnearness[pi];
fnearness[i] = sum; fnearness[i] = sum;
} }
@ -165,34 +163,35 @@ int Meshing3 :: ApplyRules
// find bounding boxes of faces // find bounding boxes of faces
triboxes.SetSize (lfaces.Size()); 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]; const MiniElement2d & face = lfaces[i];
triboxes[i].SetPoint (lpoints.Get(face[0])); triboxes[i].SetPoint (lpoints[face[0]]);
for (j = 1; j < face.GetNP(); j++) for (int j = 1; j < face.GetNP(); j++)
triboxes[i].AddPoint (lpoints.Get(face[j])); triboxes[i].AddPoint (lpoints[face[j]]);
} }
NgProfiler::StopTimer (91); NgProfiler::StopTimer (91);
NgProfiler::StartTimer (92); NgProfiler::StartTimer (92);
bool useedges = 0; bool useedges = false;
for (ri = 0; ri < rules.Size(); ri++) for (int ri = 0; ri < rules.Size(); ri++)
if (rules[ri]->GetNEd()) useedges = 1; if (rules[ri]->GetNEd()) useedges = true;
if (useedges) if (useedges)
{ {
ledges.SetSize (5 * lfacesplit); ledges.SetSize (5 * lfacesplit);
for (j = 0; j < lfacesplit; j++) for (int j = 0; j < lfacesplit; j++)
// if (fnearness[j] <= 5) // if (fnearness[j] <= 5)
{ {
const MiniElement2d & face = lfaces[j]; const MiniElement2d & face = lfaces[j];
int newp, oldp; int newp, oldp;
newp = face[face.GetNP()-1]; newp = face[face.GetNP()-1];
for (k = 0; k < face.GetNP(); k++) for (int k = 0; k < face.GetNP(); k++)
{ {
oldp = newp; oldp = newp;
newp = face[k]; newp = face[k];
@ -223,8 +222,8 @@ int Meshing3 :: ApplyRules
// check each rule: // 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; int base = (lfaces[0].GetNP() == 3) ? 100 : 200;
NgProfiler::RegionTimer regx1(base); NgProfiler::RegionTimer regx1(base);
NgProfiler::RegionTimer regx(base+ri); NgProfiler::RegionTimer regx(base+ri);
@ -232,7 +231,7 @@ int Meshing3 :: ApplyRules
// sprintf (problems.Elem(ri), ""); // sprintf (problems.Elem(ri), "");
*problems.Elem(ri) = '\0'; *problems.Elem(ri) = '\0';
rule = rules.Get(ri); vnetrule * rule = rules.Get(ri);
if (rule->GetNP(1) != lfaces[0].GetNP()) if (rule->GetNP(1) != lfaces[0].GetNP())
continue; continue;
@ -260,21 +259,21 @@ int Meshing3 :: ApplyRules
fused = 0; fused = 0;
pused = 0; pused = 0;
pmap = 0; for (auto & p : pmap) p.Invalidate();
fmapi = 0; 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; fused[0] = 1;
fmapi[0] = 1; fmapi[0] = 1;
fmapr[0] = rotind1; fmapr[0] = rotind1;
for (int j = 1; j <= lfaces[0].GetNP(); j++)
for (j = 1; j <= lfaces.Get(1).GetNP(); j++)
{ {
locpi = lfaces[0].PNumMod (j+rotind1); PointIndex locpi = lfaces[0].PNumMod (j+rotind1);
pmap.Set (rule->GetPointNr (1, j), locpi); 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 .. first nfok-1 faces are mapped properly
*/ */
nfok = 2; int nfok = 2;
NgProfiler::RegionTimer regfa(300); NgProfiler::RegionTimer regfa(300);
NgProfiler::RegionTimer regx2(base+50+ri); NgProfiler::RegionTimer regx2(base+50+ri);
while (nfok >= 2) while (nfok >= 2)
@ -293,8 +292,8 @@ int Meshing3 :: ApplyRules
// not all faces mapped // not all faces mapped
ok = 0; ok = 0;
locfi = fmapi.Get(nfok); int locfi = fmapi.Get(nfok);
locfr = fmapr.Get(nfok); int locfr = fmapr.Get(nfok);
int actfnp = rule->GetNP(nfok); int actfnp = rule->GetNP(nfok);
@ -326,28 +325,27 @@ int Meshing3 :: ApplyRules
// reference point already mapped differently ? // 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)); PointIndex locpi = pmap.Get(rule->GetPointNr (nfok, j));
if (locpi.IsValid() && locpi != locface->PNumMod(j+locfr))
if (locpi && locpi != locface->PNumMod(j+locfr))
ok = 0; ok = 0;
} }
// local point already used or point outside tolerance ? // 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; ok = 0;
else else
{ {
const Point3d & lp = lpoints.Get(locpi); const Point3d & lp = lpoints[locpi];
const Point3d & rp = rule->GetPoint(refpi); const Point3d & rp = rule->GetPoint(refpi);
if ( Dist2 (lp, rp) * rule->PointDistFactor(refpi) > minerr) if ( Dist2 (lp, rp) * rule->PointDistFactor(refpi) > minerr)
@ -370,16 +368,16 @@ int Meshing3 :: ApplyRules
fmapr.Set (nfok, locfr); fmapr.Set (nfok, locfr);
fused.Set (locfi, 1); 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 && if (rule->GetPointNr (nfok, j) <= 3 &&
pmap.Get(rule->GetPointNr(nfok, j)) != locpi) pmap.Get(rule->GetPointNr(nfok, j)) != locpi)
(*testout) << "change face1 point, mark1" << endl; (*testout) << "change face1 point, mark1" << endl;
pmap.Set(rule->GetPointNr (nfok, j), locpi); pmap.Set(rule->GetPointNr (nfok, j), locpi);
pused.Elem(locpi)++; pused[locpi]++;
} }
nfok++; nfok++;
@ -392,14 +390,15 @@ int Meshing3 :: ApplyRules
nfok--; nfok--;
fused.Set (fmapi.Get(nfok), 0); 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); int refpi = rule->GetPointNr (nfok, j);
pused.Elem(pmap.Get(refpi))--; pused[pmap.Get(refpi)]--;
if (pused.Get(pmap.Get(refpi)) == 0) 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; (*testout) << "Faces Ok" << endl;
sprintf (problems.Elem(ri), "Faces Ok"); sprintf (problems.Elem(ri), "Faces Ok");
} }
npok = 1; int npok = 1;
incnpok = 1; int incnpok = 1;
pfixed.SetSize (pmap.Size()); 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) ); pfixed.Set(i, (pmap.Get(i) != 0) );
*/
for (int i : pmap.Range())
pfixed[i] = pmap[i].IsValid();
while (npok >= 1) while (npok >= 1)
{ {
@ -444,31 +447,31 @@ int Meshing3 :: ApplyRules
else else
{ {
locpi = pmap.Elem(npok); PointIndex locpi = pmap.Elem(npok);
ok = 0; ok = 0;
if (locpi) if (locpi.IsValid())
pused.Elem(locpi)--; pused[locpi]--;
while (!ok && locpi < lpoints.Size()) while (!ok && locpi < lpoints.Size()-1+PointIndex::BASE)
{ {
ok = 1; ok = 1;
locpi++; locpi++;
if (pused.Get(locpi) || if (pused[locpi] ||
pnearness.Get(locpi) > rule->GetPNearness(npok)) pnearness[locpi] > rule->GetPNearness(npok))
{ {
ok = 0; ok = 0;
} }
else if (allowpoint.Get(locpi) != 2) else if (allowpoint[locpi] != 2)
{ {
ok = 0; ok = 0;
if (allowpoint.Get(locpi) == 1) if (allowpoint[locpi] == 1)
impossible = 0; impossible = 0;
} }
else else
{ {
const Point3d & lp = lpoints.Get(locpi); const Point3d & lp = lpoints[locpi];
const Point3d & rp = rule->GetPoint(npok); const Point3d & rp = rule->GetPoint(npok);
if ( Dist2 (lp, rp) * rule->PointDistFactor(npok) > minerr) if ( Dist2 (lp, rp) * rule->PointDistFactor(npok) > minerr)
@ -487,7 +490,7 @@ int Meshing3 :: ApplyRules
if (npok <= 3) if (npok <= 3)
(*testout) << "set face1 point, mark3" << endl; (*testout) << "set face1 point, mark3" << endl;
pused.Elem(locpi)++; pused[locpi]++;
npok++; npok++;
incnpok = 1; incnpok = 1;
} }
@ -495,7 +498,8 @@ int Meshing3 :: ApplyRules
else else
{ {
pmap.Set (npok, 0); // pmap.Set (npok, 0);
pmap.Elem(npok).Invalidate();
if (npok <= 3) if (npok <= 3)
(*testout) << "set face1 point, mark4" << endl; (*testout) << "set face1 point, mark4" << endl;
@ -516,8 +520,8 @@ int Meshing3 :: ApplyRules
if (loktestmode) if (loktestmode)
{ {
(*testout) << "Mapping found!!: Rule " << rule->Name() << endl; (*testout) << "Mapping found!!: Rule " << rule->Name() << endl;
for (i = 1; i <= pmap.Size(); i++) for (auto pi : pmap)
(*testout) << pmap.Get(i) << " "; (*testout) << pi << " ";
(*testout) << endl; (*testout) << endl;
sprintf (problems.Elem(ri), "mapping found"); sprintf (problems.Elem(ri), "mapping found");
(*testout) << rule->GetNP(1) << " = " << lfaces.Get(1).GetNP() << endl; (*testout) << rule->GetNP(1) << " = " << lfaces.Get(1).GetNP() << endl;
@ -527,7 +531,7 @@ int Meshing3 :: ApplyRules
// check mapedges: // 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), INDEX_2 in2(pmap.Get(rule->GetEdge(i).i1),
pmap.Get(rule->GetEdge(i).i2)); pmap.Get(rule->GetEdge(i).i2));
@ -537,12 +541,12 @@ int Meshing3 :: ApplyRules
// check prism edges: // check prism edges:
for (i = 1; i <= rule->GetNE(); i++) for (int i = 1; i <= rule->GetNE(); i++)
{ {
const Element & el = rule->GetElement (i); const Element & el = rule->GetElement (i);
if (el.GetType() == PRISM) 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)), INDEX_2 in2(pmap.Get(el.PNum(j)),
pmap.Get(el.PNum(j+3))); pmap.Get(el.PNum(j+3)));
@ -554,7 +558,7 @@ int Meshing3 :: ApplyRules
{ {
if (loktestmode) if (loktestmode)
(*testout) << "map pyramid, rule = " << rule->Name() << endl; (*testout) << "map pyramid, rule = " << rule->Name() << endl;
for (j = 1; j <= 2; j++) for (int j = 1; j <= 2; j++)
{ {
INDEX_2 in2; INDEX_2 in2;
if (j == 1) 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); fmapi.Set(i, 0);
@ -599,9 +603,9 @@ int Meshing3 :: ApplyRules
newu.SetSize (3 * (rule->GetNP() - rule->GetNOldP())); newu.SetSize (3 * (rule->GetNP() - rule->GetNOldP()));
allp.SetSize (3 * rule->GetNP()); 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); const Point3d & rp = rule->GetPoint(i);
oldu (3*i-3) = lp.X()-rp.X(); oldu (3*i-3) = lp.X()-rp.X();
oldu (3*i-2) = lp.Y()-rp.Y(); 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->GetNP()-rule->GetNOldP());
int idiff = 3 * 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); const Point3d & rp = rule->GetPoint(i);
allp (3*i-3) = rp.X() + newu(3*i-3 - idiff); allp (3*i-3) = rp.X() + newu(3*i-3 - idiff);
@ -644,14 +648,14 @@ int Meshing3 :: ApplyRules
{ {
const Array<Point3d> & fz = rule->GetTransFreeZone(); const Array<Point3d> & fz = rule->GetTransFreeZone();
(*testout) << "Freezone: " << endl; (*testout) << "Freezone: " << endl;
for (i = 1; i <= fz.Size(); i++) for (int i = 1; i <= fz.Size(); i++)
(*testout) << fz.Get(i) << endl; (*testout) << fz.Get(i) << endl;
} }
// check freezone: // check freezone:
for (i = 1; i <= lpoints.Size(); i++) for (int i = 1; i <= lpoints.Size(); i++)
{ {
if ( !pused.Get(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<int> lpi(4); static Array<int> lpi(4);
@ -692,7 +696,7 @@ int Meshing3 :: ApplyRules
for (li = 1; li <= lfacei.GetNP(); li++) for (li = 1; li <= lfacei.GetNP(); li++)
{ {
int lpii = 0; int lpii = 0;
int pi = lfacei.PNum(li); PointIndex pi = lfacei.PNum(li);
for (lj = 1; lj <= rule->GetNOldP(); lj++) for (lj = 1; lj <= rule->GetNOldP(); lj++)
if (pmap.Get(lj) == pi) if (pmap.Get(lj) == pi)
lpii = lj; lpii = lj;
@ -704,19 +708,19 @@ int Meshing3 :: ApplyRules
{ {
triin = rule->IsTriangleInFreeZone triin = rule->IsTriangleInFreeZone
( (
lpoints.Get(lfacei.PNum(1)), lpoints[lfacei.PNum(1)],
lpoints.Get(lfacei.PNum(2)), lpoints[lfacei.PNum(2)],
lpoints.Get(lfacei.PNum(3)), lpi, 1 lpoints[lfacei.PNum(3)], lpi, 1
); );
} }
else else
{ {
triin = rule->IsQuadInFreeZone triin = rule->IsQuadInFreeZone
( (
lpoints.Get(lfacei.PNum(1)), lpoints[lfacei.PNum(1)],
lpoints.Get(lfacei.PNum(2)), lpoints[lfacei.PNum(2)],
lpoints.Get(lfacei.PNum(3)), lpoints[lfacei.PNum(3)],
lpoints.Get(lfacei.PNum(4)), lpoints[lfacei.PNum(4)],
lpi, 1 lpi, 1
); );
} }
@ -741,7 +745,7 @@ int Meshing3 :: ApplyRules
<< lfaces.Get(i).PNum(3) << " - " << lfaces.Get(i).PNum(3) << " - "
<< lfaces.Get(i).PNum(4) << endl; << lfaces.Get(i).PNum(4) << endl;
for (int lj = 1; lj <= lfaces.Get(i).GetNP(); lj++) 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; (*testout) << endl;
@ -759,9 +763,9 @@ int Meshing3 :: ApplyRules
<< lfacei.PNum(2) << " - " << lfacei.PNum(2) << " - "
<< lfacei.PNum(3) << lfacei.PNum(3)
<< ", or " << ", or "
<< lpoints.Get(lfacei.PNum(1)) << " - " << lpoints[lfacei.PNum(1)] << " - "
<< lpoints.Get(lfacei.PNum(2)) << " - " << lpoints[lfacei.PNum(2)] << " - "
<< lpoints.Get(lfacei.PNum(3)) << lpoints[lfacei.PNum(3)]
<< endl; << endl;
(*testout) << "lpi = " << lpi.Get(1) << ", " (*testout) << "lpi = " << lpi.Get(1) << ", "
<< lpi.Get(2) << ", " << lpi.Get(3) << endl; << lpi.Get(2) << ", " << lpi.Get(3) << endl;
@ -773,10 +777,10 @@ int Meshing3 :: ApplyRules
<< lfacei.PNum(3) << " - " << lfacei.PNum(3) << " - "
<< lfacei.PNum(4) << lfacei.PNum(4)
<< ", or " << ", or "
<< lpoints.Get(lfacei.PNum(1)) << " - " << lpoints[lfacei.PNum(1)] << " - "
<< lpoints.Get(lfacei.PNum(2)) << " - " << lpoints[lfacei.PNum(2)] << " - "
<< lpoints.Get(lfacei.PNum(3)) << " - " << lpoints[lfacei.PNum(3)] << " - "
<< lpoints.Get(lfacei.PNum(4)) << lpoints[lfacei.PNum(4)]
<< endl; << endl;
sprintf (problems.Elem(ri), "triangle (%d, %d, %d) in Freezone", sprintf (problems.Elem(ri), "triangle (%d, %d, %d) in Freezone",
@ -786,13 +790,13 @@ int Meshing3 :: ApplyRules
} }
hc = 0; 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() && if (rule->GetPointNr(k, 1) <= rule->GetNOldP() &&
rule->GetPointNr(k, 2) <= rule->GetNOldP() && rule->GetPointNr(k, 2) <= rule->GetNOldP() &&
rule->GetPointNr(k, 3) <= 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)) && 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+1) == pmap.Get(rule->GetPointNr(k, 3)) &&
lfaces.Get(i).PNumMod(j+2) == pmap.Get(rule->GetPointNr(k, 2))) lfaces.Get(i).PNumMod(j+2) == pmap.Get(rule->GetPointNr(k, 2)))
@ -839,9 +843,9 @@ int Meshing3 :: ApplyRules
if (ok) if (ok)
{ {
err = 0; 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; if (hf > err) err = hf;
} }
@ -856,40 +860,37 @@ int Meshing3 :: ApplyRules
// newu = rule->GetOldUToNewU() * oldu; // newu = rule->GetOldUToNewU() * oldu;
// set new points: // set new points:
int oldnp = rule->GetNOldP();
int noldlp = lpoints.Size();
int noldlf = lfaces.Size();
oldnp = rule->GetNOldP(); for (int i = oldnp + 1; i <= rule->GetNP(); i++)
noldlp = lpoints.Size();
noldlf = lfaces.Size();
for (i = oldnp + 1; i <= rule->GetNP(); i++)
{ {
np = rule->GetPoint(i); np = rule->GetPoint(i);
np.X() += newu (3 * (i-oldnp) - 3); np.X() += newu (3 * (i-oldnp) - 3);
np.Y() += newu (3 * (i-oldnp) - 2); np.Y() += newu (3 * (i-oldnp) - 2);
np.Z() += newu (3 * (i-oldnp) - 1); np.Z() += newu (3 * (i-oldnp) - 1);
lpoints.Append (np);
pmap.Elem(i) = lpoints.Append (np); pmap.Elem(i) = lpoints.Size()-1+PointIndex::BASE;
} }
// Set new Faces: // 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)) if (!fmapi.Get(i))
{ {
MiniElement2d nface(rule->GetNP(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)); nface.PNum(j) = pmap.Get(rule->GetPointNr (i, j));
lfaces.Append (nface); lfaces.Append (nface);
} }
// Delete old Faces: // 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))); 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)) if (fmapi.Get(i))
{ {
delfaces.Append (fmapi.Get(i)); delfaces.Append (fmapi.Get(i));
@ -898,17 +899,17 @@ int Meshing3 :: ApplyRules
// check orientation // check orientation
for (i = 1; i <= rule->GetNO() && ok; i++) for (int i = 1; i <= rule->GetNO() && ok; i++)
{ {
const fourint * fouri; const fourint * fouri;
fouri = &rule->GetOrientation(i); fouri = &rule->GetOrientation(i);
Vec3d v1 (lpoints.Get(pmap.Get(fouri->i1)), Vec3d v1 (lpoints[pmap.Get(fouri->i1)],
lpoints.Get(pmap.Get(fouri->i2))); lpoints[pmap.Get(fouri->i2)]);
Vec3d v2 (lpoints.Get(pmap.Get(fouri->i1)), Vec3d v2 (lpoints[pmap.Get(fouri->i1)],
lpoints.Get(pmap.Get(fouri->i3))); lpoints[pmap.Get(fouri->i3)]);
Vec3d v3 (lpoints.Get(pmap.Get(fouri->i1)), Vec3d v3 (lpoints[pmap.Get(fouri->i1)],
lpoints.Get(pmap.Get(fouri->i4))); lpoints[pmap.Get(fouri->i4)]);
Vec3d n; Vec3d n;
Cross (v1, v2, n); Cross (v1, v2, n);
@ -927,7 +928,7 @@ int Meshing3 :: ApplyRules
// new points in free-zone ? // 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 (!rule->IsInFreeZone (lpoints.Get(pmap.Get(i))))
{ {
if (loktestmode) if (loktestmode)
@ -942,10 +943,10 @@ int Meshing3 :: ApplyRules
// insert new elements // insert new elements
for (i = 1; i <= rule->GetNE(); i++) for (int i = 1; i <= rule->GetNE(); i++)
{ {
elements.Append (rule->GetElement(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)); elements.Elem(i).PNum(j) = pmap.Get(elements.Get(i).PNum(j));
} }
@ -953,9 +954,9 @@ int Meshing3 :: ApplyRules
// Calculate Element badness // Calculate Element badness
teterr = 0; 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; if (hf > teterr) teterr = hf;
} }
@ -978,29 +979,29 @@ int Meshing3 :: ApplyRules
double oldlen = 0; double oldlen = 0;
double newlen = 0; double newlen = 0;
for (i = 1; i <= rule->GetNDelF(); i++) for (int i = 1; i <= rule->GetNDelF(); i++)
{ {
const Element2d & face = const Element2d & face =
rule->GetFace (rule->GetDelFace(i)); rule->GetFace (rule->GetDelFace(i));
for (j = 1; j <= 3; j++) for (int j = 1; j <= 3; j++)
{ {
const Point3d & p1 = const Point3d & p1 =
lpoints.Get(pmap.Get(face.PNumMod(j))); lpoints[pmap.Get(face.PNumMod(j))];
const Point3d & p2 = const Point3d & p2 =
lpoints.Get(pmap.Get(face.PNumMod(j+1))); lpoints[pmap.Get(face.PNumMod(j+1))];
oldlen += Dist(p1, p2); 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); const Element2d & face = rule->GetFace (i);
for (j = 1; j <= 3; j++) for (int j = 1; j <= 3; j++)
{ {
const Point3d & p1 = const Point3d & p1 =
lpoints.Get(pmap.Get(face.PNumMod(j))); lpoints[pmap.Get(face.PNumMod(j))];
const Point3d & p2 = const Point3d & p2 =
lpoints.Get(pmap.Get(face.PNumMod(j+1))); lpoints[pmap.Get(face.PNumMod(j+1))];
newlen += Dist(p1, p2); newlen += Dist(p1, p2);
} }
} }
@ -1057,7 +1058,7 @@ int Meshing3 :: ApplyRules
if (testmode) if (testmode)
{ {
for (i = 1; i <= rule->GetNOldP(); i++) for (int i = 1; i <= rule->GetNOldP(); i++)
{ {
(*testout) << "P" << i << ": Ref: " (*testout) << "P" << i << ": Ref: "
<< rule->GetPoint (i) << " is: " << rule->GetPoint (i) << " is: "
@ -1066,19 +1067,19 @@ int Meshing3 :: ApplyRules
} }
tempnewpoints.SetSize (0); 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)); tempnewpoints.Append (lpoints.Get(i));
tempnewfaces.SetSize (0); 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)); tempnewfaces.Append (lfaces.Get(i));
tempdelfaces.SetSize (0); tempdelfaces.SetSize (0);
for (i = 1; i <= delfaces.Size(); i++) for (int i = 1; i <= delfaces.Size(); i++)
tempdelfaces.Append (delfaces.Get(i)); tempdelfaces.Append (delfaces.Get(i));
tempelements.SetSize (0); tempelements.SetSize (0);
for (i = 1; i <= elements.Size(); i++) for (int i = 1; i <= elements.Size(); i++)
tempelements.Append (elements.Get(i)); tempelements.Append (elements.Get(i));
} }
@ -1096,15 +1097,13 @@ int Meshing3 :: ApplyRules
nfok = rule->GetNOldF(); 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); int refpi = rule->GetPointNr (nfok, j);
pused.Elem(pmap.Get(refpi))--; pused[pmap.Get(refpi)]--;
if (pused.Get(pmap.Get(refpi)) == 0) if (pused[pmap.Get(refpi)] == 0)
{ pmap.Elem(refpi).Invalidate();
pmap.Set(refpi, 0);
}
} }
} }
@ -1115,15 +1114,32 @@ int Meshing3 :: ApplyRules
if (found) if (found)
{ {
/*
for (i = 1; i <= tempnewpoints.Size(); i++) for (i = 1; i <= tempnewpoints.Size(); i++)
lpoints.Append (tempnewpoints.Get(i)); lpoints.Append (tempnewpoints.Get(i));
*/
for (Point3d p : tempnewpoints)
lpoints.Append(p);
/*
for (i = 1; i <= tempnewfaces.Size(); i++) for (i = 1; i <= tempnewfaces.Size(); i++)
if (tempnewfaces.Get(i).PNum(1)) if (tempnewfaces.Get(i).PNum(1))
lfaces.Append (tempnewfaces.Get(i)); 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++) for (i = 1; i <= tempdelfaces.Size(); i++)
delfaces.Append (tempdelfaces.Get(i)); delfaces.Append (tempdelfaces.Get(i));
*/
for (int i : tempdelfaces.Range())
delfaces.Append (tempdelfaces[i]);
/*
for (i = 1; i <= tempelements.Size(); i++) for (i = 1; i <= tempelements.Size(); i++)
elements.Append (tempelements.Get(i)); elements.Append (tempelements.Get(i));
*/
for (int i : tempelements.Range())
elements.Append (tempelements[i]);
} }
retminerr = minerr; retminerr = minerr;

26
libsrc/meshing/smoothing3.cpp
View File

@ -1422,19 +1422,18 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
if(lochfunc) if(lochfunc)
{ {
for(int i=1; i<=points.Size(); i++) for (PointIndex pi : points.Range())
pointh[i] = GetH(points.Get(i)); pointh[pi] = GetH(points[pi]);
} }
else else
{ {
pointh = 0; pointh = 0;
for(int i=0; i<GetNE(); i++) for (Element & el : VolumeElements())
{ {
const Element & el = VolumeElement(i+1);
double h = pow(el.Volume(points),1./3.); double h = pow(el.Volume(points),1./3.);
for(int j=1; j<=el.GetNV(); j++) for (PointIndex pi : el.PNums())
if(h > pointh[el.PNum(j)]) if (h > pointh[pi])
pointh[el.PNum(j)] = h; pointh[pi] = h;
} }
} }
@ -1455,21 +1454,15 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
const char * savetask = multithread.task; const char * savetask = multithread.task;
multithread.task = "Smooth Mesh"; 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 ( (*this)[pi].Type() == INNERPOINT && perrs[pi] > 0.01 * badmax)
{ {
if (multithread.terminate) if (multithread.terminate)
throw NgException ("Meshing stopped"); throw NgException ("Meshing stopped");
multithread.percent = 100.0 * (pi+1-PointIndex::BASE) / points.Size(); 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); if ( (pi+1-PointIndex::BASE) % printmod == 0) PrintDot (printdot);
double lh = pointh[pi]; double lh = pointh[pi];
@ -1503,7 +1496,6 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
} }
PrintDot ('\n'); PrintDot ('\n');
delete pf; delete pf;
multithread.task = savetask; multithread.task = savetask;

184
libsrc/meshing/topology.cpp
View File

@ -144,7 +144,7 @@ namespace netgen
{ swap (face.I2(), face.I3()); facedir += 2; } { swap (face.I2(), face.I3()); facedir += 2; }
if (face.I1() > face.I2()) if (face.I1() > face.I2())
{ swap (face.I1(), face.I2()); facedir += 4; } { swap (face.I1(), face.I2()); facedir += 4; }
if (face.I1() != v) continue; if (face.I1() != v) continue;
func (face, elnr, j, true, facedir); func (face, elnr, j, true, facedir);
@ -608,15 +608,15 @@ namespace netgen
{ {
case 3: case 3:
edges[elnr][loc_edge].nr = edgenum; edges[elnr][loc_edge].nr = edgenum;
edges[elnr][loc_edge].orient = edgedir; // edges[elnr][loc_edge].orient = edgedir;
break; break;
case 2: case 2:
surfedges[elnr][loc_edge].nr = edgenum; surfedges[elnr][loc_edge].nr = edgenum;
surfedges[elnr][loc_edge].orient = edgedir; // surfedges[elnr][loc_edge].orient = edgedir;
break; break;
case 1: case 1:
segedges[elnr].nr = edgenum; segedges[elnr].nr = edgenum;
segedges[elnr].orient = edgedir; // segedges[elnr].orient = edgedir;
break; break;
} }
}); });
@ -783,12 +783,12 @@ namespace netgen
if (volume) if (volume)
{ {
faces[elnr][j].fnr = facenum; faces[elnr][j].fnr = facenum;
faces[elnr][j].forient = facedir; // faces[elnr][j].forient = facedir;
} }
else else
{ {
surffaces[elnr].fnr = facenum; surffaces[elnr].fnr = facenum;
surffaces[elnr].forient = facedir; // surffaces[elnr].forient = facedir;
} }
}); });
} }
@ -1346,7 +1346,8 @@ namespace netgen
eorient.SetSize (ned); eorient.SetSize (ned);
for (int i = 1; i <= ned; i++) 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] > 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<int> & forient) const void MeshTopology :: GetElementFaceOrientations (int elnr, Array<int> & forient) const
@ -1354,8 +1355,9 @@ namespace netgen
int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType()); int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType());
forient.SetSize (nfa); forient.SetSize (nfa);
for (int i = 1; i <= nfa; i++) for (int i = 1; i <= nfa; i++)
forient.Elem(i) = faces.Get(elnr)[i-1].forient; // 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]-1) % 8;
forient.Elem(i) = GetElementFaceOrientation(elnr, i-1);
} }
@ -1377,7 +1379,8 @@ namespace netgen
*/ */
if (edges.Get(elnr)[i].nr == -1) return i; if (edges.Get(elnr)[i].nr == -1) return i;
eledges[i] = edges.Get(elnr)[i].nr+1; 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 else
@ -1432,7 +1435,8 @@ namespace netgen
*/ */
if (faces.Get(elnr)[i].fnr == -1) return i; if (faces.Get(elnr)[i].fnr == -1) return i;
elfaces[i] = faces.Get(elnr)[i].fnr+1; 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 else
@ -1484,13 +1488,15 @@ namespace netgen
eorient.SetSize (ned); eorient.SetSize (ned);
for (int i = 0; i < ned; i++) for (int i = 0; i < ned; i++)
// eorient[i] = (surfedges.Get(elnr)[i] > 0) ? 1 : -1; // 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 int MeshTopology :: GetSurfaceElementFaceOrientation (int elnr) const
{ {
// return (surffaces.Get(elnr)-1) % 8; // 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 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; if (surfedges.Get(elnr)[i].nr == -1) return i;
eledges[i] = surfedges.Get(elnr)[i].nr+1; 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; eledges[0] = segedges.Get(elnr).nr+1;
if (orient) 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; 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<int> & vertices) const void MeshTopology :: GetFaceVertices (int fnr, Array<int> & vertices) const
{ {
vertices.SetSize(4); vertices.SetSize(4);

18
libsrc/meshing/topology.hpp
View File

@ -17,14 +17,14 @@ namespace netgen
struct T_EDGE struct T_EDGE
{ {
int orient:1; // int orient:1;
int nr:31; // 0-based int nr; // 0-based
}; };
struct T_FACE struct T_FACE
{ {
int forient:3; // int forient:3;
int fnr:29; // 0-based int fnr; // 0-based
}; };
@ -87,7 +87,8 @@ public:
void GetSegmentEdge (int segnr, int & enr, int & orient) const void GetSegmentEdge (int segnr, int & enr, int & orient) const
{ {
enr = segedges.Get(segnr).nr+1; enr = segedges.Get(segnr).nr+1;
orient = segedges.Get(segnr).orient; // orient = segedges.Get(segnr).orient;
orient = GetSegmentEdgeOrientation(segnr);
} }
void GetElementEdges (int elnr, Array<int> & edges) const; void GetElementEdges (int elnr, Array<int> & edges) const;
@ -98,6 +99,13 @@ public:
int GetElementEdges (int elnr, int * edges, int * orient) const; int GetElementEdges (int elnr, int * edges, int * orient) const;
int GetElementFaces (int elnr, int * faces, 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<int> & vertices) const; void GetFaceVertices (int fnr, Array<int> & vertices) const;
void GetFaceVertices (int fnr, int * vertices) const; void GetFaceVertices (int fnr, int * vertices) const;
void GetEdgeVertices (int enr, int & v1, int & v2) const; void GetEdgeVertices (int enr, int & v1, int & v2) const;

3
libsrc/occ/occgeom.hpp
View File

@ -441,8 +441,7 @@ namespace netgen
// External access to the mesh generation functions within the OCC // External access to the mesh generation functions within the OCC
// subsystem (Not sure if this is the best way to implement this....!!) // subsystem (Not sure if this is the best way to implement this....!!)
extern int OCCGenerateMesh (OCCGeometry & occgeometry, shared_ptr<Mesh> & mesh, extern int OCCGenerateMesh (OCCGeometry & occgeometry, shared_ptr<Mesh> & mesh,
MeshingParameters & mparam, MeshingParameters & mparam);
int perfstepsstart, int perfstepsend);
extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh); extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh);

3
libsrc/stlgeom/stlgeom.cpp
View File

@ -455,7 +455,8 @@ int STLGeometry :: AddEdge(int ap1, int ap2)
STLEdge e(ap1,ap2); STLEdge e(ap1,ap2);
e.SetLeftTrig(GetLeftTrig(ap1,ap2)); e.SetLeftTrig(GetLeftTrig(ap1,ap2));
e.SetRightTrig(GetRightTrig(ap1,ap2)); e.SetRightTrig(GetRightTrig(ap1,ap2));
return edges.Append(e); edges.Append(e);
return edges.Size();
} }
void STLGeometry :: STLDoctorConfirmEdge() void STLGeometry :: STLDoctorConfirmEdge()

6
libsrc/stlgeom/stlgeom.hpp
View File

@ -302,11 +302,11 @@ namespace netgen
DLL_HEADER int GetNodeOfSelTrig() const; 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);} const Vec3d & GetNormal(int nr) const {return normals.Get(nr);}
void SetNormal(int nr, const Vec3d& n) {normals.Elem(nr) = n;} 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); int AddEdge(int p1, int p2);
STLEdge GetEdge(int nr) {return edges.Get(nr);} STLEdge GetEdge(int nr) {return edges.Get(nr);}
@ -434,7 +434,7 @@ namespace netgen
int ProjectOnWholeSurface (Point<3> & p3d) const; int ProjectOnWholeSurface (Point<3> & p3d) const;
int GetNLines() const {return lines.Size();} 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);} STLLine* GetLine(int nr) const {return lines.Get(nr);}
int GetLineP(int lnr, int pnr) const {return lines.Get(lnr)->PNum(pnr);} int GetLineP(int lnr, int pnr) const {return lines.Get(lnr)->PNum(pnr);}
int GetLineNP(int nr) const {return lines.Get(nr)->NP();} int GetLineNP(int nr) const {return lines.Get(nr)->NP();}

3
libsrc/stlgeom/stltool.cpp
View File

@ -19,7 +19,8 @@ int AddPointIfNotExists(Array<Point3d>& ap, const Point3d& p, double eps)
for (int i = 1; i <= ap.Size(); i++) for (int i = 1; i <= ap.Size(); i++)
if (Dist2(ap.Get(i),p) <= eps2 ) if (Dist2(ap.Get(i),p) <= eps2 )
return i; return i;
return ap.Append(p); ap.Append(p);
return ap.Size();
} }
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

3
libsrc/stlgeom/stltopology.cpp
View File

@ -592,7 +592,8 @@ void STLTopology :: FindNeighbourTrigs()
} }
else else
{ {
enr = topedges.Append (STLTopEdge (pi1, pi2, i, 0)); topedges.Append (STLTopEdge (pi1, pi2, i, 0));
enr = topedges.Size();
ht_topedges->Set (i2, enr); ht_topedges->Set (i2, enr);
trig.EdgeNum(j) = enr; trig.EdgeNum(j) = enr;
} }

2
libsrc/stlgeom/stltopology.hpp
View File

@ -293,7 +293,7 @@ public:
int GetNP() const { return points.Size(); } 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); } const Point<3> & GetPoint(int nr) const { return points.Get(nr); }
int GetPointNum (const Point<3> & p); int GetPointNum (const Point<3> & p);
void SetPoint(int nr, const Point<3> & p) { points.Elem(nr) = p; } void SetPoint(int nr, const Point<3> & p) { points.Elem(nr) = p; }

4
libsrc/visualization/mvdraw.cpp
View File

@ -63,8 +63,8 @@ namespace netgen
int VisualScene :: selface; int VisualScene :: selface;
int VisualScene :: selelement; int VisualScene :: selelement;
int VisualScene :: selpoint; PointIndex VisualScene :: selpoint;
int VisualScene :: selpoint2; PointIndex VisualScene :: selpoint2;
int VisualScene :: locpi; int VisualScene :: locpi;
int VisualScene :: seledge; int VisualScene :: seledge;

8
libsrc/visualization/mvdraw.hpp
View File

@ -26,8 +26,8 @@ namespace netgen
static int DLL_HEADER selface; static int DLL_HEADER selface;
static int selelement; static int selelement;
static int DLL_HEADER selpoint; static PointIndex DLL_HEADER selpoint;
static int selpoint2; static PointIndex selpoint2;
static int locpi; static int locpi;
static int DLL_HEADER seledge; static int DLL_HEADER seledge;
@ -235,8 +235,8 @@ namespace netgen
const Point3d & center, const Point3d & center,
const double rad, const double rad,
const int displaylist, const int displaylist,
int & selelement, int & selface, int & seledge, int & selpoint, int & selelement, int & selface, int & seledge, PointIndex & selpoint,
int & selpoint2, int & locpi); PointIndex & selpoint2, int & locpi);
} }

2
libsrc/visualization/vispar.hpp
View File

@ -76,7 +76,7 @@ public:
int drawcurveprojedge; int drawcurveprojedge;
int centerpoint; PointIndex centerpoint;
int drawelement; int drawelement;
// stl: // stl:

27
libsrc/visualization/vsmesh.cpp
View File

@ -440,12 +440,12 @@ namespace netgen
char buf[30]; char buf[30];
if (vispar.drawpointnumbers) 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()); glRasterPos3d (p.X(), p.Y(), p.Z());
sprintf (buf, "%d", i); sprintf (buf, "%d", int(pi));
// glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf); // glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf);
MyOpenGLText (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 // copy to be thread-safe
Element el = mesh->VolumeElement (i); Element el = mesh->VolumeElement (ei);
if ( (el.GetNP() == 4) || (el.GetNP() == 10)) if ( (el.GetNP() == 4) || (el.GetNP() == 10))
{ {
glBegin (GL_LINES); glBegin (GL_LINES);
@ -747,17 +747,16 @@ namespace netgen
} }
} }
for (SurfaceElementIndex sei : mesh->SurfaceElements().Range())
for (int i = 1; i <= mesh->GetNSE(); i++)
{ {
Element2d el = mesh->SurfaceElement(i); Element2d el = mesh->SurfaceElement(sei); // copy to be thread-safe
if (!el.BadElement()) if (!el.BadElement())
continue; continue;
int drawel = 1; bool drawel = true;
for (int j = 1; j <= el.GetNP(); j++) for (int j = 1; j <= el.GetNP(); j++)
if (!el.PNum(j)) if (!el.PNum(j).IsValid())
drawel = 0; drawel = false;
if (!drawel) if (!drawel)
continue; continue;
@ -3343,8 +3342,8 @@ namespace netgen
const Point3d & center, const Point3d & center,
const double rad, const double rad,
const int displaylist, const int displaylist,
int & selelement, int & selface, int & seledge, int & selpoint, int & selelement, int & selface, int & seledge, PointIndex & selpoint,
int & selpoint2, int & locpi) PointIndex & selpoint2, int & locpi)
{ {
auto mesh = vsmesh.GetMesh(); auto mesh = vsmesh.GetMesh();

4
libsrc/visualization/vssolution.cpp
View File

@ -4216,8 +4216,8 @@ namespace netgen
Point<3> p2 = locgrid[teti[pi2]]; Point<3> p2 = locgrid[teti[pi2]];
cppt.lami = p2 + edgelam[ednr] * (p1-p2); cppt.lami = p2 + edgelam[ednr] * (p1-p2);
pts.Append (cppt);
pnr = pts.Append (cppt)-1; pnr = pts.Size()-1;
edges.Set (pair, pnr); edges.Set (pair, pnr);
} }