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

improved usage of PointIndex

Browse Source
This commit is contained in:
Joachim Schoeberl 2013-04-02 20:29:53 +00:00
parent 15a1e07092
commit c46329ebf7
23 changed files with 334 additions and 394 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
libsrc/meshing/adfront3.cpp
View File

@ -86,8 +86,7 @@ AdFront3 :: ~AdFront3 ()
void AdFront3 :: GetPoints (Array<Point<3> > & apoints) const void AdFront3 :: GetPoints (Array<Point<3> > & apoints) const
{ {
for (PointIndex pi = PointIndex::BASE; for (PointIndex pi = points.Begin(); pi < points.End(); pi++)
pi < points.Size()+PointIndex::BASE; pi++)
apoints.Append (points[pi].P()); apoints.Append (points[pi].P());
} }
@ -106,7 +105,8 @@ PointIndex AdFront3 :: AddPoint (const Point<3> & p, PointIndex globind)
else else
{ {
points.Append (FrontPoint3 (p, globind)); points.Append (FrontPoint3 (p, globind));
return points.Size()-1+PointIndex::BASE; return --points.End();
// return points.Size()-1+PointIndex::BASE;
} }
} }
@ -301,9 +301,8 @@ void AdFront3 :: RebuildInternalTables ()
int np = points.Size(); int np = points.Size();
for (int i = PointIndex::BASE; for (PointIndex pi = points.Begin(); pi < points.End(); pi++)
i < np+PointIndex::BASE; i++) points[pi].cluster = pi;
points[i].cluster = i;
NgProfiler::StopTimer (timer_a); NgProfiler::StopTimer (timer_a);
NgProfiler::StartTimer (timer_b); NgProfiler::StartTimer (timer_b);
@ -400,9 +399,8 @@ void AdFront3 :: RebuildInternalTables ()
{ {
for (int i = 1; i <= faces.Size(); i++) for (int i = 1; i <= faces.Size(); i++)
faces.Elem(i).cluster = 1; faces.Elem(i).cluster = 1;
for (int i = PointIndex::BASE; for (PointIndex pi = points.Begin(); pi < points.End(); pi++)
i < points.Size()+PointIndex::BASE; i++) points[pi].cluster = 1;
points[i].cluster = 1;
} }
if (hashon) if (hashon)
@ -508,7 +506,6 @@ int AdFront3 :: GetLocals (int fstind,
INDEX i, j; INDEX i, j;
PointIndex pstind; PointIndex pstind;
INDEX pi;
Point3d midp, p0; Point3d midp, p0;
// static Array<int, PointIndex::BASE> invpindex; // static Array<int, PointIndex::BASE> invpindex;
@ -591,7 +588,7 @@ int AdFront3 :: GetLocals (int fstind,
for (i = 1; i <= locfaces.Size(); i++) for (i = 1; i <= locfaces.Size(); i++)
for (j = 1; j <= locfaces.Get(i).GetNP(); j++) for (j = 1; j <= locfaces.Get(i).GetNP(); j++)
{ {
pi = locfaces.Get(i).PNum(j); PointIndex pi = locfaces.Get(i).PNum(j);
invpindex[pi] = -1; invpindex[pi] = -1;
} }
@ -599,7 +596,7 @@ int AdFront3 :: GetLocals (int fstind,
{ {
for (j = 1; j <= locfaces.Get(i).GetNP(); j++) for (j = 1; j <= locfaces.Get(i).GetNP(); j++)
{ {
pi = locfaces.Get(i).PNum(j); PointIndex pi = locfaces.Get(i).PNum(j);
if (invpindex[pi] == -1) if (invpindex[pi] == -1)
{ {
pindex.Append (pi); pindex.Append (pi);
@ -703,12 +700,12 @@ void AdFront3 :: GetGroup (int fi,
invpindex.SetSize (points.Size()); invpindex.SetSize (points.Size());
for (i = 1; i <= points.Size(); i++) for (PointIndex pi = points.Begin(); pi < points.End(); pi++)
if (points.Get(i).Valid()) if (points.Get(i).Valid())
{ {
grouppoints.Append (points.Get(i).P()); grouppoints.Append (points[pi].P());
pindex.Append (i); pindex.Append (pi);
invpindex.Elem(i) = pindex.Size(); invpindex[pi] = pindex.Size();
} }
for (i = 1; i <= faces.Size(); i++) for (i = 1; i <= faces.Size(); i++)

4
libsrc/meshing/adfront3.hpp
View File

@ -176,7 +176,7 @@ public:
class AdFront3 class AdFront3
{ {
/// ///
Array<FrontPoint3, PointIndex::BASE> points; Array<FrontPoint3, PointIndex::BASE, PointIndex> points;
/// ///
Array<FrontFace> faces; Array<FrontFace> faces;
/// ///
@ -208,7 +208,7 @@ int rebuildcounter;
int lasti; int lasti;
/// minimal selection-value of baseelements /// minimal selection-value of baseelements
int minval; int minval;
Array<int, PointIndex::BASE> invpindex; Array<int, PointIndex::BASE, PointIndex> invpindex;
Array<char> pingroup; Array<char> pingroup;
/// ///

343
libsrc/meshing/bisect.cpp
View File

@ -6,18 +6,17 @@
namespace netgen namespace netgen
{ {
//#include "../interface/writeuser.hpp"
class MarkedTet; class MarkedTet;
class MarkedPrism; class MarkedPrism;
class MarkedIdentification; class MarkedIdentification;
class MarkedTri; class MarkedTri;
class MarkedQuad; class MarkedQuad;
typedef MoveableArray<MarkedTet> T_MTETS; typedef Array<MarkedTet> T_MTETS;
typedef MoveableArray<MarkedPrism> T_MPRISMS; typedef Array<MarkedPrism> T_MPRISMS;
typedef MoveableArray<MarkedIdentification> T_MIDS; typedef Array<MarkedIdentification> T_MIDS;
typedef MoveableArray<MarkedTri> T_MTRIS; typedef Array<MarkedTri> T_MTRIS;
typedef MoveableArray<MarkedQuad> T_MQUADS; typedef Array<MarkedQuad> T_MQUADS;
@ -855,8 +854,7 @@ namespace netgen
INDEX_2_CLOSED_HASHTABLE<int> & edgenumber, INDEX_2_CLOSED_HASHTABLE<int> & edgenumber,
MarkedTet & mt) MarkedTet & mt)
{ {
int i, j, k; for (int i = 0; i < 4; i++)
for (i = 0; i < 4; i++)
mt.pnums[i] = el[i]; mt.pnums[i] = el[i];
mt.marked = 0; mt.marked = 0;
@ -867,8 +865,8 @@ namespace netgen
int val = 0; int val = 0;
// find marked edge of tet: // find marked edge of tet:
for (i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
for (j = i+1; j < 4; j++) for (int j = i+1; j < 4; j++)
{ {
INDEX_2 i2(mt.pnums[i], mt.pnums[j]); INDEX_2 i2(mt.pnums[i], mt.pnums[j]);
i2.Sort(); i2.Sort();
@ -883,11 +881,11 @@ namespace netgen
// find marked edges of faces: // find marked edges of faces:
for (k = 0; k < 4; k++) for (int k = 0; k < 4; k++)
{ {
val = 0; val = 0;
for (i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
for (j = i+1; j < 4; j++) for (int j = i+1; j < 4; j++)
if (i != k && j != k) if (i != k && j != k)
{ {
INDEX_2 i2(mt.pnums[i], mt.pnums[j]); INDEX_2 i2(mt.pnums[i], mt.pnums[j]);
@ -910,19 +908,17 @@ namespace netgen
INDEX_2_CLOSED_HASHTABLE<int> & edgenumber, INDEX_2_CLOSED_HASHTABLE<int> & edgenumber,
MarkedPrism & mp) MarkedPrism & mp)
{ {
int i, j;
if (el.GetType() == PRISM || if (el.GetType() == PRISM ||
el.GetType() == PRISM12) el.GetType() == PRISM12)
{ {
for (i = 0; i < 6; i++) for (int i = 0; i < 6; i++)
mp.pnums[i] = el[i]; mp.pnums[i] = el[i];
} }
else if (el.GetType() == PYRAMID) else if (el.GetType() == PYRAMID)
{ {
static int map[6] = static int map[6] =
{ 1, 2, 5, 4, 3, 5 }; { 1, 2, 5, 4, 3, 5 };
for (i = 0; i < 6; i++) for (int i = 0; i < 6; i++)
mp.pnums[i] = el.PNum(map[i]); mp.pnums[i] = el.PNum(map[i]);
} }
else if (el.GetType() == TET || else if (el.GetType() == TET ||
@ -930,7 +926,7 @@ namespace netgen
{ {
static int map[6] = static int map[6] =
{ 1, 4, 3, 2, 4, 3 }; { 1, 4, 3, 2, 4, 3 };
for (i = 0; i < 6; i++) for (int i = 0; i < 6; i++)
mp.pnums[i] = el.PNum(map[i]); mp.pnums[i] = el.PNum(map[i]);
} }
@ -947,8 +943,8 @@ namespace netgen
mp.order = 1; mp.order = 1;
int val = 0; int val = 0;
for (i = 0; i < 2; i++) for (int i = 0; i < 2; i++)
for (j = i+1; j < 3; j++) for (int j = i+1; j < 3; j++)
{ {
INDEX_2 i2(mp.pnums[i], mp.pnums[j]); INDEX_2 i2(mp.pnums[i], mp.pnums[j]);
i2.Sort(); i2.Sort();
@ -1016,8 +1012,7 @@ namespace netgen
INDEX_2_CLOSED_HASHTABLE<int> & edgenumber, INDEX_2_CLOSED_HASHTABLE<int> & edgenumber,
MarkedTri & mt) MarkedTri & mt)
{ {
int i, j; for (int i = 0; i < 3; i++)
for (i = 0; i < 3; i++)
{ {
mt.pnums[i] = el[i]; mt.pnums[i] = el[i];
mt.pgeominfo[i] = el.GeomInfoPi (i+1); mt.pgeominfo[i] = el.GeomInfoPi (i+1);
@ -1030,8 +1025,8 @@ namespace netgen
mt.order = 1; mt.order = 1;
int val = 0; int val = 0;
for (i = 0; i < 2; i++) for (int i = 0; i < 2; i++)
for (j = i+1; j < 3; j++) for (int j = i+1; j < 3; j++)
{ {
INDEX_2 i2(mt.pnums[i], mt.pnums[j]); INDEX_2 i2(mt.pnums[i], mt.pnums[j]);
i2.Sort(); i2.Sort();
@ -1073,8 +1068,7 @@ namespace netgen
INDEX_2_CLOSED_HASHTABLE<int> & edgenumber, INDEX_2_CLOSED_HASHTABLE<int> & edgenumber,
MarkedQuad & mq) MarkedQuad & mq)
{ {
int i; for (int i = 0; i < 4; i++)
for (i = 0; i < 4; i++)
mq.pnums[i] = el[i]; mq.pnums[i] = el[i];
Swap (mq.pnums[2], mq.pnums[3]); Swap (mq.pnums[2], mq.pnums[3]);
@ -1217,8 +1211,6 @@ namespace netgen
*testout << "bisect tet " << oldtet << endl; *testout << "bisect tet " << oldtet << endl;
#endif #endif
int i, j, k;
// points vis a vis from tet-edge // points vis a vis from tet-edge
int vis1, vis2; int vis1, vis2;
@ -1233,10 +1225,10 @@ namespace netgen
// is tet of type P ? // is tet of type P ?
int istypep = 0; int istypep = 0;
for (i = 0; i < 4; i++) for (int i = 0; i < 4; i++)
{ {
int cnt = 0; int cnt = 0;
for (j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
if (oldtet.faceedges[j] == i) if (oldtet.faceedges[j] == i)
cnt++; cnt++;
if (cnt == 3) if (cnt == 3)
@ -1245,7 +1237,7 @@ namespace netgen
for (i = 0; i < 4; i++) for (int i = 0; i < 4; i++)
{ {
newtet1.pnums[i] = oldtet.pnums[i]; newtet1.pnums[i] = oldtet.pnums[i];
newtet2.pnums[i] = oldtet.pnums[i]; newtet2.pnums[i] = oldtet.pnums[i];
@ -1263,7 +1255,7 @@ namespace netgen
*testout << "newtet2,before = " << newtet2 << endl; *testout << "newtet2,before = " << newtet2 << endl;
#endif #endif
for (i = 0; i < 4; i++) for (int i = 0; i < 4; i++)
{ {
if (i == oldtet.tetedge1) if (i == oldtet.tetedge1)
{ {
@ -1272,10 +1264,10 @@ namespace netgen
newtet2.faceedges[vis1] = i; // cut faces newtet2.faceedges[vis1] = i; // cut faces
newtet2.faceedges[vis2] = i; newtet2.faceedges[vis2] = i;
j = 0; int j = 0;
while (j == i || j == oldtet.faceedges[i]) while (j == i || j == oldtet.faceedges[i])
j++; j++;
k = 6 - i - oldtet.faceedges[i] - j; int k = 6 - i - oldtet.faceedges[i] - j;
newtet2.tetedge1 = j; // tet-edge newtet2.tetedge1 = j; // tet-edge
newtet2.tetedge2 = k; newtet2.tetedge2 = k;
@ -1310,10 +1302,10 @@ namespace netgen
newtet1.faceedges[i] = oldtet.faceedges[i]; // inherited face newtet1.faceedges[i] = oldtet.faceedges[i]; // inherited face
newtet1.faceedges[vis1] = i; newtet1.faceedges[vis1] = i;
newtet1.faceedges[vis2] = i; newtet1.faceedges[vis2] = i;
j = 0; int j = 0;
while (j == i || j == oldtet.faceedges[i]) while (j == i || j == oldtet.faceedges[i])
j++; j++;
k = 6 - i - oldtet.faceedges[i] - j; int k = 6 - i - oldtet.faceedges[i] - j;
newtet1.tetedge1 = j; newtet1.tetedge1 = j;
newtet1.tetedge2 = k; newtet1.tetedge2 = k;
@ -1353,19 +1345,16 @@ namespace netgen
void BTBisectPrism (const MarkedPrism & oldprism, int newp1, int newp2, void BTBisectPrism (const MarkedPrism & oldprism, int newp1, int newp2,
MarkedPrism & newprism1, MarkedPrism & newprism2) MarkedPrism & newprism1, MarkedPrism & newprism2)
{ {
int i; for (int i = 0; i < 6; i++)
for (i = 0; i < 6; i++)
{ {
newprism1.pnums[i] = oldprism.pnums[i]; newprism1.pnums[i] = oldprism.pnums[i];
newprism2.pnums[i] = oldprism.pnums[i]; newprism2.pnums[i] = oldprism.pnums[i];
} }
int pe1, pe2; int pe1 = 0;
pe1 = 0;
if (pe1 == oldprism.markededge) if (pe1 == oldprism.markededge)
pe1++; pe1++;
pe2 = 3 - oldprism.markededge - pe1; int pe2 = 3 - oldprism.markededge - pe1;
newprism1.pnums[pe2] = newp1; newprism1.pnums[pe2] = newp1;
newprism1.pnums[pe2+3] = newp2; newprism1.pnums[pe2+3] = newp2;
@ -1390,7 +1379,7 @@ namespace netgen
void BTBisectIdentification (const MarkedIdentification & oldid, void BTBisectIdentification (const MarkedIdentification & oldid,
Array<int> & newp, Array<PointIndex> & newp,
MarkedIdentification & newid1, MarkedIdentification & newid1,
MarkedIdentification & newid2) MarkedIdentification & newid2)
{ {
@ -1440,9 +1429,7 @@ namespace netgen
void BTBisectTri (const MarkedTri & oldtri, int newp, const PointGeomInfo & newpgi, void BTBisectTri (const MarkedTri & oldtri, int newp, const PointGeomInfo & newpgi,
MarkedTri & newtri1, MarkedTri & newtri2) MarkedTri & newtri1, MarkedTri & newtri2)
{ {
int i; for (int i = 0; i < 3; i++)
for (i = 0; i < 3; i++)
{ {
newtri1.pnums[i] = oldtri.pnums[i]; newtri1.pnums[i] = oldtri.pnums[i];
newtri1.pgeominfo[i] = oldtri.pgeominfo[i]; newtri1.pgeominfo[i] = oldtri.pgeominfo[i];
@ -1450,11 +1437,10 @@ namespace netgen
newtri2.pgeominfo[i] = oldtri.pgeominfo[i]; newtri2.pgeominfo[i] = oldtri.pgeominfo[i];
} }
int pe1, pe2; int pe1 = 0;
pe1 = 0;
if (pe1 == oldtri.markededge) if (pe1 == oldtri.markededge)
pe1++; pe1++;
pe2 = 3 - oldtri.markededge - pe1; int pe2 = 3 - oldtri.markededge - pe1;
newtri1.pnums[pe2] = newp; newtri1.pnums[pe2] = newp;
newtri1.pgeominfo[pe2] = newpgi; newtri1.pgeominfo[pe2] = newpgi;
@ -1477,8 +1463,6 @@ namespace netgen
newtri1.order = oldtri.order; newtri1.order = oldtri.order;
newtri2.incorder = 0; newtri2.incorder = 0;
newtri2.order = oldtri.order; newtri2.order = oldtri.order;
} }
@ -1487,9 +1471,7 @@ namespace netgen
int newp2, const PointGeomInfo & npgi2, int newp2, const PointGeomInfo & npgi2,
MarkedQuad & newquad1, MarkedQuad & newquad2) MarkedQuad & newquad1, MarkedQuad & newquad2)
{ {
int i; for (int i = 0; i < 4; i++)
for (i = 0; i < 4; i++)
{ {
newquad1.pnums[i] = oldquad.pnums[i]; newquad1.pnums[i] = oldquad.pnums[i];
newquad1.pgeominfo[i] = oldquad.pgeominfo[i]; newquad1.pgeominfo[i] = oldquad.pgeominfo[i];
@ -1581,12 +1563,10 @@ namespace netgen
int MarkHangingIdentifications(T_MIDS & mids, int MarkHangingIdentifications(T_MIDS & mids,
const INDEX_2_CLOSED_HASHTABLE<int> & cutedges) const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges)
{ {
int i, j;
int hanging = 0; int hanging = 0;
for (i = 1; i <= mids.Size(); i++) for (int i = 1; i <= mids.Size(); i++)
{ {
if (mids.Elem(i).marked) if (mids.Elem(i).marked)
{ {
@ -1595,8 +1575,7 @@ namespace netgen
} }
const int np = mids.Get(i).np; const int np = mids.Get(i).np;
for(int j = 0; j < np; j++)
for(j = 0; j < np; j++)
{ {
INDEX_2 edge1(mids.Get(i).pnums[j], INDEX_2 edge1(mids.Get(i).pnums[j],
mids.Get(i).pnums[(j+1) % np]); mids.Get(i).pnums[(j+1) % np]);
@ -1678,12 +1657,10 @@ namespace netgen
int MarkHangingTets (T_MTETS & mtets, int MarkHangingTets (T_MTETS & mtets,
const INDEX_2_CLOSED_HASHTABLE<int> & cutedges) const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges)
{ {
int i, j, k;
int hanging = 0; int hanging = 0;
for (i = 1; i <= mtets.Size(); i++) for (int i = 1; i <= mtets.Size(); i++)
{ {
MarkedTet & teti = mtets.Elem(i); MarkedTet & teti = mtets.Elem(i);
@ -1693,8 +1670,8 @@ namespace netgen
continue; continue;
} }
for (j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
for (k = j+1; k < 4; k++) for (int k = j+1; k < 4; k++)
{ {
INDEX_2 edge(teti.pnums[j], INDEX_2 edge(teti.pnums[j],
teti.pnums[k]); teti.pnums[k]);
@ -1713,12 +1690,10 @@ namespace netgen
int MarkHangingPrisms (T_MPRISMS & mprisms, int MarkHangingPrisms (T_MPRISMS & mprisms,
const INDEX_2_CLOSED_HASHTABLE<int> & cutedges) const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges)
{ {
int i, j, k;
int hanging = 0; int hanging = 0;
for (i = 1; i <= mprisms.Size(); i++) for (int i = 1; i <= mprisms.Size(); i++)
{ {
if (mprisms.Elem(i).marked) if (mprisms.Elem(i).marked)
{ {
@ -1726,8 +1701,8 @@ namespace netgen
continue; continue;
} }
for (j = 0; j < 2; j++) for (int j = 0; j < 2; j++)
for (k = j+1; k < 3; k++) for (int k = j+1; k < 3; k++)
{ {
INDEX_2 edge1(mprisms.Get(i).pnums[j], INDEX_2 edge1(mprisms.Get(i).pnums[j],
mprisms.Get(i).pnums[k]); mprisms.Get(i).pnums[k]);
@ -1749,20 +1724,18 @@ namespace netgen
int MarkHangingTris (T_MTRIS & mtris, int MarkHangingTris (T_MTRIS & mtris,
const INDEX_2_CLOSED_HASHTABLE<int> & cutedges) const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges)
{ {
int i, j, k;
int hanging = 0; int hanging = 0;
for (i = 1; i <= mtris.Size(); i++) for (int i = 1; i <= mtris.Size(); i++)
{ {
if (mtris.Get(i).marked) if (mtris.Get(i).marked)
{ {
hanging = 1; hanging = 1;
continue; continue;
} }
for (j = 0; j < 2; j++) for (int j = 0; j < 2; j++)
for (k = j+1; k < 3; k++) for (int k = j+1; k < 3; k++)
{ {
INDEX_2 edge(mtris.Get(i).pnums[j], INDEX_2 edge(mtris.Get(i).pnums[j],
mtris.Get(i).pnums[k]); mtris.Get(i).pnums[k]);
@ -1780,12 +1753,10 @@ namespace netgen
int MarkHangingQuads (T_MQUADS & mquads, int MarkHangingQuads (T_MQUADS & mquads,
const INDEX_2_CLOSED_HASHTABLE<int> & cutedges) const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges)
{ {
int i;
int hanging = 0; int hanging = 0;
for (i = 1; i <= mquads.Size(); i++) for (int i = 1; i <= mquads.Size(); i++)
{ {
if (mquads.Elem(i).marked) if (mquads.Elem(i).marked)
{ {
@ -1834,11 +1805,10 @@ namespace netgen
void ConnectToNodeRec (int node, int tonode, void ConnectToNodeRec (int node, int tonode,
const TABLE<int> & conto, Array<int> & connecttonode) const TABLE<int> & conto, Array<int> & connecttonode)
{ {
int i, n2;
// (*testout) << "connect " << node << " to " << tonode << endl; // (*testout) << "connect " << node << " to " << tonode << endl;
for (i = 1; i <= conto.EntrySize(node); i++) for (int i = 1; i <= conto.EntrySize(node); i++)
{ {
n2 = conto.Get(node, i); int n2 = conto.Get(node, i);
if (!connecttonode.Get(n2)) if (!connecttonode.Get(n2))
{ {
connecttonode.Elem(n2) = tonode; connecttonode.Elem(n2) = tonode;
@ -1944,16 +1914,15 @@ namespace netgen
const Array< Array<int,PointIndex::BASE>* > & idmaps, const Array< Array<int,PointIndex::BASE>* > & idmaps,
const string & refinfofile) const string & refinfofile)
{ {
mtets.SetName ("bisection, tets"); // mtets.SetName ("bisection, tets");
mprisms.SetName ("bisection, prisms"); // mprisms.SetName ("bisection, prisms");
mtris.SetName ("bisection, trigs"); // mtris.SetName ("bisection, trigs");
mquads.SetName ("bisection, quads"); // nmquads.SetName ("bisection, quads");
mids.SetName ("bisection, identifications"); // mids.SetName ("bisection, identifications");
//int np = mesh.GetNP(); //int np = mesh.GetNP();
int ne = mesh.GetNE(); int ne = mesh.GetNE();
int nse = mesh.GetNSE(); int nse = mesh.GetNSE();
int i, j, k, l, m;
/* /*
if (mtets.Size() + mprisms.Size() == mesh.GetNE()) if (mtets.Size() + mprisms.Size() == mesh.GetNE())
@ -1983,13 +1952,13 @@ namespace netgen
INDEX_2_HASHTABLE<int> shortedges(100); INDEX_2_HASHTABLE<int> shortedges(100);
for (i = 1; i <= ne; i++) for (int i = 1; i <= ne; i++)
{ {
const Element & el = mesh.VolumeElement(i); const Element & el = mesh.VolumeElement(i);
if (el.GetType() == PRISM || if (el.GetType() == PRISM ||
el.GetType() == PRISM12) el.GetType() == PRISM12)
{ {
for (j = 1; j <= 3; j++) for (int j = 1; j <= 3; j++)
{ {
INDEX_2 se(el.PNum(j), el.PNum(j+3)); INDEX_2 se(el.PNum(j), el.PNum(j+3));
se.Sort(); se.Sort();
@ -2006,7 +1975,7 @@ namespace netgen
BTSortEdges (mesh, idmaps, edgenumber); BTSortEdges (mesh, idmaps, edgenumber);
for (i = 1; i <= ne; i++) for (int i = 1; i <= ne; i++)
{ {
const Element & el = mesh.VolumeElement(i); const Element & el = mesh.VolumeElement(i);
@ -2018,8 +1987,8 @@ namespace netgen
// if tet has short edge, it is handled as degenerated prism // if tet has short edge, it is handled as degenerated prism
int foundse = 0; int foundse = 0;
for (j = 1; j <= 3; j++) for (int j = 1; j <= 3; j++)
for (k = j+1; k <= 4; k++) for (int k = j+1; k <= 4; k++)
{ {
INDEX_2 se(el.PNum(j), el.PNum(k)); INDEX_2 se(el.PNum(j), el.PNum(k));
se.Sort(); se.Sort();
@ -2040,8 +2009,8 @@ namespace netgen
pi[2] = p3; pi[2] = p3;
pi[3] = p4; pi[3] = p4;
int cnt = 0; int cnt = 0;
for (l = 1; l <= 4; l++) for (int l = 1; l <= 4; l++)
for (m = 0; m < 3; m++) for (int m = 0; m < 3; m++)
if (pi[m] > pi[m+1]) if (pi[m] > pi[m+1])
{ {
Swap (pi[m], pi[m+1]); Swap (pi[m], pi[m+1]);
@ -2110,7 +2079,7 @@ namespace netgen
} }
} }
for (i = 1; i <= nse; i++) for (int i = 1; i <= nse; i++)
{ {
const Element2d & el = mesh.SurfaceElement(i); const Element2d & el = mesh.SurfaceElement(i);
if (el.GetType() == TRIG || if (el.GetType() == TRIG ||
@ -2128,7 +2097,7 @@ namespace netgen
} }
MarkedIdentification mi; MarkedIdentification mi;
for(j=0; j<idmaps.Size(); j++) for(int j=0; j<idmaps.Size(); j++)
if(BTDefineMarkedId(el, edgenumber, *idmaps[j], mi)) if(BTDefineMarkedId(el, edgenumber, *idmaps[j], mi))
mids.Append(mi); mids.Append(mi);
} }
@ -2138,10 +2107,10 @@ namespace netgen
mesh.mlparentelement.SetSize(ne); mesh.mlparentelement.SetSize(ne);
for (i = 1; i <= ne; i++) for (int i = 1; i <= ne; i++)
mesh.mlparentelement.Elem(i) = 0; mesh.mlparentelement.Elem(i) = 0;
mesh.mlparentsurfaceelement.SetSize(nse); mesh.mlparentsurfaceelement.SetSize(nse);
for (i = 1; i <= nse; i++) for (int i = 1; i <= nse; i++)
mesh.mlparentsurfaceelement.Elem(i) = 0; mesh.mlparentsurfaceelement.Elem(i) = 0;
if (printmessage_importance>0) if (printmessage_importance>0)
@ -2417,12 +2386,18 @@ namespace netgen
// const Array < Array<Element2d>* > & surfelements_before, // const Array < Array<Element2d>* > & surfelements_before,
// const Array < Array<int>* > & markedsurfelts_num) // const Array < Array<int>* > & markedsurfelts_num)
{ {
/*
T_MTETS mtets_old; mtets_old.Copy(mtets); T_MTETS mtets_old; mtets_old.Copy(mtets);
T_MPRISMS mprisms_old; mprisms_old.Copy(mprisms); T_MPRISMS mprisms_old; mprisms_old.Copy(mprisms);
T_MIDS mids_old; mids_old.Copy(mids); T_MIDS mids_old; mids_old.Copy(mids);
T_MTRIS mtris_old; mtris_old.Copy(mtris); T_MTRIS mtris_old; mtris_old.Copy(mtris);
T_MQUADS mquads_old; mquads_old.Copy(mquads); T_MQUADS mquads_old; mquads_old.Copy(mquads);
*/
T_MTETS mtets_old (mtets);
T_MPRISMS mprisms_old (mprisms);
T_MIDS mids_old (mids);
T_MTRIS mtris_old (mtris);
T_MQUADS mquads_old (mquads);
@ -2704,7 +2679,7 @@ namespace netgen
// int ne = mesh.GetNE(); // int ne = mesh.GetNE();
// int nse = mesh.GetNSE(); // int nse = mesh.GetNSE();
int i, j, l; // int i, j, l;
// int initnp = np; // int initnp = np;
// int maxsteps = 3; // int maxsteps = 3;
@ -2745,22 +2720,22 @@ namespace netgen
v_order = 0; v_order = 0;
for (ElementIndex ei = 0; ei < ne; ei++) for (ElementIndex ei = 0; ei < ne; ei++)
for (j = 0; j < mesh[ei].GetNP(); j++) for (int j = 0; j < mesh[ei].GetNP(); j++)
if (mesh[ei].GetOrder() > v_order[mesh[ei][j]]) if (mesh[ei].GetOrder() > v_order[mesh[ei][j]])
v_order[mesh[ei][j]] = mesh[ei].GetOrder(); v_order[mesh[ei][j]] = mesh[ei].GetOrder();
for (SurfaceElementIndex sei = 0; sei < nse; sei++) for (SurfaceElementIndex sei = 0; sei < nse; sei++)
for (j = 0; j < mesh[sei].GetNP(); j++) for (int j = 0; j < mesh[sei].GetNP(); j++)
if (mesh[sei].GetOrder() > v_order[mesh[sei][j]]) if (mesh[sei].GetOrder() > v_order[mesh[sei][j]])
v_order[mesh[sei][j]] = mesh[sei].GetOrder(); v_order[mesh[sei][j]] = mesh[sei].GetOrder();
for (ElementIndex ei = 0; ei < ne; ei++) for (ElementIndex ei = 0; ei < ne; ei++)
for (j = 0; j < mesh[ei].GetNP(); j++) for (int j = 0; j < mesh[ei].GetNP(); j++)
if (mesh[ei].GetOrder() < v_order[mesh[ei][j]]-1) if (mesh[ei].GetOrder() < v_order[mesh[ei][j]]-1)
mesh[ei].SetOrder(v_order[mesh[ei][j]]-1); mesh[ei].SetOrder(v_order[mesh[ei][j]]-1);
for (SurfaceElementIndex sei = 0; sei < nse; sei++) for (SurfaceElementIndex sei = 0; sei < nse; sei++)
for (j = 0; j < mesh[sei].GetNP(); j++) for (int j = 0; j < mesh[sei].GetNP(); j++)
if (mesh[sei].GetOrder() < v_order[mesh[sei][j]]-1) if (mesh[sei].GetOrder() < v_order[mesh[sei][j]]-1)
mesh[sei].SetOrder(v_order[mesh[sei][j]]-1); mesh[sei].SetOrder(v_order[mesh[sei][j]]-1);
@ -2773,11 +2748,11 @@ namespace netgen
// INDEX_2_HASHTABLE<int> cutedges(10 + 5 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size())); // INDEX_2_HASHTABLE<int> cutedges(10 + 5 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size()));
INDEX_2_CLOSED_HASHTABLE<int> cutedges(10 + 9 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size())); INDEX_2_CLOSED_HASHTABLE<PointIndex> cutedges(10 + 9 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size()));
bool noprojection = false; bool noprojection = false;
for (l = 1; l <= 1; l++) for (int l = 1; l <= 1; l++)
{ {
int marked = 0; int marked = 0;
if (opt.refinementfilename) if (opt.refinementfilename)
@ -2790,15 +2765,15 @@ namespace netgen
if(st == "refinementinfo") if(st == "refinementinfo")
// new version // new version
{ {
for(i=1; i<=mtets.Size(); i++) for(int i=1; i<=mtets.Size(); i++)
mtets.Elem(i).marked = 0; mtets.Elem(i).marked = 0;
for(i=1; i<=mprisms.Size(); i++) for(int i=1; i<=mprisms.Size(); i++)
mprisms.Elem(i).marked = 0; mprisms.Elem(i).marked = 0;
for(i=1; i<=mtris.Size(); i++) for(int i=1; i<=mtris.Size(); i++)
mtris.Elem(i).marked = 0; mtris.Elem(i).marked = 0;
for(i=1; i<=mquads.Size(); i++) for(int i=1; i<=mquads.Size(); i++)
mquads.Elem(i).marked = 0; mquads.Elem(i).marked = 0;
for(i=1; i<=mprisms.Size(); i++) for(int i=1; i<=mprisms.Size(); i++)
mids.Elem(i).marked = 0; mids.Elem(i).marked = 0;
inf >> st; inf >> st;
@ -2844,7 +2819,7 @@ namespace netgen
else if(st == "orthobrick") else if(st == "orthobrick")
{ {
double bounds[6]; double bounds[6];
for(i=0; i<6; i++) for(int i=0; i<6; i++)
inf >> bounds[i]; inf >> bounds[i];
int cnt = 0; int cnt = 0;
@ -2855,14 +2830,14 @@ namespace netgen
// //
Point<3> center(0,0,0); Point<3> center(0,0,0);
for(i=0; i<el.GetNP(); i++) for(int i=0; i<el.GetNP(); i++)
{ {
const MeshPoint & point = mesh[el[i]]; const MeshPoint & point = mesh[el[i]];
center(0) += point(0); center(0) += point(0);
center(1) += point(1); center(1) += point(1);
center(2) += point(2); center(2) += point(2);
} }
for(i=0; i<3; i++) for(int i=0; i<3; i++)
center(i) *= 1./double(el.GetNP()); center(i) *= 1./double(el.GetNP());
if(bounds[0] <= center(0) && center(0) <= bounds[3] && if(bounds[0] <= center(0) && center(0) <= bounds[3] &&
bounds[1] <= center(1) && center(1) <= bounds[4] && bounds[1] <= center(1) && center(1) <= bounds[4] &&
@ -2919,7 +2894,7 @@ namespace netgen
inf.open(opt.refinementfilename); inf.open(opt.refinementfilename);
char ch; char ch;
for (i = 1; i <= mtets.Size(); i++) for (int i = 1; i <= mtets.Size(); i++)
{ {
inf >> ch; inf >> ch;
if(!inf) if(!inf)
@ -2940,7 +2915,7 @@ namespace netgen
{ {
int cnttet = 0; int cnttet = 0;
int cntprism = 0; int cntprism = 0;
for (i = 1; i <= mesh.GetNE(); i++) for (int i = 1; i <= mesh.GetNE(); i++)
{ {
if (mesh.VolumeElement(i).GetType() == TET || if (mesh.VolumeElement(i).GetType() == TET ||
mesh.VolumeElement(i).GetType() == TET10) mesh.VolumeElement(i).GetType() == TET10)
@ -2963,7 +2938,7 @@ namespace netgen
} }
} }
else else
for (i = 1; i <= mtets.Size(); i++) for (int i = 1; i <= mtets.Size(); i++)
{ {
mtets.Elem(i).marked = mtets.Elem(i).marked =
3 * mesh.VolumeElement(i).TestRefinementFlag(); 3 * mesh.VolumeElement(i).TestRefinementFlag();
@ -2989,7 +2964,7 @@ namespace netgen
int cnttrig = 0; int cnttrig = 0;
int cntquad = 0; int cntquad = 0;
for (i = 1; i <= mesh.GetNSE(); i++) for (int i = 1; i <= mesh.GetNSE(); i++)
{ {
if (mesh.SurfaceElement(i).GetType() == TRIG || if (mesh.SurfaceElement(i).GetType() == TRIG ||
mesh.SurfaceElement(i).GetType() == TRIG6) mesh.SurfaceElement(i).GetType() == TRIG6)
@ -3064,26 +3039,26 @@ namespace netgen
{ {
PrintMessage(3,"refine p"); PrintMessage(3,"refine p");
for (i = 1; i <= mtets.Size(); i++) for (int i = 1; i <= mtets.Size(); i++)
mtets.Elem(i).incorder = mtets.Elem(i).marked ? 1 : 0; mtets.Elem(i).incorder = mtets.Elem(i).marked ? 1 : 0;
for (i = 1; i <= mtets.Size(); i++) for (int i = 1; i <= mtets.Size(); i++)
if (mtets.Elem(i).incorder) if (mtets.Elem(i).incorder)
mtets.Elem(i).marked = 0; mtets.Elem(i).marked = 0;
for (i = 1; i <= mprisms.Size(); i++) for (int i = 1; i <= mprisms.Size(); i++)
mprisms.Elem(i).incorder = mprisms.Elem(i).marked ? 1 : 0; mprisms.Elem(i).incorder = mprisms.Elem(i).marked ? 1 : 0;
for (i = 1; i <= mprisms.Size(); i++) for (int i = 1; i <= mprisms.Size(); i++)
if (mprisms.Elem(i).incorder) if (mprisms.Elem(i).incorder)
mprisms.Elem(i).marked = 0; mprisms.Elem(i).marked = 0;
for (i = 1; i <= mtris.Size(); i++) for (int i = 1; i <= mtris.Size(); i++)
mtris.Elem(i).incorder = mtris.Elem(i).marked ? 1 : 0; mtris.Elem(i).incorder = mtris.Elem(i).marked ? 1 : 0;
for (i = 1; i <= mtris.Size(); i++) for (int i = 1; i <= mtris.Size(); i++)
{ {
if (mtris.Elem(i).incorder) if (mtris.Elem(i).incorder)
mtris.Elem(i).marked = 0; mtris.Elem(i).marked = 0;
@ -3098,7 +3073,7 @@ namespace netgen
if (mesh.GetDimension() == 3) if (mesh.GetDimension() == 3)
{ {
for (i = 1; i <= mesh.GetNSeg(); i++) for (int i = 1; i <= mesh.GetNSeg(); i++)
{ {
const Segment & seg = mesh.LineSegment(i); const Segment & seg = mesh.LineSegment(i);
singv.Set (seg[0]); singv.Set (seg[0]);
@ -3118,10 +3093,10 @@ namespace netgen
// vertices with 2 different bnds // vertices with 2 different bnds
Array<int> bndind(np); Array<int> bndind(np);
bndind = 0; bndind = 0;
for (i = 1; i <= mesh.GetNSeg(); i++) for (int i = 1; i <= mesh.GetNSeg(); i++)
{ {
const Segment & seg = mesh.LineSegment(i); const Segment & seg = mesh.LineSegment(i);
for (j = 0; j < 2; j++) for (int j = 0; j < 2; j++)
{ {
int pi = (j == 0) ? seg[0] : seg[1]; int pi = (j == 0) ? seg[0] : seg[1];
if (bndind.Elem(pi) == 0) if (bndind.Elem(pi) == 0)
@ -3134,47 +3109,47 @@ namespace netgen
for (i = 1; i <= mtets.Size(); i++) for (int i = 1; i <= mtets.Size(); i++)
mtets.Elem(i).incorder = 1; mtets.Elem(i).incorder = 1;
for (i = 1; i <= mtets.Size(); i++) for (int i = 1; i <= mtets.Size(); i++)
{ {
if (!mtets.Elem(i).marked) if (!mtets.Elem(i).marked)
mtets.Elem(i).incorder = 0; mtets.Elem(i).incorder = 0;
for (j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
if (singv.Test (mtets.Elem(i).pnums[j])) if (singv.Test (mtets.Elem(i).pnums[j]))
mtets.Elem(i).incorder = 0; mtets.Elem(i).incorder = 0;
} }
for (i = 1; i <= mtets.Size(); i++) for (int i = 1; i <= mtets.Size(); i++)
if (mtets.Elem(i).incorder) if (mtets.Elem(i).incorder)
mtets.Elem(i).marked = 0; mtets.Elem(i).marked = 0;
for (i = 1; i <= mprisms.Size(); i++) for (int i = 1; i <= mprisms.Size(); i++)
mprisms.Elem(i).incorder = 1; mprisms.Elem(i).incorder = 1;
for (i = 1; i <= mprisms.Size(); i++) for (int i = 1; i <= mprisms.Size(); i++)
{ {
if (!mprisms.Elem(i).marked) if (!mprisms.Elem(i).marked)
mprisms.Elem(i).incorder = 0; mprisms.Elem(i).incorder = 0;
for (j = 0; j < 6; j++) for (int j = 0; j < 6; j++)
if (singv.Test (mprisms.Elem(i).pnums[j])) if (singv.Test (mprisms.Elem(i).pnums[j]))
mprisms.Elem(i).incorder = 0; mprisms.Elem(i).incorder = 0;
} }
for (i = 1; i <= mprisms.Size(); i++) for (int i = 1; i <= mprisms.Size(); i++)
if (mprisms.Elem(i).incorder) if (mprisms.Elem(i).incorder)
mprisms.Elem(i).marked = 0; mprisms.Elem(i).marked = 0;
for (i = 1; i <= mtris.Size(); i++) for (int i = 1; i <= mtris.Size(); i++)
mtris.Elem(i).incorder = 1; mtris.Elem(i).incorder = 1;
for (i = 1; i <= mtris.Size(); i++) for (int i = 1; i <= mtris.Size(); i++)
{ {
if (!mtris.Elem(i).marked) if (!mtris.Elem(i).marked)
mtris.Elem(i).incorder = 0; mtris.Elem(i).incorder = 0;
for (j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
if (singv.Test (mtris.Elem(i).pnums[j])) if (singv.Test (mtris.Elem(i).pnums[j]))
mtris.Elem(i).incorder = 0; mtris.Elem(i).incorder = 0;
} }
for (i = 1; i <= mtris.Size(); i++) for (int i = 1; i <= mtris.Size(); i++)
{ {
if (mtris.Elem(i).incorder) if (mtris.Elem(i).incorder)
mtris.Elem(i).marked = 0; mtris.Elem(i).marked = 0;
@ -3195,12 +3170,12 @@ namespace netgen
// refine volume elements // refine volume elements
int nel = mtets.Size(); int nel = mtets.Size();
for (i = 1; i <= nel; i++) for (int i = 1; i <= nel; i++)
if (mtets.Elem(i).marked) if (mtets.Elem(i).marked)
{ {
MarkedTet oldtet; MarkedTet oldtet;
MarkedTet newtet1, newtet2; MarkedTet newtet1, newtet2;
int newp; PointIndex newp;
oldtet = mtets.Get(i); oldtet = mtets.Get(i);
@ -3236,12 +3211,12 @@ namespace netgen
} }
int npr = mprisms.Size(); int npr = mprisms.Size();
for (i = 1; i <= npr; i++) for (int i = 1; i <= npr; i++)
if (mprisms.Elem(i).marked) if (mprisms.Elem(i).marked)
{ {
MarkedPrism oldprism; MarkedPrism oldprism;
MarkedPrism newprism1, newprism2; MarkedPrism newprism1, newprism2;
int newp1, newp2; PointIndex newp1, newp2;
oldprism = mprisms.Get(i); oldprism = mprisms.Get(i);
int pi1 = 0; int pi1 = 0;
@ -3284,11 +3259,11 @@ namespace netgen
} }
int nid = mids.Size(); int nid = mids.Size();
for (i = 1; i <= nid; i++) for (int i = 1; i <= nid; i++)
if (mids.Elem(i).marked) if (mids.Elem(i).marked)
{ {
MarkedIdentification oldid,newid1,newid2; MarkedIdentification oldid,newid1,newid2;
Array<int> newp; Array<PointIndex> newp;
oldid = mids.Get(i); oldid = mids.Get(i);
@ -3305,7 +3280,7 @@ namespace netgen
edges.Append(INDEX_2(oldid.pnums[(oldid.markededge+2)%oldid.np + oldid.np], edges.Append(INDEX_2(oldid.pnums[(oldid.markededge+2)%oldid.np + oldid.np],
oldid.pnums[(oldid.markededge+3)%oldid.np + oldid.np])); oldid.pnums[(oldid.markededge+3)%oldid.np + oldid.np]));
} }
for (j = 0; j < edges.Size(); j++) for (int j = 0; j < edges.Size(); j++)
{ {
edges[j].Sort(); edges[j].Sort();
@ -3337,7 +3312,7 @@ namespace netgen
int nsel = mtris.Size(); int nsel = mtris.Size();
for (i = 1; i <= nsel; i++) for (int i = 1; i <= nsel; i++)
if (mtris.Elem(i).marked) if (mtris.Elem(i).marked)
{ {
MarkedTri oldtri; MarkedTri oldtri;
@ -3389,12 +3364,12 @@ namespace netgen
} }
int nquad = mquads.Size(); int nquad = mquads.Size();
for (i = 1; i <= nquad; i++) for (int i = 1; i <= nquad; i++)
if (mquads.Elem(i).marked) if (mquads.Elem(i).marked)
{ {
MarkedQuad oldquad; MarkedQuad oldquad;
MarkedQuad newquad1, newquad2; MarkedQuad newquad1, newquad2;
int newp1, newp2; PointIndex newp1, newp2;
oldquad = mquads.Get(i); oldquad = mquads.Get(i);
/* /*
@ -3488,7 +3463,7 @@ namespace netgen
hangingedge = 0; hangingedge = 0;
int nseg = mesh.GetNSeg (); int nseg = mesh.GetNSeg ();
for (i = 1; i <= nseg; i++) for (int i = 1; i <= nseg; i++)
{ {
Segment & seg = mesh.LineSegment (i); Segment & seg = mesh.LineSegment (i);
INDEX_2 edge(seg[0], seg[1]); INDEX_2 edge(seg[0], seg[1]);
@ -3561,34 +3536,34 @@ namespace netgen
v_order = 0; v_order = 0;
if (mesh.GetDimension() == 3) if (mesh.GetDimension() == 3)
{ {
for (i = 1; i <= mtets.Size(); i++) for (int i = 1; i <= mtets.Size(); i++)
if (mtets.Elem(i).incorder) if (mtets.Elem(i).incorder)
mtets.Elem(i).order++; mtets.Elem(i).order++;
for (i = 0; i < mtets.Size(); i++) for (int i = 0; i < mtets.Size(); i++)
for (j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
if (int(mtets[i].order) > v_order.Elem(mtets[i].pnums[j])) if (int(mtets[i].order) > v_order.Elem(mtets[i].pnums[j]))
v_order.Elem(mtets[i].pnums[j]) = mtets[i].order; v_order.Elem(mtets[i].pnums[j]) = mtets[i].order;
for (i = 0; i < mtets.Size(); i++) for (int i = 0; i < mtets.Size(); i++)
for (j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
if (int(mtets[i].order) < v_order.Elem(mtets[i].pnums[j])-1) if (int(mtets[i].order) < v_order.Elem(mtets[i].pnums[j])-1)
mtets[i].order = v_order.Elem(mtets[i].pnums[j])-1; mtets[i].order = v_order.Elem(mtets[i].pnums[j])-1;
} }
else else
{ {
for (i = 1; i <= mtris.Size(); i++) for (int i = 1; i <= mtris.Size(); i++)
if (mtris.Elem(i).incorder) if (mtris.Elem(i).incorder)
{ {
mtris.Elem(i).order++; mtris.Elem(i).order++;
} }
for (i = 0; i < mtris.Size(); i++) for (int i = 0; i < mtris.Size(); i++)
for (j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
if (int(mtris[i].order) > v_order.Elem(mtris[i].pnums[j])) if (int(mtris[i].order) > v_order.Elem(mtris[i].pnums[j]))
v_order.Elem(mtris[i].pnums[j]) = mtris[i].order; v_order.Elem(mtris[i].pnums[j]) = mtris[i].order;
for (i = 0; i < mtris.Size(); i++) for (int i = 0; i < mtris.Size(); i++)
{ {
for (j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
if (int(mtris[i].order) < v_order.Elem(mtris[i].pnums[j])-1) if (int(mtris[i].order) < v_order.Elem(mtris[i].pnums[j])-1)
mtris[i].order = v_order.Elem(mtris[i].pnums[j])-1; mtris[i].order = v_order.Elem(mtris[i].pnums[j])-1;
} }
@ -3604,20 +3579,20 @@ namespace netgen
mesh.ClearVolumeElements(); mesh.ClearVolumeElements();
mesh.VolumeElements().SetAllocSize (mtets.Size()+mprisms.Size()); mesh.VolumeElements().SetAllocSize (mtets.Size()+mprisms.Size());
for (i = 1; i <= mtets.Size(); i++) for (int i = 1; i <= mtets.Size(); i++)
{ {
Element el(TET); Element el(TET);
el.SetIndex (mtets.Get(i).matindex); el.SetIndex (mtets.Get(i).matindex);
for (j = 1; j <= 4; j++) for (int j = 1; j <= 4; j++)
el.PNum(j) = mtets.Get(i).pnums[j-1]; el.PNum(j) = mtets.Get(i).pnums[j-1];
el.SetOrder (mtets.Get(i).order); el.SetOrder (mtets.Get(i).order);
mesh.AddVolumeElement (el); mesh.AddVolumeElement (el);
} }
for (i = 1; i <= mprisms.Size(); i++) for (int i = 1; i <= mprisms.Size(); i++)
{ {
Element el(PRISM); Element el(PRISM);
el.SetIndex (mprisms.Get(i).matindex); el.SetIndex (mprisms.Get(i).matindex);
for (j = 1; j <= 6; j++) for (int j = 1; j <= 6; j++)
el.PNum(j) = mprisms.Get(i).pnums[j-1]; el.PNum(j) = mprisms.Get(i).pnums[j-1];
el.SetOrder (mprisms.Get(i).order); el.SetOrder (mprisms.Get(i).order);
@ -3638,7 +3613,7 @@ namespace netgen
{ {
case 1: case 1:
{ {
for (j = 1; j <= 5; j++) for (int j = 1; j <= 5; j++)
el.PNum(j) = mprisms.Get(i).pnums[map[j-1]-1]; el.PNum(j) = mprisms.Get(i).pnums[map[j-1]-1];
el.SetType (PYRAMID); el.SetType (PYRAMID);
@ -3652,7 +3627,7 @@ namespace netgen
if (!deg1) map = tetmap1; if (!deg1) map = tetmap1;
if (!deg2) map = tetmap2; if (!deg2) map = tetmap2;
if (!deg3) map = tetmap3; if (!deg3) map = tetmap3;
for (j = 1; j <= 4; j++) for (int j = 1; j <= 4; j++)
el.PNum(j) = mprisms.Get(i).pnums[map[j-1]-1]; el.PNum(j) = mprisms.Get(i).pnums[map[j-1]-1];
/* /*
if (!deg1) el.PNum(4) = el.PNum(4); if (!deg1) el.PNum(4) = el.PNum(4);
@ -3669,23 +3644,23 @@ namespace netgen
} }
mesh.ClearSurfaceElements(); mesh.ClearSurfaceElements();
for (i = 1; i <= mtris.Size(); i++) for (int i = 1; i <= mtris.Size(); i++)
{ {
Element2d el(TRIG); Element2d el(TRIG);
el.SetIndex (mtris.Get(i).surfid); el.SetIndex (mtris.Get(i).surfid);
el.SetOrder (mtris.Get(i).order); el.SetOrder (mtris.Get(i).order);
for (j = 1; j <= 3; j++) for (int j = 1; j <= 3; j++)
{ {
el.PNum(j) = mtris.Get(i).pnums[j-1]; el.PNum(j) = mtris.Get(i).pnums[j-1];
el.GeomInfoPi(j) = mtris.Get(i).pgeominfo[j-1]; el.GeomInfoPi(j) = mtris.Get(i).pgeominfo[j-1];
} }
mesh.AddSurfaceElement (el); mesh.AddSurfaceElement (el);
} }
for (i = 1; i <= mquads.Size(); i++) for (int i = 1; i <= mquads.Size(); i++)
{ {
Element2d el(QUAD); Element2d el(QUAD);
el.SetIndex (mquads.Get(i).surfid); el.SetIndex (mquads.Get(i).surfid);
for (j = 1; j <= 4; j++) for (int j = 1; j <= 4; j++)
el.PNum(j) = mquads.Get(i).pnums[j-1]; el.PNum(j) = mquads.Get(i).pnums[j-1];
Swap (el.PNum(3), el.PNum(4)); Swap (el.PNum(3), el.PNum(4));
mesh.AddSurfaceElement (el); mesh.AddSurfaceElement (el);
@ -3699,7 +3674,7 @@ namespace netgen
if (mesh.mglevels <= 2) if (mesh.mglevels <= 2)
{ {
PrintMessage(4,"RESETTING mlbetweennodes"); PrintMessage(4,"RESETTING mlbetweennodes");
for (i = 1; i <= np; i++) for (int i = 1; i <= np; i++)
{ {
mesh.mlbetweennodes.Elem(i).I1() = 0; mesh.mlbetweennodes.Elem(i).I1() = 0;
mesh.mlbetweennodes.Elem(i).I2() = 0; mesh.mlbetweennodes.Elem(i).I2() = 0;
@ -3720,11 +3695,11 @@ namespace netgen
BitArray isnewpoint(np); BitArray isnewpoint(np);
isnewpoint.Clear(); isnewpoint.Clear();
for (i = 1; i <= cutedges.Size(); i++) for (int i = 1; i <= cutedges.Size(); i++)
if (cutedges.UsedPos(i)) if (cutedges.UsedPos(i))
{ {
INDEX_2 edge; INDEX_2 edge;
int newpi; PointIndex newpi;
cutedges.GetData (i, edge, newpi); cutedges.GetData (i, edge, newpi);
isnewpoint.Set(newpi); isnewpoint.Set(newpi);
mesh.mlbetweennodes.Elem(newpi) = edge; mesh.mlbetweennodes.Elem(newpi) = edge;
@ -3801,7 +3776,7 @@ namespace netgen
// update identification tables // update identification tables
for (i = 1; i <= mesh.GetIdentifications().GetMaxNr(); i++) for (int i = 1; i <= mesh.GetIdentifications().GetMaxNr(); i++)
{ {
Array<int,PointIndex::BASE> identmap; Array<int,PointIndex::BASE> identmap;
@ -3826,18 +3801,18 @@ namespace netgen
} }
*/ */
for (j = 1; j <= cutedges.Size(); j++) for (int j = 1; j <= cutedges.Size(); j++)
if (cutedges.UsedPos(j)) if (cutedges.UsedPos(j))
{ {
INDEX_2 i2; INDEX_2 i2;
int newpi; PointIndex newpi;
cutedges.GetData (j, i2, newpi); cutedges.GetData (j, i2, newpi);
INDEX_2 oi2(identmap.Get(i2.I1()), INDEX_2 oi2(identmap.Get(i2.I1()),
identmap.Get(i2.I2())); identmap.Get(i2.I2()));
oi2.Sort(); oi2.Sort();
if (cutedges.Used (oi2)) if (cutedges.Used (oi2))
{ {
int onewpi = cutedges.Get(oi2); PointIndex onewpi = cutedges.Get(oi2);
mesh.GetIdentifications().Add (newpi, onewpi, i); mesh.GetIdentifications().Add (newpi, onewpi, i);
} }
} }
@ -3913,8 +3888,8 @@ namespace netgen
// update id-maps // update id-maps
j=0; int j=0;
for(i=1; i<=mesh.GetIdentifications().GetMaxNr(); i++) for(int i=1; i<=mesh.GetIdentifications().GetMaxNr(); i++)
{ {
if(mesh.GetIdentifications().GetType(i) == Identifications::PERIODIC) if(mesh.GetIdentifications().GetType(i) == Identifications::PERIODIC)
{ {
@ -3969,7 +3944,7 @@ namespace netgen
for(i=0; i<idmaps.Size(); i++) for(int i=0; i<idmaps.Size(); i++)
delete idmaps[i]; delete idmaps[i];
idmaps.DeleteAll(); idmaps.DeleteAll();

4
libsrc/meshing/boundarylayer.cpp
View File

@ -535,9 +535,9 @@ namespace netgen
// Lock all the prism points so that the rest of the mesh can be // Lock all the prism points so that the rest of the mesh can be
// optimised without invalidating the entire mesh // optimised without invalidating the entire mesh
for (i = 1; i <= np; i++) for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++)
{ {
if(bndnodes.Test(i)) mesh.AddLockedPoint(i); if(bndnodes.Test(i)) mesh.AddLockedPoint(pi);
} }
// Now, actually pull back the old surface points to create // Now, actually pull back the old surface points to create

73
libsrc/meshing/delaunay.cpp
View File

@ -266,7 +266,6 @@ namespace netgen
Point3d tpmin, tpmax; Point3d tpmin, tpmax;
tettree.GetIntersecting (newp, newp, treesearch); tettree.GetIntersecting (newp, newp, treesearch);
double quot,minquot(1e20); double quot,minquot(1e20);
for (int j = 0; j < treesearch.Size(); j++) for (int j = 0; j < treesearch.Size(); j++)
@ -343,7 +342,6 @@ namespace netgen
int changed = 1; int changed = 1;
int nstarti = 1, starti; int nstarti = 1, starti;
while (changed) while (changed)
{ {
changed = 0; changed = 0;
@ -417,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 ( PointIndex (i3.I1())); const Point3d & p1 = mesh.Point ( i3.I1());
const Point3d & p2 = mesh.Point ( PointIndex (i3.I2())); const Point3d & p2 = mesh.Point ( i3.I2());
const Point3d & p3 = mesh.Point ( PointIndex (i3.I3())); const Point3d & p3 = mesh.Point ( i3.I3());
Vec3d v1(p1, p2); Vec3d v1(p1, p2);
@ -444,7 +442,6 @@ namespace netgen
} }
} }
} // while (changed) } // while (changed)
// (*testout) << "newels: " << endl; // (*testout) << "newels: " << endl;
Array<Element> newels; Array<Element> newels;
@ -510,7 +507,6 @@ namespace netgen
freelist.Append (celind); freelist.Append (celind);
} }
int hasclose = 0; int hasclose = 0;
for (int j = 1; j <= closesphere.GetArray().Size(); j++) for (int j = 1; j <= closesphere.GetArray().Size(); j++)
{ {
@ -605,9 +601,6 @@ namespace netgen
Array<DelaunayTet> & tempels, Array<DelaunayTet> & tempels,
int oldnp, DelaunayTet & startel, Point3d & pmin, Point3d & pmax) int oldnp, DelaunayTet & startel, Point3d & pmin, Point3d & pmax)
{ {
int i, j, k;
const Point<3> * pp[4];
Array<Point<3> > centers; Array<Point<3> > centers;
Array<double> radi2; Array<double> radi2;
@ -616,17 +609,17 @@ namespace netgen
// new: local box // new: local box
mesh.GetBox (pmax, pmin); // lower bound for pmax, upper for pmin mesh.GetBox (pmax, pmin); // lower bound for pmax, upper for pmin
for (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 (j = 0; j < face.GetNP(); j++) for (int j = 0; j < face.GetNP(); j++)
{ {
pmin.SetToMin (mesh.Point (face[j])); pmin.SetToMin (mesh.Point (face[j]));
pmax.SetToMax (mesh.Point (face[j])); pmax.SetToMax (mesh.Point (face[j]));
} }
} }
for (i = 0; i < mesh.LockedPoints().Size(); i++) for (int i = 0; i < mesh.LockedPoints().Size(); i++)
{ {
pmin.SetToMin (mesh.Point (mesh.LockedPoints()[i])); pmin.SetToMin (mesh.Point (mesh.LockedPoints()[i]));
pmax.SetToMax (mesh.Point (mesh.LockedPoints()[i])); pmax.SetToMax (mesh.Point (mesh.LockedPoints()[i]));
@ -661,18 +654,18 @@ namespace netgen
// flag points to use for Delaunay: // flag points to use for Delaunay:
BitArrayChar<PointIndex::BASE> usep(np); BitArrayChar<PointIndex::BASE> usep(np);
usep.Clear(); usep.Clear();
for (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 (j = 0; j < face.GetNP(); j++) for (int j = 0; j < face.GetNP(); j++)
usep.Set (face[j]); usep.Set (face[j]);
} }
for (i = oldnp + PointIndex::BASE; for (int i = oldnp + PointIndex::BASE;
i < np + PointIndex::BASE; i++) i < np + PointIndex::BASE; i++)
usep.Set (i); usep.Set (i);
for (i = 0; i < mesh.LockedPoints().Size(); i++) for (int i = 0; i < mesh.LockedPoints().Size(); i++)
usep.Set (mesh.LockedPoints()[i]); usep.Set (mesh.LockedPoints()[i]);
@ -697,7 +690,7 @@ namespace netgen
tpmin = tpmax = mesh.Point(startel[0]); tpmin = tpmax = mesh.Point(startel[0]);
for (k = 1; k < 4; k++) for (int k = 1; k < 4; k++)
{ {
tpmin.SetToMin (mesh.Point (startel[k])); tpmin.SetToMin (mesh.Point (startel[k]));
tpmax.SetToMax (mesh.Point (startel[k])); tpmax.SetToMax (mesh.Point (startel[k]));
@ -705,14 +698,11 @@ namespace netgen
tpmax = tpmax + 0.01 * (tpmax - tpmin); tpmax = tpmax + 0.01 * (tpmax - tpmin);
tettree.Insert (tpmin, tpmax, 1); tettree.Insert (tpmin, tpmax, 1);
Point<3> pc; Point<3> pc;
for (k = 0; k < 4; k++) const Point<3> * pp[4];
{ for (int k = 0; k < 4; k++)
pp[k] = &mesh.Point (startel[k]); pp[k] = &mesh.Point (startel[k]);
}
CalcSphereCenter (&pp[0], pc); CalcSphereCenter (&pp[0], pc);
centers.Append (pc); centers.Append (pc);
@ -722,50 +712,50 @@ namespace netgen
IndexSet insphere(mesh.GetNP()); IndexSet insphere(mesh.GetNP());
IndexSet closesphere(mesh.GetNP()); IndexSet closesphere(mesh.GetNP());
// "random" reordering of points (speeds a factor 3 - 5 !!!) // "random" reordering of points (speeds a factor 3 - 5 !!!)
Array<PointIndex, PointIndex::BASE, PointIndex> mixed(np);
Array<int> mixed(np);
int prims[] = { 11, 13, 17, 19, 23, 29, 31, 37 }; int prims[] = { 11, 13, 17, 19, 23, 29, 31, 37 };
int prim; int prim;
i = 0; {
int i = 0;
while (np % prims[i] == 0) i++; while (np % prims[i] == 0) i++;
prim = prims[i]; prim = prims[i];
}
for (i = 1; i <= np; i++) for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End()-4; pi++)
mixed.Elem(i) = (prim * i) % np + PointIndex::BASE; mixed[pi] = PointIndex ( (prim * pi) % np + PointIndex::BASE );
for (i = 1; i <= np; i++) for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End()-4; pi++)
{ {
if (i % 1000 == 0) if (pi % 1000 == 0)
{ {
if (i % 10000 == 0) if (pi % 10000 == 0)
PrintDot ('+'); PrintDot ('+');
else else
PrintDot ('.'); PrintDot ('.');
} }
multithread.percent = 100.0 * i / np; multithread.percent = 100.0 * pi / np;
if (multithread.terminate) if (multithread.terminate)
break; break;
PointIndex newpi = mixed.Get(i); PointIndex newpi = mixed[pi];
if (!usep.Test(newpi)) if (!usep.Test(newpi))
continue; continue;
cntp++; cntp++;
const Point3d & newp = mesh.Point(newpi); const MeshPoint & newp = mesh[newpi];
AddDelaunayPoint (newpi, newp, tempels, mesh, AddDelaunayPoint (newpi, newp, tempels, mesh,
tettree, meshnb, centers, radi2, tettree, meshnb, centers, radi2,
connected, treesearch, freelist, list, insphere, closesphere); connected, treesearch, freelist, list, insphere, closesphere);
} }
for (i = tempels.Size(); i >= 1; i--) for (int i = tempels.Size(); i >= 1; i--)
if (tempels.Get(i)[0] <= 0) if (tempels.Get(i)[0] <= 0)
tempels.DeleteElement (i); tempels.DeleteElement (i);
@ -819,7 +809,7 @@ namespace netgen
// improve delaunay - mesh by swapping !!!! // improve delaunay - mesh by swapping !!!!
Mesh tempmesh; Mesh tempmesh;
for (PointIndex pi = PointIndex::BASE; pi < mesh.GetNP()+PointIndex::BASE; pi++) for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++)
tempmesh.AddPoint (mesh[pi]); tempmesh.AddPoint (mesh[pi]);
for (int i = 1; i <= tempels.Size(); i++) for (int i = 1; i <= tempels.Size(); i++)
@ -874,8 +864,8 @@ namespace netgen
// for (i = mesh.GetNP() - 3; i <= mesh.GetNP(); i++) // for (i = mesh.GetNP() - 3; i <= mesh.GetNP(); i++)
// tempmesh.AddLockedPoint (i); // tempmesh.AddLockedPoint (i);
for (PointIndex pi = PointIndex::BASE; for (PointIndex pi = tempmesh.Points().Begin();
pi < tempmesh.GetNP() + PointIndex::BASE; pi++) pi < tempmesh.Points().End(); pi++)
tempmesh.AddLockedPoint (pi); tempmesh.AddLockedPoint (pi);
// tempmesh.PrintMemInfo(cout); // tempmesh.PrintMemInfo(cout);
@ -1288,8 +1278,7 @@ namespace netgen
INDEX_3_CLOSED_HASHTABLE<INDEX_2> faceht(100); INDEX_3_CLOSED_HASHTABLE<INDEX_2> faceht(100);
Element2d hel(TRIG); Element2d hel(TRIG);
for (PointIndex pi = PointIndex::BASE; for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++)
pi < mesh.GetNP()+PointIndex::BASE; pi++)
{ {
faceht.SetSize (4 * elsonpoint[pi].Size()); faceht.SetSize (4 * elsonpoint[pi].Size());
for (int ii = 0; ii < elsonpoint[pi].Size(); ii++) for (int ii = 0; ii < elsonpoint[pi].Size(); ii++)

4
libsrc/meshing/delaunay2d.cpp
View File

@ -94,7 +94,7 @@ namespace netgen
{ {
if (meshbox.IsIn (npoints.Get(i))) if (meshbox.IsIn (npoints.Get(i)))
{ {
int gpnum = mesh.AddPoint (npoints.Get(i)); PointIndex gpnum = mesh.AddPoint (npoints.Get(i));
adfront->AddPoint (npoints.Get(i), gpnum); adfront->AddPoint (npoints.Get(i), gpnum);
if (debugparam.slowchecks) if (debugparam.slowchecks)
@ -149,7 +149,7 @@ namespace netgen
{ {
if (meshbox.IsIn (npoints.Get(i))) if (meshbox.IsIn (npoints.Get(i)))
{ {
int gpnum = mesh.AddPoint (npoints.Get(i)); PointIndex gpnum = mesh.AddPoint (npoints.Get(i));
adfront->AddPoint (npoints.Get(i), gpnum); adfront->AddPoint (npoints.Get(i), gpnum);
} }
} }

2
libsrc/meshing/geomsearch.cpp
View File

@ -19,7 +19,7 @@ namespace netgen
} }
} }
void GeomSearch3d :: Init (Array <FrontPoint3,PointIndex::BASE> *pointsi, Array <FrontFace> *facesi) void GeomSearch3d :: Init (Array <FrontPoint3,PointIndex::BASE, PointIndex> *pointsi, Array <FrontFace> *facesi)
{ {
points = pointsi; points = pointsi;
faces = facesi; faces = facesi;

4
libsrc/meshing/geomsearch.hpp
View File

@ -22,7 +22,7 @@ public:
virtual ~GeomSearch3d(); virtual ~GeomSearch3d();
/// ///
void Init (Array <FrontPoint3,PointIndex::BASE> *pointsi, Array <FrontFace> *facesi); void Init (Array <FrontPoint3,PointIndex::BASE, PointIndex> *pointsi, Array <FrontFace> *facesi);
///get elements max extension ///get elements max extension
void ElemMaxExt(Point3d& minp, Point3d& maxp, const MiniElement2d& elem); void ElemMaxExt(Point3d& minp, Point3d& maxp, const MiniElement2d& elem);
@ -47,7 +47,7 @@ public:
private: private:
Array <FrontFace> *faces; // Pointers to Arrays in Adfront Array <FrontFace> *faces; // Pointers to Arrays in Adfront
Array <FrontPoint3,PointIndex::BASE> *points; Array <FrontPoint3,PointIndex::BASE, PointIndex> *points;
Array <Array <int>*> hashtable; Array <Array <int>*> hashtable;

3
libsrc/meshing/improve2.cpp
View File

@ -533,8 +533,7 @@ namespace netgen
Array<Vec<3>,PointIndex::BASE> normals(np); Array<Vec<3>,PointIndex::BASE> normals(np);
for (PointIndex pi = PointIndex::BASE; for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++)
pi < np + PointIndex::BASE; pi++)
{ {
if (elementsonnode[pi].Size()) if (elementsonnode[pi].Size())
{ {

2
libsrc/meshing/improve2gen.cpp
View File

@ -308,7 +308,7 @@ namespace netgen
if (el0.GetNP() != rel0.GetNP()) continue; if (el0.GetNP() != rel0.GetNP()) continue;
pmap = -1; pmap = PointIndex (-1);
for (int k = 0; k < el0.GetNP(); k++) for (int k = 0; k < el0.GetNP(); k++)
{ {

25
libsrc/meshing/meshclass.hpp
View File

@ -13,7 +13,6 @@
namespace netgen namespace netgen
{ {
enum resthtype { RESTRICTH_FACE, RESTRICTH_EDGE, enum resthtype { RESTRICTH_FACE, RESTRICTH_EDGE,
RESTRICTH_SURFACEELEMENT, RESTRICTH_POINT, RESTRICTH_SEGMENT }; RESTRICTH_SURFACEELEMENT, RESTRICTH_POINT, RESTRICTH_SEGMENT };
@ -25,16 +24,9 @@ namespace netgen
{ {
public: public:
typedef ::netgen::T_POINTS T_POINTS; typedef ::netgen::T_POINTS T_POINTS;
// typedef MoveableArray<MeshPoint,PointIndex::BASE> T_POINTS;
// typedef MoveableArray<Element> T_VOLELEMENTS;
// typedef MoveableArray<Element2d> T_SURFELEMENTS;
// typedef Array<MeshPoint,PointIndex::BASE> T_POINTS;
typedef Array<Element> T_VOLELEMENTS; typedef Array<Element> T_VOLELEMENTS;
typedef Array<Element2d> T_SURFELEMENTS; typedef Array<Element2d> T_SURFELEMENTS;
private: private:
/// point coordinates /// point coordinates
T_POINTS points; T_POINTS points;
@ -105,7 +97,7 @@ namespace netgen
/// geometric search tree for interval intersection search /// geometric search tree for interval intersection search
Box3dTree * elementsearchtree; Box3dTree * elementsearchtree;
/// time stamp for tree /// time stamp for tree
int elementsearchtreets; mutable int elementsearchtreets;
/// element -> face, element -> edge etc ... /// element -> face, element -> edge etc ...
class MeshTopology * topology; class MeshTopology * topology;
@ -216,12 +208,7 @@ namespace netgen
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);
/*
#ifdef PARALLEL
PointIndex AddPoint (const Point3d & p, bool aisghost, int layer = 1);
PointIndex AddPoint (const Point3d & p, bool aisghost, int layer, POINTTYPE type);
#endif
*/
int GetNP () const { return points.Size(); } int GetNP () const { return points.Size(); }
MeshPoint & Point(int i) { return points.Elem(i); } MeshPoint & Point(int i) { return points.Elem(i); }
@ -544,22 +531,22 @@ namespace netgen
void SetPointSearchStartElement(const int el) const {ps_startelement = el;} void SetPointSearchStartElement(const int el) const {ps_startelement = el;}
/// gives element of point, barycentric coordinates /// gives element of point, barycentric coordinates
int GetElementOfPoint (const Point3d & p, int GetElementOfPoint (const netgen::Point<3> & p,
double * lami, double * lami,
bool build_searchtree = 0, bool build_searchtree = 0,
const int index = -1, const int index = -1,
const bool allowindex = true) const; const bool allowindex = true) const;
int GetElementOfPoint (const Point3d & p, int GetElementOfPoint (const netgen::Point<3> & p,
double * lami, double * lami,
const Array<int> * const indices, const Array<int> * const indices,
bool build_searchtree = 0, bool build_searchtree = 0,
const bool allowindex = true) const; const bool allowindex = true) const;
int GetSurfaceElementOfPoint (const Point3d & p, int GetSurfaceElementOfPoint (const netgen::Point<3> & p,
double * lami, double * lami,
bool build_searchtree = 0, bool build_searchtree = 0,
const int index = -1, const int index = -1,
const bool allowindex = true) const; const bool allowindex = true) const;
int GetSurfaceElementOfPoint (const Point3d & p, int GetSurfaceElementOfPoint (const netgen::Point<3> & p,
double * lami, double * lami,
const Array<int> * const indices, const Array<int> * const indices,
bool build_searchtree = 0, bool build_searchtree = 0,

11
libsrc/meshing/meshfunc.cpp
View File

@ -130,8 +130,7 @@ namespace netgen
mpquad.check_impossible = qstep == 1; // for prisms only (air domain in trafo) mpquad.check_impossible = qstep == 1; // for prisms only (air domain in trafo)
for (PointIndex pi = PointIndex::BASE; for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
pi < mesh3d.GetNP()+PointIndex::BASE; pi++)
meshing.AddPoint (mesh3d[pi], pi); meshing.AddPoint (mesh3d[pi], pi);
mesh3d.GetIdentifications().GetPairs (0, connectednodes); mesh3d.GetIdentifications().GetPairs (0, connectednodes);
@ -174,11 +173,9 @@ namespace netgen
mesh3d.FindOpenElements(k); mesh3d.FindOpenElements(k);
for (PointIndex pi = PointIndex::BASE; for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
pi < mesh3d.GetNP()+PointIndex::BASE; pi++)
meshing.AddPoint (mesh3d[pi], pi); meshing.AddPoint (mesh3d[pi], pi);
for (int i = 1; i <= mesh3d.GetNOpenElements(); i++) for (int i = 1; i <= mesh3d.GetNOpenElements(); i++)
meshing.AddBoundaryElement (mesh3d.OpenElement(i)); meshing.AddBoundaryElement (mesh3d.OpenElement(i));
@ -217,9 +214,7 @@ namespace netgen
Array<int, PointIndex::BASE> glob2loc(mesh3d.GetNP()); Array<int, PointIndex::BASE> glob2loc(mesh3d.GetNP());
glob2loc = -1; glob2loc = -1;
for (PointIndex pi = PointIndex::BASE; for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
pi < mesh3d.GetNP()+PointIndex::BASE; pi++)
if (domain_bbox.IsIn (mesh3d[pi])) if (domain_bbox.IsIn (mesh3d[pi]))
glob2loc[pi] = glob2loc[pi] =
meshing.AddPoint (mesh3d[pi], pi); meshing.AddPoint (mesh3d[pi], pi);

3
libsrc/meshing/meshing2.cpp
View File

@ -221,7 +221,6 @@ namespace netgen
int z1, z2, oldnp(-1); int z1, z2, oldnp(-1);
bool found; bool found;
int rulenr(-1); int rulenr(-1);
int globind;
Point<3> p1, p2; Point<3> p1, p2;
const PointGeomInfo * blgeominfo1; const PointGeomInfo * blgeominfo1;
@ -1243,7 +1242,7 @@ namespace netgen
for (int i = oldnp+1; i <= locpoints.Size(); i++) for (int i = oldnp+1; i <= locpoints.Size(); i++)
{ {
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);
} }

12
libsrc/meshing/meshing3.cpp
View File

@ -190,7 +190,6 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
int found; int found;
referencetransform trans; referencetransform trans;
int rotind; int rotind;
INDEX globind;
Point3d inp; Point3d inp;
float err; float err;
@ -657,7 +656,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
for (i = oldnp+1; i <= locpoints.Size(); i++) for (i = oldnp+1; i <= locpoints.Size(); i++)
{ {
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);
} }
@ -747,7 +746,7 @@ void Meshing3 :: BlockFill (Mesh & mesh, double gh)
for (i = 1; i <= adfront->GetNP(); i++) for (i = 1; i <= adfront->GetNP(); i++)
{ {
const Point3d & p = adfront->GetPoint(i); const Point3d & p = adfront->GetPoint(PointIndex(i));
if (i == 1) if (i == 1)
{ {
xmin = xmax = p.X(); xmin = xmax = p.X();
@ -777,7 +776,8 @@ void Meshing3 :: BlockFill (Mesh & mesh, double gh)
PrintMessage (5, "n1 = ", n1, " n2 = ", n2, " n3 = ", n3); PrintMessage (5, "n1 = ", n1, " n2 = ", n2, " n3 = ", n3);
Array<blocktyp> inner(n); Array<blocktyp> inner(n);
Array<int> pointnr(n), frontpointnr(n); Array<PointIndex> pointnr(n);
Array<int> frontpointnr(n);
// initialize inner to 1 // initialize inner to 1
@ -1187,7 +1187,7 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
{ {
if (meshbox.IsIn (npoints.Get(i))) if (meshbox.IsIn (npoints.Get(i)))
{ {
int gpnum = mesh.AddPoint (npoints.Get(i)); PointIndex gpnum = mesh.AddPoint (npoints.Get(i));
adfront->AddPoint (npoints.Get(i), gpnum); adfront->AddPoint (npoints.Get(i), gpnum);
if (debugparam.slowchecks) if (debugparam.slowchecks)
@ -1255,7 +1255,7 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
{ {
if (meshbox.IsIn (npoints.Get(i))) if (meshbox.IsIn (npoints.Get(i)))
{ {
int gpnum = mesh.AddPoint (npoints.Get(i)); PointIndex gpnum = mesh.AddPoint (npoints.Get(i));
adfront->AddPoint (npoints.Get(i), gpnum); adfront->AddPoint (npoints.Get(i), gpnum);
} }
} }

23
libsrc/meshing/meshtype.hpp
View File

@ -120,11 +120,16 @@ namespace netgen
PointIndex () { ; } PointIndex () { ; }
PointIndex (int ai) : i(ai) { ; } PointIndex (int ai) : i(ai) { ; }
PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; } PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; }
PointIndex & operator= (int ai) { i = ai; return *this; } // PointIndex & operator= (int ai) { i = ai; return *this; }
operator int () const { return i; } operator int () const { return i; }
int GetInt () const { return i; } // int GetInt () const { return i; }
PointIndex operator++ (int) { int hi = i; i++; return PointIndex(hi); } // PointIndex operator+ (int i2) { return PointIndex (i+i2); }
PointIndex operator-- (int) { int hi = i; i--; return PointIndex(hi); } // PointIndex operator++ (int) { int hi = i; i++; return PointIndex(hi); }
// PointIndex operator-- (int) { int hi = i; i--; return PointIndex(hi); }
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; }
#ifdef BASE0 #ifdef BASE0
enum { BASE = 0 }; enum { BASE = 0 };
@ -140,7 +145,7 @@ namespace netgen
inline ostream & operator<< (ostream & ost, const PointIndex & pi) inline ostream & operator<< (ostream & ost, const PointIndex & pi)
{ {
return (ost << pi.GetInt()); return (ost << int(pi));
} }
@ -275,7 +280,7 @@ namespace netgen
typedef Array<MeshPoint, PointIndex::BASE> T_POINTS; typedef Array<MeshPoint, PointIndex::BASE, PointIndex> T_POINTS;
@ -371,6 +376,9 @@ namespace netgen
/// ///
const PointIndex & operator[] (int i) const { return pnum[i]; } const PointIndex & operator[] (int i) const { return pnum[i]; }
FlatArray<const PointIndex> PNums () const
{ return FlatArray<const PointIndex> (np, &pnum[0]); }
/// ///
PointIndex & PNum (int i) { return pnum[i-1]; } PointIndex & PNum (int i) { return pnum[i-1]; }
/// ///
@ -618,6 +626,9 @@ namespace netgen
/// ///
const PointIndex & operator[] (int i) const { return pnum[i]; } const PointIndex & operator[] (int i) const { return pnum[i]; }
FlatArray<const PointIndex> PNums () const
{ return FlatArray<const PointIndex> (np, &pnum[0]); }
/// ///
PointIndex & PNum (int i) { return pnum[i-1]; } PointIndex & PNum (int i) { return pnum[i-1]; }
/// ///

22
libsrc/meshing/refine.cpp
View File

@ -18,7 +18,7 @@ namespace netgen
mesh.SetNP(mesh.GetNV()); mesh.SetNP(mesh.GetNV());
INDEX_2_HASHTABLE<int> between(mesh.GetNP() + 5); INDEX_2_HASHTABLE<PointIndex> between(mesh.GetNP() + 5);
int oldne, oldns, oldnf; int oldne, oldns, oldnf;
@ -87,7 +87,7 @@ namespace netgen
case TRIG: case TRIG:
case TRIG6: case TRIG6:
{ {
ArrayMem<int,6> pnums(6); ArrayMem<PointIndex,6> pnums(6);
ArrayMem<PointGeomInfo,6> pgis(6); ArrayMem<PointGeomInfo,6> pgis(6);
static int betw[3][3] = static int betw[3][3] =
@ -162,7 +162,7 @@ namespace netgen
case QUAD6: case QUAD6:
case QUAD8: case QUAD8:
{ {
ArrayMem<int,9> pnums(9); ArrayMem<PointIndex,9> pnums(9);
ArrayMem<PointGeomInfo,9> pgis(9); ArrayMem<PointGeomInfo,9> pgis(9);
static int betw[5][3] = static int betw[5][3] =
@ -253,7 +253,7 @@ namespace netgen
case TET: case TET:
case TET10: case TET10:
{ {
ArrayMem<int,10> pnums(10); ArrayMem<PointIndex,10> pnums(10);
static int betw[6][3] = static int betw[6][3] =
{ { 1, 2, 5 }, { { 1, 2, 5 },
{ 1, 3, 6 }, { 1, 3, 6 },
@ -334,7 +334,7 @@ namespace netgen
} }
case HEX: case HEX:
{ {
ArrayMem<int,27> pnums(27); ArrayMem<PointIndex,27> pnums(27);
static int betw[13][3] = static int betw[13][3] =
{ { 1, 2, 9 }, { { 1, 2, 9 },
{ 3, 4, 10 }, { 3, 4, 10 },
@ -384,7 +384,7 @@ namespace netgen
{ 17, 20, 27 }, { 17, 20, 27 },
}; };
pnums = -1; pnums = PointIndex(-1);
for (j = 1; j <= 8; j++) for (j = 1; j <= 8; j++)
pnums.Elem(j) = el.PNum(j); pnums.Elem(j) = el.PNum(j);
@ -460,7 +460,7 @@ namespace netgen
} }
case PRISM: case PRISM:
{ {
ArrayMem<int,18> pnums(18); ArrayMem<PointIndex,18> pnums(18);
static int betw[9][3] = static int betw[9][3] =
{ { 3, 1, 7 }, { { 3, 1, 7 },
{ 1, 2, 8 }, { 1, 2, 8 },
@ -494,7 +494,7 @@ namespace netgen
}; };
//int elrev = el.flags.reverse; //int elrev = el.flags.reverse;
pnums = -1; pnums = PointIndex(-1);
for (j = 1; j <= 6; j++) for (j = 1; j <= 6; j++)
pnums.Elem(j) = el.PNum(j); pnums.Elem(j) = el.PNum(j);
@ -594,14 +594,14 @@ namespace netgen
for (int k = 1; k <= between.GetBagSize(j); k++) for (int k = 1; k <= between.GetBagSize(j); k++)
{ {
INDEX_2 i2; INDEX_2 i2;
int newpi; PointIndex newpi;
between.GetData (j, k, i2, newpi); between.GetData (j, k, i2, newpi);
INDEX_2 oi2(identmap.Get(i2.I1()), INDEX_2 oi2(identmap.Get(i2.I1()),
identmap.Get(i2.I2())); identmap.Get(i2.I2()));
oi2.Sort(); oi2.Sort();
if (between.Used (oi2)) if (between.Used (oi2))
{ {
int onewpi = between.Get(oi2); PointIndex onewpi = between.Get(oi2);
mesh.GetIdentifications().Add (newpi, onewpi, i); mesh.GetIdentifications().Add (newpi, onewpi, i);
} }
} }
@ -638,7 +638,7 @@ namespace netgen
for (int j = 1; j <= between.GetBagSize(i); j++) for (int j = 1; j <= between.GetBagSize(i); j++)
{ {
INDEX_2 parent; INDEX_2 parent;
int child; PointIndex child;
between.GetData (i, j, parent, child); between.GetData (i, j, parent, child);
can.Elem(child) = Center (can.Elem(parent.I1()), can.Elem(child) = Center (can.Elem(parent.I1()),
can.Elem(parent.I2())); can.Elem(parent.I2()));

10
libsrc/meshing/secondorder.cpp
View File

@ -18,7 +18,7 @@ namespace netgen
mesh.ComputeNVertices(); mesh.ComputeNVertices();
mesh.SetNP(mesh.GetNV()); mesh.SetNP(mesh.GetNV());
INDEX_2_HASHTABLE<int> between(mesh.GetNP() + 5); INDEX_2_HASHTABLE<PointIndex> between(mesh.GetNP() + 5);
bool thinlayers = 0; bool thinlayers = 0;
@ -235,18 +235,18 @@ namespace netgen
Array<int,PointIndex::BASE> identmap; Array<int,PointIndex::BASE> identmap;
mesh.GetIdentifications().GetMap (i, identmap); mesh.GetIdentifications().GetMap (i, identmap);
for (INDEX_2_HASHTABLE<int>::Iterator it = between.Begin(); for (INDEX_2_HASHTABLE<PointIndex>::Iterator it = between.Begin();
it != between.End(); it++) it != between.End(); it++)
{ {
INDEX_2 i2; INDEX_2 i2;
int newpi; PointIndex newpi;
between.GetData (it, i2, newpi); between.GetData (it, i2, newpi);
INDEX_2 oi2(identmap.Get(i2.I1()), INDEX_2 oi2(identmap.Get(i2.I1()),
identmap.Get(i2.I2())); identmap.Get(i2.I2()));
oi2.Sort(); oi2.Sort();
if (between.Used (oi2)) if (between.Used (oi2))
{ {
int onewpi = between.Get(oi2); PointIndex onewpi = between.Get(oi2);
mesh.GetIdentifications().Add (newpi, onewpi, i); mesh.GetIdentifications().Add (newpi, onewpi, i);
} }
} }
@ -288,7 +288,7 @@ namespace netgen
} }
*/ */
for (INDEX_2_HASHTABLE<int>::Iterator it = between.Begin(); for (INDEX_2_HASHTABLE<PointIndex>::Iterator it = between.Begin();
it != between.End(); it++) it != between.End(); it++)
{ {
mesh.mlbetweennodes[between.GetData (it)] = between.GetHash(it); mesh.mlbetweennodes[between.GetData (it)] = between.GetHash(it);

3
libsrc/meshing/smoothing2.5.cpp
View File

@ -141,8 +141,7 @@ namespace netgen
Array<SurfaceElementIndex> locelements(0); Array<SurfaceElementIndex> locelements(0);
Array<int> locrots(0); Array<int> locrots(0);
for (PointIndex pi = PointIndex::BASE; for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++)
pi < mesh.GetNP()+PointIndex::BASE; pi++)
{ {
if (mesh[pi].Type() != SURFACEPOINT) if (mesh[pi].Type() != SURFACEPOINT)
continue; continue;

53
libsrc/meshing/smoothing3.cpp
View File

@ -106,9 +106,9 @@ namespace netgen
{ {
const INDEX_3 & el = faces[j]; const INDEX_3 & el = faces[j];
double bad = CalcTetBadness (points[el.I1()], double bad = CalcTetBadness (points[PointIndex (el.I1())],
points[el.I3()], points[PointIndex (el.I3())],
points[el.I2()], points[PointIndex (el.I2())],
pp, 0, mp); pp, 0, mp);
badness += bad; badness += bad;
} }
@ -202,9 +202,9 @@ namespace netgen
for (int i = 1; i <= nf; i++) for (int i = 1; i <= nf; i++)
{ {
const Point3d & p1 = points[faces.Get(i).I1()]; const Point3d & p1 = points[PointIndex(faces.Get(i).I1())];
const Point3d & p2 = points[faces.Get(i).I2()]; const Point3d & p2 = points[PointIndex(faces.Get(i).I2())];
const Point3d & p3 = points[faces.Get(i).I3()]; const Point3d & p3 = points[PointIndex(faces.Get(i).I3())];
Vec3d v1 (p1, p2); Vec3d v1 (p1, p2);
Vec3d v2 (p1, p3); Vec3d v2 (p1, p3);
Vec3d n; Vec3d n;
@ -500,7 +500,7 @@ namespace netgen
int ne = elementsonpoint[actpind].Size(); int ne = elementsonpoint[actpind].Size();
int i, j; int i, j;
int pi1, pi2, pi3; PointIndex pi1, pi2, pi3;
m.SetSize (ne, 4); m.SetSize (ne, 4);
@ -1383,7 +1383,7 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
bad1 += hbad; bad1 += hbad;
} }
for (PointIndex i = PointIndex::BASE; i < np+PointIndex::BASE; i++) for (int i = perrs.Begin(); i < perrs.End(); i++)
if (perrs[i] > badmax) if (perrs[i] > badmax)
badmax = perrs[i]; badmax = perrs[i];
badmax = 0; badmax = 0;
@ -1456,14 +1456,13 @@ 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 i = PointIndex::BASE; for (PointIndex pi = points.Begin(); pi < points.End(); pi++)
i < points.Size()+PointIndex::BASE; i++) if ( (*this)[pi].Type() == INNERPOINT && perrs[pi] > 0.01 * badmax)
if ( (*this)[i].Type() == INNERPOINT && perrs[i] > 0.01 * badmax)
{ {
if (multithread.terminate) if (multithread.terminate)
throw NgException ("Meshing stopped"); throw NgException ("Meshing stopped");
multithread.percent = 100.0 * (i+1-PointIndex::BASE) / points.Size(); multithread.percent = 100.0 * (pi+1-PointIndex::BASE) / points.Size();
/* /*
if (points.Size() < 1000) if (points.Size() < 1000)
PrintDot (); PrintDot ();
@ -1471,14 +1470,14 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
if ( (i+1-PointIndex::BASE) % 10 == 0) if ( (i+1-PointIndex::BASE) % 10 == 0)
PrintDot ('+'); PrintDot ('+');
*/ */
if ( (i+1-PointIndex::BASE) % printmod == 0) PrintDot (printdot); if ( (pi+1-PointIndex::BASE) % printmod == 0) PrintDot (printdot);
double lh = pointh[i]; double lh = pointh[pi];
pf->SetLocalH (lh); pf->SetLocalH (lh);
par.typx = lh; par.typx = lh;
freeminf.SetPoint (points[i]); freeminf.SetPoint (points[pi]);
pf->SetPointIndex (i); pf->SetPointIndex (pi);
x = 0; x = 0;
int pok; int pok;
@ -1488,8 +1487,8 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
{ {
pok = pf->MovePointToInner (); pok = pf->MovePointToInner ();
freeminf.SetPoint (points[i]); freeminf.SetPoint (points[pi]);
pf->SetPointIndex (i); pf->SetPointIndex (pi);
} }
if (pok) if (pok)
@ -1497,9 +1496,9 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
//*testout << "start BFGS, pok" << endl; //*testout << "start BFGS, pok" << endl;
BFGS (x, freeminf, par); BFGS (x, freeminf, par);
//*testout << "BFGS complete, pok" << endl; //*testout << "BFGS complete, pok" << endl;
points[i](0) += x(0); points[pi](0) += x(0);
points[i](1) += x(1); points[pi](1) += x(1);
points[i](2) += x(2); points[pi](2) += x(2);
} }
} }
PrintDot ('\n'); PrintDot ('\n');
@ -1581,9 +1580,9 @@ void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp,
const char * savetask = multithread.task; const char * savetask = multithread.task;
multithread.task = "Smooth Mesh Jacobian"; multithread.task = "Smooth Mesh Jacobian";
for (i = 1; i <= points.Size(); i++) for (PointIndex pi = points.Begin(); i < points.End(); pi++)
{ {
if ((*this)[PointIndex(i)].Type() != INNERPOINT) if ((*this)[pi].Type() != INNERPOINT)
continue; continue;
if(usepoint && !usepoint->Test(i)) if(usepoint && !usepoint->Test(i))
@ -1613,7 +1612,7 @@ void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp,
double lh = pointh[i]; double lh = pointh[i];
par.typx = lh; par.typx = lh;
pf.SetPointIndex (i); pf.SetPointIndex (pi);
x = 0; x = 0;
int pok = (pf.Func (x) < 1e10); int pok = (pf.Func (x) < 1e10);
@ -1735,7 +1734,7 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp,
const char * savetask = multithread.task; const char * savetask = multithread.task;
multithread.task = "Smooth Mesh Jacobian"; multithread.task = "Smooth Mesh Jacobian";
for (i = 1; i <= points.Size(); i++) for (PointIndex pi = points.Begin(); pi <= points.End(); pi++)
if ( usepoint.Test(i) ) if ( usepoint.Test(i) )
{ {
//(*testout) << "improvejac, p = " << i << endl; //(*testout) << "improvejac, p = " << i << endl;
@ -1762,9 +1761,9 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp,
double lh = pointh[i];//GetH(points.Get(i)); double lh = pointh[i];//GetH(points.Get(i));
par.typx = lh; par.typx = lh;
pf.SetPointIndex (i); pf.SetPointIndex (pi);
int brother = -1; PointIndex brother (-1);
if(usesum) if(usesum)
{ {
for(j=0; brother == -1 && j<used_idmaps->Size(); j++) for(j=0; brother == -1 && j<used_idmaps->Size(); j++)

4
libsrc/meshing/specials.cpp
View File

@ -126,7 +126,7 @@ void CutOffAndCombine (Mesh & mesh, const Mesh & othermesh)
for (j = 1; j <= 3; j++) for (j = 1; j <= 3; j++)
locked.Clear (mesh.OpenElement(i).PNum(j)); locked.Clear (mesh.OpenElement(i).PNum(j));
for (i = 1; i <= locked.Size(); i++) for (PointIndex i (1); i <= locked.Size(); i++)
if (locked.Test(i)) if (locked.Test(i))
{ {
mesh.AddLockedPoint (i); mesh.AddLockedPoint (i);
@ -135,7 +135,7 @@ void CutOffAndCombine (Mesh & mesh, const Mesh & othermesh)
Array<int> pmat(onp); Array<PointIndex> pmat(onp);
for (i = 1; i <= onp; i++) for (i = 1; i <= onp; i++)
pmat.Elem(i) = mesh.AddPoint (othermesh.Point(i)); pmat.Elem(i) = mesh.AddPoint (othermesh.Point(i));

10
libsrc/meshing/topology.cpp
View File

@ -48,16 +48,6 @@ namespace netgen
vert2surfelement = 0; vert2surfelement = 0;
vert2segment = 0; vert2segment = 0;
timestamp = -1; timestamp = -1;
edge2vert.SetName ("edge2vert");
face2vert.SetName ("face2vert");
edges.SetName ("el2edge");
faces.SetName ("el2face");
surfedges.SetName ("surfel2edge");
segedges.SetName ("segment2edge");
surffaces.SetName ("surfel2face");
surf2volelement.SetName ("surfel2el");
face2surfel.SetName ("face2surfel");
} }
MeshTopology :: ~MeshTopology () MeshTopology :: ~MeshTopology ()

18
libsrc/meshing/topology.hpp
View File

@ -19,15 +19,15 @@ class MeshTopology
bool buildedges; bool buildedges;
bool buildfaces; bool buildfaces;
MoveableArray<INDEX_2> edge2vert; Array<INDEX_2> edge2vert;
MoveableArray<INDEX_4> face2vert; Array<INDEX_4> face2vert;
MoveableArray<int[12]> edges; Array<int[12]> edges;
MoveableArray<int[6]> faces; Array<int[6]> faces;
MoveableArray<int[4]> surfedges; Array<int[4]> surfedges;
MoveableArray<int> segedges; Array<int> segedges;
MoveableArray<int> surffaces; Array<int> surffaces;
MoveableArray<INDEX_2> surf2volelement; Array<INDEX_2> surf2volelement;
MoveableArray<int> face2surfel; Array<int> face2surfel;
TABLE<ElementIndex,PointIndex::BASE> *vert2element; TABLE<ElementIndex,PointIndex::BASE> *vert2element;
TABLE<int,PointIndex::BASE> *vert2surfelement; TABLE<int,PointIndex::BASE> *vert2surfelement;
TABLE<int,PointIndex::BASE> *vert2segment; TABLE<int,PointIndex::BASE> *vert2segment;

10
libsrc/meshing/zrefine.cpp
View File

@ -147,11 +147,11 @@ namespace netgen
for (j = 0; j <= 1; j++) for (j = 0; j <= 1; j++)
{ {
int pi1 = el.PNum( (j+0) % 4 + 1); PointIndex pi1 = el.PNum( (j+0) % 4 + 1);
int pi2 = el.PNum( (j+1) % 4 + 1); PointIndex pi2 = el.PNum( (j+1) % 4 + 1);
int pi3 = el.PNum( (j+2) % 4 + 1); PointIndex pi3 = el.PNum( (j+2) % 4 + 1);
int pi4 = el.PNum( (j+3) % 4 + 1); PointIndex pi4 = el.PNum( (j+3) % 4 + 1);
int pi5 = el.PNum(5); PointIndex pi5 = el.PNum(5);
INDEX_2 edge1(pi1, pi4); INDEX_2 edge1(pi1, pi4);
INDEX_2 edge2(pi2, pi3); INDEX_2 edge2(pi2, pi3);