mirror of
https://github.com/NGSolve/netgen.git
synced 2024-12-25 05:20:34 +05:00
more use of PointIndex
This commit is contained in:
Joachim Schöberl
parent
4c32c2ac25
commit
bdbc415589
@ -484,9 +484,9 @@ int AdFront3 :: SelectBaseElement ()
|
||||
|
||||
|
||||
int AdFront3 :: GetLocals (int fstind,
|
||||
Array<Point3d > & locpoints,
|
||||
Array<Point3d, PointIndex::BASE> & locpoints,
|
||||
Array<MiniElement2d> & locfaces, // local index
|
||||
Array<PointIndex> & pindex,
|
||||
Array<PointIndex, PointIndex::BASE> & pindex,
|
||||
Array<INDEX> & findex,
|
||||
INDEX_2_HASHTABLE<int> & getconnectedpairs,
|
||||
float xh,
|
||||
@ -600,12 +600,13 @@ int AdFront3 :: GetLocals (int fstind,
|
||||
if (invpindex[pi] == -1)
|
||||
{
|
||||
pindex.Append (pi);
|
||||
invpindex[pi] = pindex.Size(); // -1+PointIndex::BASE;
|
||||
locfaces.Elem(i).PNum(j) = locpoints.Append (points[pi].P());
|
||||
}
|
||||
else
|
||||
locfaces.Elem(i).PNum(j) = invpindex[pi];
|
||||
|
||||
locpoints.Append (points[pi].P());
|
||||
invpindex[pi] = pindex.Size()-1+PointIndex::BASE;
|
||||
}
|
||||
// locfaces.Elem(i).PNum(j) = locpoints.Append (points[pi].P());
|
||||
// }
|
||||
// else
|
||||
locfaces.Elem(i).PNum(j) = invpindex[pi];
|
||||
}
|
||||
}
|
||||
|
||||
@ -655,9 +656,9 @@ int AdFront3 :: GetLocals (int fstind,
|
||||
|
||||
// returns all points connected with fi
|
||||
void AdFront3 :: GetGroup (int fi,
|
||||
Array<MeshPoint> & grouppoints,
|
||||
Array<MeshPoint, PointIndex::BASE> & grouppoints,
|
||||
Array<MiniElement2d> & groupelements,
|
||||
Array<PointIndex> & pindex,
|
||||
Array<PointIndex, PointIndex::BASE> & pindex,
|
||||
Array<INDEX> & findex)
|
||||
{
|
||||
// static Array<char> pingroup;
|
||||
|
||||
@ -260,9 +260,9 @@ public:
|
||||
|
||||
///
|
||||
int GetLocals (int baseelement,
|
||||
Array<Point3d > & locpoints,
|
||||
Array<Point3d, PointIndex::BASE> & locpoints,
|
||||
Array<MiniElement2d> & locfaces, // local index
|
||||
Array<PointIndex> & pindex,
|
||||
Array<PointIndex, PointIndex::BASE> & pindex,
|
||||
Array<INDEX> & findex,
|
||||
INDEX_2_HASHTABLE<int> & connectedpairs,
|
||||
float xh,
|
||||
@ -271,9 +271,9 @@ public:
|
||||
|
||||
///
|
||||
void GetGroup (int fi,
|
||||
Array<MeshPoint> & grouppoints,
|
||||
Array<MeshPoint, PointIndex::BASE> & grouppoints,
|
||||
Array<MiniElement2d> & groupelements,
|
||||
Array<PointIndex> & pindex,
|
||||
Array<PointIndex, PointIndex::BASE> & pindex,
|
||||
Array<INDEX> & findex);
|
||||
|
||||
///
|
||||
|
||||
@ -590,7 +590,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
|
||||
void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
const BitArray * working_elements)
|
||||
{
|
||||
PointIndex pi1(0), pi2(0), pi3(0), pi4(0), pi5(0), pi6(0);
|
||||
PointIndex pi3(0), pi4(0), pi5(0), pi6(0);
|
||||
int cnt = 0;
|
||||
|
||||
Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET);
|
||||
@ -636,8 +636,12 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
|
||||
// find elements on node
|
||||
for (ElementIndex ei = 0; ei < ne; ei++)
|
||||
for (PointIndex pi : mesh[ei].PNums())
|
||||
elementsonnode.Add (pi, ei);
|
||||
/*
|
||||
for (int j = 0; j < mesh[ei].GetNP(); j++)
|
||||
elementsonnode.Add (mesh[ei][j], ei);
|
||||
*/
|
||||
|
||||
|
||||
// INDEX_2_HASHTABLE<int> edgeused(2 * ne + 5);
|
||||
@ -648,7 +652,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
if (multithread.terminate)
|
||||
break;
|
||||
|
||||
if(mesh.GetDimension()==3 && mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex())
|
||||
if (mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex())
|
||||
continue;
|
||||
|
||||
multithread.percent = 100.0 * (ei+1) / ne;
|
||||
@ -678,18 +682,14 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
const Element & elemi = mesh[ei];
|
||||
if (elemi.IsDeleted()) continue;
|
||||
|
||||
|
||||
// (*testout) << "check element " << elemi << endl;
|
||||
|
||||
int mattyp = elemi.GetIndex();
|
||||
|
||||
static const int tetedges[6][2] =
|
||||
{ { 0, 1 }, { 0, 2 }, { 0, 3 },
|
||||
{ 1, 2 }, { 1, 3 }, { 2, 3 } };
|
||||
|
||||
pi1 = elemi[tetedges[j][0]];
|
||||
pi2 = elemi[tetedges[j][1]];
|
||||
|
||||
PointIndex pi1 = elemi[tetedges[j][0]];
|
||||
PointIndex pi2 = elemi[tetedges[j][1]];
|
||||
|
||||
if (pi2 < pi1) Swap (pi1, pi2);
|
||||
|
||||
@ -718,21 +718,22 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
|
||||
if (has1 && has2)
|
||||
{ // only once
|
||||
for (int l = 0; l < hasbothpoints.Size(); l++)
|
||||
if (hasbothpoints[l] == elnr)
|
||||
has1 = 0;
|
||||
if (hasbothpoints.Contains (elnr))
|
||||
has1 = false;
|
||||
|
||||
if (has1)
|
||||
hasbothpoints.Append (elnr);
|
||||
{
|
||||
hasbothpoints.Append (elnr);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
bool puretet = 1;
|
||||
for (int k = 0; k < hasbothpoints.Size(); k++)
|
||||
if (mesh[hasbothpoints[k]].GetType () != TET)
|
||||
puretet = 0;
|
||||
if (!puretet)
|
||||
continue;
|
||||
bool puretet = true;
|
||||
for (ElementIndex ei : hasbothpoints)
|
||||
if (mesh[ei].GetType () != TET)
|
||||
puretet = false;
|
||||
|
||||
if (!puretet) continue;
|
||||
|
||||
int nsuround = hasbothpoints.Size();
|
||||
|
||||
@ -763,10 +764,10 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
for (int k = 0; k < 3; k++) // JS, 201212
|
||||
{
|
||||
const Element & elemk = mesh[hasbothpoints[k]];
|
||||
bool has1 = 0;
|
||||
bool has1 = false;
|
||||
for (int l = 0; l < 4; l++)
|
||||
if (elemk[l] == pi4)
|
||||
has1 = 1;
|
||||
has1 = true;
|
||||
if (has1)
|
||||
{
|
||||
for (int l = 0; l < 4; l++)
|
||||
@ -775,11 +776,10 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
}
|
||||
}
|
||||
|
||||
if(pi5 == 0)
|
||||
if (!pi5.IsValid())
|
||||
throw NgException("Illegal state observed in SwapImprove");
|
||||
|
||||
|
||||
|
||||
el32[0] = pi1;
|
||||
el32[1] = pi2;
|
||||
el32[2] = pi4;
|
||||
@ -882,7 +882,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
mesh[hasbothpoints[0]] = el21;
|
||||
mesh[hasbothpoints[1]] = el22;
|
||||
for (int l = 0; l < 4; l++)
|
||||
mesh[hasbothpoints[2]][l] = 0;
|
||||
mesh[hasbothpoints[2]][l].Invalidate();
|
||||
mesh[hasbothpoints[2]].Delete();
|
||||
|
||||
for (int k = 0; k < 2; k++)
|
||||
@ -912,14 +912,11 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
el1[3] = pi4;
|
||||
}
|
||||
|
||||
pi5 = 0;
|
||||
pi5.Invalidate();
|
||||
for (int k = 0; k < 4; k++)
|
||||
{
|
||||
const Element & elem = mesh[hasbothpoints[k]];
|
||||
bool has1 = 0;
|
||||
for (int l = 0; l < 4; l++)
|
||||
if (elem[l] == pi4)
|
||||
has1 = 1;
|
||||
bool has1 = elem.PNums().Contains(pi4);
|
||||
if (has1)
|
||||
{
|
||||
for (int l = 0; l < 4; l++)
|
||||
@ -928,14 +925,11 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
}
|
||||
}
|
||||
|
||||
pi6 = 0;
|
||||
pi6.Invalidate();
|
||||
for (int k = 0; k < 4; k++)
|
||||
{
|
||||
const Element & elem = mesh[hasbothpoints[k]];
|
||||
bool has1 = 0;
|
||||
for (int l = 0; l < 4; l++)
|
||||
if (elem[l] == pi3)
|
||||
has1 = 1;
|
||||
bool has1 = elem.PNums().Contains(pi3);
|
||||
if (has1)
|
||||
{
|
||||
for (int l = 0; l < 4; l++)
|
||||
@ -1202,7 +1196,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
Element hel(TET);
|
||||
|
||||
ArrayMem<PointIndex, 50> suroundpts(nsuround);
|
||||
ArrayMem<char, 50> tetused(nsuround);
|
||||
ArrayMem<bool, 50> tetused(nsuround);
|
||||
|
||||
Element & elem = mesh[hasbothpoints[0]];
|
||||
|
||||
@ -1228,33 +1222,32 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
|
||||
|
||||
// suroundpts.SetSize (nsuround);
|
||||
suroundpts = -17;
|
||||
suroundpts[0] = pi3;
|
||||
suroundpts[1] = pi4;
|
||||
|
||||
tetused = 0;
|
||||
tetused[0] = 1;
|
||||
tetused = false;
|
||||
tetused[0] = true;
|
||||
|
||||
for (int l = 2; l < nsuround; l++)
|
||||
{
|
||||
int oldpi = suroundpts[l-1];
|
||||
int newpi = 0;
|
||||
PointIndex oldpi = suroundpts[l-1];
|
||||
PointIndex newpi;
|
||||
newpi.Invalidate();
|
||||
|
||||
|
||||
for (int k = 0; k < nsuround && !newpi; k++)
|
||||
for (int k = 0; k < nsuround && !newpi.IsValid(); k++)
|
||||
if (!tetused[k])
|
||||
{
|
||||
const Element & nel = mesh[hasbothpoints[k]];
|
||||
|
||||
for (int k2 = 0; k2 < 4 && !newpi; k2++)
|
||||
for (int k2 = 0; k2 < 4 && !newpi.IsValid(); k2++)
|
||||
if (nel[k2] == oldpi)
|
||||
{
|
||||
newpi =
|
||||
nel[0] + nel[1] + nel[2] + nel[3]
|
||||
- pi1 - pi2 - oldpi;
|
||||
|
||||
tetused[k] = 1;
|
||||
tetused[k] = true;
|
||||
suroundpts[l] = newpi;
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -1406,8 +1399,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
|
||||
*/
|
||||
rel.Delete();
|
||||
for (int k1 = 0; k1 < 4; k1++)
|
||||
rel[k1] = 0;
|
||||
|
||||
rel[k1].Invalidate();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -176,14 +176,14 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
|
||||
NgProfiler::RegionTimer reg (meshing3_timer);
|
||||
|
||||
|
||||
Array<Point3d > locpoints; // local points
|
||||
Array<MiniElement2d> locfaces; // local faces
|
||||
Array<PointIndex> pindex; // mapping from local to front point numbering
|
||||
Array<int> allowpoint; // point is allowd ?
|
||||
Array<INDEX> findex; // mapping from local to front face numbering
|
||||
//INDEX_2_HASHTABLE<int> connectedpairs(100); // connecgted pairs for prism meshing
|
||||
Array<Point3d, PointIndex::BASE> locpoints; // local points
|
||||
Array<MiniElement2d> locfaces; // local faces
|
||||
Array<PointIndex, PointIndex::BASE> pindex; // mapping from local to front point numbering
|
||||
Array<int, PointIndex::BASE> allowpoint; // point is allowd ?
|
||||
Array<INDEX> findex; // mapping from local to front face numbering
|
||||
//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<Element> locelements; // new generated elements
|
||||
|
||||
@ -206,9 +206,9 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
|
||||
|
||||
|
||||
// for star-shaped domain meshing
|
||||
Array<MeshPoint> grouppoints;
|
||||
Array<MeshPoint, PointIndex::BASE> grouppoints;
|
||||
Array<MiniElement2d> groupfaces;
|
||||
Array<PointIndex> grouppindex;
|
||||
Array<PointIndex, PointIndex::BASE> grouppindex;
|
||||
Array<INDEX> groupfindex;
|
||||
|
||||
|
||||
@ -269,9 +269,9 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
|
||||
}
|
||||
|
||||
const MiniElement2d & bel = adfront->GetFace (baseelem);
|
||||
const Point3d & p1 = adfront->GetPoint (bel.PNum(1));
|
||||
const Point3d & p2 = adfront->GetPoint (bel.PNum(2));
|
||||
const Point3d & p3 = adfront->GetPoint (bel.PNum(3));
|
||||
const Point3d & p1 = adfront->GetPoint (bel[0]);
|
||||
const Point3d & p2 = adfront->GetPoint (bel[1]);
|
||||
const Point3d & p3 = adfront->GetPoint (bel[2]);
|
||||
|
||||
// (*testout) << endl << "base = " << bel << endl;
|
||||
|
||||
@ -302,9 +302,6 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
|
||||
|
||||
// (*testout) << "locfaces = " << endl << locfaces << endl;
|
||||
|
||||
int pi1 = pindex.Get(locfaces[0].PNum(1));
|
||||
int pi2 = pindex.Get(locfaces[0].PNum(2));
|
||||
int pi3 = pindex.Get(locfaces[0].PNum(3));
|
||||
|
||||
//loktestmode = 1;
|
||||
testmode = loktestmode; //changed
|
||||
@ -316,6 +313,9 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
|
||||
if (loktestmode)
|
||||
{
|
||||
(*testout) << "baseel = " << baseelem << ", ind = " << findex.Get(1) << endl;
|
||||
int pi1 = pindex[locfaces[0].PNum(1)];
|
||||
int pi2 = pindex[locfaces[0].PNum(2)];
|
||||
int pi3 = pindex[locfaces[0].PNum(3)];
|
||||
(*testout) << "pi = " << pi1 << ", " << pi2 << ", " << pi3 << endl;
|
||||
}
|
||||
|
||||
@ -369,7 +369,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
|
||||
grouppindex, groupfindex);
|
||||
|
||||
bool onlytri = 1;
|
||||
for (i = 0; i < groupfaces.Size(); i++)
|
||||
for (auto i : groupfaces.Range())
|
||||
if (groupfaces[i].GetNP() != 3)
|
||||
onlytri = 0;
|
||||
|
||||
@ -441,31 +441,28 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
|
||||
{
|
||||
// set transformatino to reference coordinates
|
||||
|
||||
if (locfaces.Get(1).GetNP() == 3)
|
||||
if (locfaces[0].GetNP() == 3)
|
||||
{
|
||||
trans.Set (locpoints.Get(locfaces.Get(1).PNumMod(1+rotind)),
|
||||
locpoints.Get(locfaces.Get(1).PNumMod(2+rotind)),
|
||||
locpoints.Get(locfaces.Get(1).PNumMod(3+rotind)), hshould);
|
||||
trans.Set (locpoints[locfaces[0].PNumMod(1+rotind)],
|
||||
locpoints[locfaces[0].PNumMod(2+rotind)],
|
||||
locpoints[locfaces[0].PNumMod(3+rotind)], hshould);
|
||||
}
|
||||
else
|
||||
{
|
||||
trans.Set (locpoints.Get(locfaces.Get(1).PNumMod(1+rotind)),
|
||||
locpoints.Get(locfaces.Get(1).PNumMod(2+rotind)),
|
||||
locpoints.Get(locfaces.Get(1).PNumMod(4+rotind)), hshould);
|
||||
trans.Set (locpoints[locfaces[0].PNumMod(1+rotind)],
|
||||
locpoints[locfaces[0].PNumMod(2+rotind)],
|
||||
locpoints[locfaces[0].PNumMod(4+rotind)], hshould);
|
||||
}
|
||||
|
||||
trans.ToPlain (locpoints, plainpoints);
|
||||
// trans.ToPlain (locpoints, plainpoints);
|
||||
|
||||
plainpoints.SetSize (locpoints.Size());
|
||||
for (auto i : locpoints.Range())
|
||||
trans.ToPlain (locpoints[i], plainpoints[i]);
|
||||
|
||||
for (i = 1; i <= allowpoint.Size(); i++)
|
||||
{
|
||||
if (plainpoints.Get(i).Z() > 0)
|
||||
{
|
||||
//if(loktestmode)
|
||||
// (*testout) << "plainpoints.Get(i).Z() = " << plainpoints.Get(i).Z() << " > 0" << endl;
|
||||
allowpoint.Elem(i) = 0;
|
||||
}
|
||||
}
|
||||
for (auto i : allowpoint.Range())
|
||||
if (plainpoints[i].Z() > 0)
|
||||
allowpoint[i] = false;
|
||||
|
||||
stat.cnttrials++;
|
||||
|
||||
@ -508,7 +505,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
|
||||
if (found) stat.cntsucc++;
|
||||
|
||||
locpoints.SetSize (plainpoints.Size());
|
||||
for (i = oldnp+1; i <= plainpoints.Size(); i++)
|
||||
for (int i = oldnp+1; i <= plainpoints.Size(); i++)
|
||||
trans.FromPlain (plainpoints.Elem(i), locpoints.Elem(i));
|
||||
|
||||
|
||||
@ -516,13 +513,13 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
|
||||
// avoid meshing from large to small mesh-size
|
||||
if (uselocalh && found && stat.qualclass <= 3)
|
||||
{
|
||||
for (i = 1; i <= locelements.Size(); i++)
|
||||
for (int i = 1; i <= locelements.Size(); i++)
|
||||
{
|
||||
Point3d pmin = locpoints.Get(locelements.Get(i).PNum(1));
|
||||
Point3d pmin = locpoints[locelements.Get(i).PNum(1)];
|
||||
Point3d pmax = pmin;
|
||||
for (j = 2; j <= 4; j++)
|
||||
{
|
||||
const Point3d & hp = locpoints.Get(locelements.Get(i).PNum(j));
|
||||
const Point3d & hp = locpoints[locelements.Get(i).PNum(j)];
|
||||
pmin.SetToMin (hp);
|
||||
pmax.SetToMax (hp);
|
||||
}
|
||||
@ -533,10 +530,10 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
|
||||
}
|
||||
if (found)
|
||||
{
|
||||
for (i = 1; i <= locelements.Size(); i++)
|
||||
for (j = 1; j <= 4; j++)
|
||||
for (int i = 1; i <= locelements.Size(); i++)
|
||||
for (int j = 1; j <= 4; j++)
|
||||
{
|
||||
const Point3d & hp = locpoints.Get(locelements.Get(i).PNum(j));
|
||||
const Point3d & hp = locpoints[locelements.Get(i).PNum(j)];
|
||||
if (Dist (hp, pmid) > hinner)
|
||||
found = 0;
|
||||
}
|
||||
@ -655,25 +652,27 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
|
||||
|
||||
pindex.SetSize(locpoints.Size());
|
||||
|
||||
for (i = oldnp+1; i <= locpoints.Size(); i++)
|
||||
for (int i = oldnp+1; i <= locpoints.Size(); i++)
|
||||
{
|
||||
PointIndex globind = mesh.AddPoint (locpoints.Get(i));
|
||||
pindex.Elem(i) = adfront -> AddPoint (locpoints.Get(i), globind);
|
||||
}
|
||||
|
||||
for (i = 1; i <= locelements.Size(); i++)
|
||||
cout << "bbb121212" << endl;
|
||||
|
||||
for (int i = 1; i <= locelements.Size(); i++)
|
||||
{
|
||||
Point3d * hp1, * hp2, * hp3, * hp4;
|
||||
hp1 = &locpoints.Elem(locelements.Get(i).PNum(1));
|
||||
hp2 = &locpoints.Elem(locelements.Get(i).PNum(2));
|
||||
hp3 = &locpoints.Elem(locelements.Get(i).PNum(3));
|
||||
hp4 = &locpoints.Elem(locelements.Get(i).PNum(4));
|
||||
hp1 = &locpoints[locelements.Get(i).PNum(1)];
|
||||
hp2 = &locpoints[locelements.Get(i).PNum(2)];
|
||||
hp3 = &locpoints[locelements.Get(i).PNum(3)];
|
||||
hp4 = &locpoints[locelements.Get(i).PNum(4)];
|
||||
|
||||
tetvol += (1.0 / 6.0) * ( Cross ( *hp2 - *hp1, *hp3 - *hp1) * (*hp4 - *hp1) );
|
||||
|
||||
for (j = 1; j <= locelements.Get(i).NP(); j++)
|
||||
locelements.Elem(i).PNum(j) =
|
||||
adfront -> GetGlobalIndex (pindex.Get(locelements.Get(i).PNum(j)));
|
||||
adfront -> GetGlobalIndex (pindex[locelements.Get(i).PNum(j)]);
|
||||
|
||||
mesh.AddVolumeElement (locelements.Get(i));
|
||||
stat.cntelem++;
|
||||
@ -683,7 +682,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
|
||||
{
|
||||
for (j = 1; j <= locfaces.Get(i).GetNP(); j++)
|
||||
locfaces.Elem(i).PNum(j) =
|
||||
pindex.Get(locfaces.Get(i).PNum(j));
|
||||
pindex[locfaces.Get(i).PNum(j)];
|
||||
// (*testout) << "add face " << locfaces.Get(i) << endl;
|
||||
adfront->AddFace (locfaces.Get(i));
|
||||
}
|
||||
|
||||
@ -42,7 +42,8 @@ public:
|
||||
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,
|
||||
INDEX_2_HASHTABLE<int> & connectedpairs,
|
||||
Array<Element> & elements,
|
||||
|
||||
@ -646,6 +646,8 @@ namespace netgen
|
||||
dshape(3,1) = (1-p(0));
|
||||
break;
|
||||
}
|
||||
default:
|
||||
throw NgException ("illegal element type in GetDShapeNew");
|
||||
}
|
||||
}
|
||||
|
||||
@ -1750,6 +1752,8 @@ namespace netgen
|
||||
{
|
||||
case TET: pp = &eltetqp[0][0]; break;
|
||||
case TET10: pp = &eltet10qp[ip-1][0]; break;
|
||||
default:
|
||||
throw NgException ("illegal element shape in GetIntegrationPoint");
|
||||
}
|
||||
|
||||
p(0) = pp[0];
|
||||
|
||||
@ -8,29 +8,29 @@ extern double minother;
|
||||
extern double minwithoutother;
|
||||
|
||||
|
||||
static double CalcElementBadness (const Array<Point3d> & points,
|
||||
const Element & elem)
|
||||
static double CalcElementBadness (const Array<Point3d, PointIndex::BASE> & points,
|
||||
const Element & elem)
|
||||
{
|
||||
double vol, l, l4, l5, l6;
|
||||
if (elem.GetNP() != 4)
|
||||
{
|
||||
if (elem.GetNP() == 5)
|
||||
{
|
||||
double z = points.Get(elem.PNum(5)).Z();
|
||||
double z = points[elem.PNum(5)].Z();
|
||||
if (z > -1e-8) return 1e8;
|
||||
return (-1 / z) - z; // - 2;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
Vec3d v1 = points.Get(elem.PNum(2)) - points.Get(elem.PNum(1));
|
||||
Vec3d v2 = points.Get(elem.PNum(3)) - points.Get(elem.PNum(1));
|
||||
Vec3d v3 = points.Get(elem.PNum(4)) - points.Get(elem.PNum(1));
|
||||
Vec3d v1 = points[elem.PNum(2)] - points[elem.PNum(1)];
|
||||
Vec3d v2 = points[elem.PNum(3)] - points[elem.PNum(1)];
|
||||
Vec3d v3 = points[elem.PNum(4)] - points[elem.PNum(1)];
|
||||
|
||||
vol = - (Cross (v1, v2) * v3);
|
||||
l4 = Dist (points.Get(elem.PNum(2)), points.Get(elem.PNum(3)));
|
||||
l5 = Dist (points.Get(elem.PNum(2)), points.Get(elem.PNum(4)));
|
||||
l6 = Dist (points.Get(elem.PNum(3)), points.Get(elem.PNum(4)));
|
||||
l4 = Dist (points[elem.PNum(2)], points[elem.PNum(3)]);
|
||||
l5 = Dist (points[elem.PNum(2)], points[elem.PNum(4)]);
|
||||
l6 = Dist (points[elem.PNum(3)], points[elem.PNum(4)]);
|
||||
|
||||
l = v1.Length() + v2.Length() + v3.Length() + l4 + l5 + l6;
|
||||
|
||||
@ -49,8 +49,8 @@ static double CalcElementBadness (const Array<Point3d> & points,
|
||||
|
||||
int Meshing3 :: ApplyRules
|
||||
(
|
||||
Array<Point3d> & 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<Point3d, PointIndex::BASE> & lpoints, // in: local points, out: old+new local points
|
||||
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
|
||||
INDEX lfacesplit, // for local faces in outer radius
|
||||
INDEX_2_HASHTABLE<int> & connectedpairs, // connected pairs for prism-meshing
|
||||
@ -65,25 +65,23 @@ int Meshing3 :: ApplyRules
|
||||
{
|
||||
NgProfiler::RegionTimer regtot(97);
|
||||
|
||||
int i, j, k, ri, nfok, npok, incnpok, refpi, locpi, locfi, locfr;
|
||||
float hf, err, minerr, teterr, minteterr;
|
||||
float err, minerr, teterr, minteterr;
|
||||
char ok, found, hc;
|
||||
vnetrule * rule;
|
||||
// vnetrule * rule;
|
||||
Vector oldu, newu, newu1, newu2, allp;
|
||||
Vec3d ui;
|
||||
Point3d np;
|
||||
int oldnp, noldlp, noldlf;
|
||||
const MiniElement2d * locface = NULL;
|
||||
int loktestmode;
|
||||
|
||||
|
||||
Array<int> pused; // point is already mapped
|
||||
Array<char> fused; // face is already mapped
|
||||
Array<int> pmap; // map of reference point to local point
|
||||
Array<char> pfixed; // point mapped by face-map
|
||||
Array<int> fmapi; // face in reference is mapped to face nr ...
|
||||
Array<int> fmapr; // face in reference is rotated to map
|
||||
Array<Point3d> transfreezone; // transformed free-zone
|
||||
Array<int, PointIndex::BASE> pused; // point is already mapped, number of uses
|
||||
Array<char> fused; // face is already mapped
|
||||
Array<PointIndex> pmap; // map of reference point to local point
|
||||
Array<bool, PointIndex::BASE> pfixed; // point mapped by face-map
|
||||
Array<int> fmapi; // face in reference is mapped to face nr ...
|
||||
Array<int> fmapr; // face in reference is rotated to map
|
||||
Array<Point3d> transfreezone; // transformed free-zone
|
||||
INDEX_2_CLOSED_HASHTABLE<int> ledges(100); // edges in local environment
|
||||
|
||||
Array<Point3d> tempnewpoints;
|
||||
@ -108,8 +106,9 @@ int Meshing3 :: ApplyRules
|
||||
fnearness.SetSize (lfacesplit);
|
||||
|
||||
pnearness = INT_MAX/10;
|
||||
for (j = 0; j < lfaces[0].GetNP(); j++)
|
||||
pnearness[lfaces[0][j]] = 0;
|
||||
|
||||
for (PointIndex pi : lfaces[0].PNums())
|
||||
pnearness[pi] = 0;
|
||||
|
||||
NgProfiler::RegionTimer reg2(98);
|
||||
|
||||
@ -118,24 +117,24 @@ int Meshing3 :: ApplyRules
|
||||
for (int loop = 0; loop < 2; loop++)
|
||||
{
|
||||
|
||||
for (i = 0; i < lfacesplit; i++)
|
||||
for (int i = 0; i < lfacesplit; i++)
|
||||
{
|
||||
const MiniElement2d & hface = lfaces[i];
|
||||
|
||||
int minn = INT_MAX-1;
|
||||
for (j = 0; j < hface.GetNP(); j++)
|
||||
for (PointIndex pi : hface.PNums())
|
||||
{
|
||||
int hi = pnearness[hface[j]];
|
||||
int hi = pnearness[pi];
|
||||
if (hi < minn) minn = hi;
|
||||
}
|
||||
if (minn < INT_MAX/10)
|
||||
for (j = 0; j < hface.GetNP(); j++)
|
||||
if (pnearness[hface[j]] > minn+1)
|
||||
pnearness[hface[j]] = minn+1;
|
||||
for (PointIndex pi : hface.PNums())
|
||||
if (pnearness[pi] > minn+1)
|
||||
pnearness[pi] = minn+1;
|
||||
}
|
||||
|
||||
for (i = 1; i <= connectedpairs.GetNBags(); i++)
|
||||
for (j = 1; j <= connectedpairs.GetBagSize(i); j++)
|
||||
for (int i = 1; i <= connectedpairs.GetNBags(); i++)
|
||||
for (int j = 1; j <= connectedpairs.GetBagSize(i); j++)
|
||||
{
|
||||
INDEX_2 edge;
|
||||
int val;
|
||||
@ -147,14 +146,13 @@ int Meshing3 :: ApplyRules
|
||||
if (pnearness[edge.I2()] > pnearness[edge.I1()] + 1)
|
||||
pnearness[edge.I2()] = pnearness[edge.I1()] + 1;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
for (i = 0; i < fnearness.Size(); i++)
|
||||
for (int i : fnearness.Range())
|
||||
{
|
||||
int sum = 0;
|
||||
for (j = 0; j < lfaces[i].GetNP(); j++)
|
||||
sum += pnearness[lfaces[i][j]];
|
||||
for (PointIndex pi : lfaces[i].PNums())
|
||||
sum += pnearness[pi];
|
||||
fnearness[i] = sum;
|
||||
}
|
||||
|
||||
@ -165,34 +163,35 @@ int Meshing3 :: ApplyRules
|
||||
// find bounding boxes of faces
|
||||
|
||||
triboxes.SetSize (lfaces.Size());
|
||||
for (i = 0; i < lfaces.Size(); i++)
|
||||
// for (int i = 0; i < lfaces.Size(); i++)
|
||||
for (auto i : lfaces.Range())
|
||||
{
|
||||
const MiniElement2d & face = lfaces[i];
|
||||
triboxes[i].SetPoint (lpoints.Get(face[0]));
|
||||
for (j = 1; j < face.GetNP(); j++)
|
||||
triboxes[i].AddPoint (lpoints.Get(face[j]));
|
||||
triboxes[i].SetPoint (lpoints[face[0]]);
|
||||
for (int j = 1; j < face.GetNP(); j++)
|
||||
triboxes[i].AddPoint (lpoints[face[j]]);
|
||||
}
|
||||
|
||||
NgProfiler::StopTimer (91);
|
||||
NgProfiler::StartTimer (92);
|
||||
|
||||
|
||||
bool useedges = 0;
|
||||
for (ri = 0; ri < rules.Size(); ri++)
|
||||
if (rules[ri]->GetNEd()) useedges = 1;
|
||||
bool useedges = false;
|
||||
for (int ri = 0; ri < rules.Size(); ri++)
|
||||
if (rules[ri]->GetNEd()) useedges = true;
|
||||
|
||||
if (useedges)
|
||||
{
|
||||
ledges.SetSize (5 * lfacesplit);
|
||||
|
||||
for (j = 0; j < lfacesplit; j++)
|
||||
for (int j = 0; j < lfacesplit; j++)
|
||||
// if (fnearness[j] <= 5)
|
||||
{
|
||||
const MiniElement2d & face = lfaces[j];
|
||||
int newp, oldp;
|
||||
|
||||
newp = face[face.GetNP()-1];
|
||||
for (k = 0; k < face.GetNP(); k++)
|
||||
for (int k = 0; k < face.GetNP(); k++)
|
||||
{
|
||||
oldp = newp;
|
||||
newp = face[k];
|
||||
@ -223,7 +222,7 @@ int Meshing3 :: ApplyRules
|
||||
|
||||
// check each rule:
|
||||
|
||||
for (ri = 1; ri <= rules.Size(); ri++)
|
||||
for (int ri = 1; ri <= rules.Size(); ri++)
|
||||
{
|
||||
int base = (lfaces[0].GetNP() == 3) ? 100 : 200;
|
||||
NgProfiler::RegionTimer regx1(base);
|
||||
@ -232,7 +231,7 @@ int Meshing3 :: ApplyRules
|
||||
// sprintf (problems.Elem(ri), "");
|
||||
*problems.Elem(ri) = '\0';
|
||||
|
||||
rule = rules.Get(ri);
|
||||
vnetrule * rule = rules.Get(ri);
|
||||
|
||||
if (rule->GetNP(1) != lfaces[0].GetNP())
|
||||
continue;
|
||||
@ -260,21 +259,21 @@ int Meshing3 :: ApplyRules
|
||||
|
||||
fused = 0;
|
||||
pused = 0;
|
||||
pmap = 0;
|
||||
for (auto & p : pmap) p.Invalidate();
|
||||
fmapi = 0;
|
||||
for (i = 1; i <= fmapr.Size(); i++)
|
||||
fmapr.Set(i, rule->GetNP(i));
|
||||
|
||||
for (int i : fmapr.Range())
|
||||
fmapr[i] = rule->GetNP(i+1);
|
||||
|
||||
fused[0] = 1;
|
||||
fmapi[0] = 1;
|
||||
fmapr[0] = rotind1;
|
||||
|
||||
|
||||
for (j = 1; j <= lfaces.Get(1).GetNP(); j++)
|
||||
for (int j = 1; j <= lfaces[0].GetNP(); j++)
|
||||
{
|
||||
locpi = lfaces[0].PNumMod (j+rotind1);
|
||||
PointIndex locpi = lfaces[0].PNumMod (j+rotind1);
|
||||
pmap.Set (rule->GetPointNr (1, j), locpi);
|
||||
pused.Elem(locpi)++;
|
||||
pused[locpi]++;
|
||||
}
|
||||
|
||||
/*
|
||||
@ -282,7 +281,7 @@ int Meshing3 :: ApplyRules
|
||||
nfok .. first nfok-1 faces are mapped properly
|
||||
*/
|
||||
|
||||
nfok = 2;
|
||||
int nfok = 2;
|
||||
NgProfiler::RegionTimer regfa(300);
|
||||
NgProfiler::RegionTimer regx2(base+50+ri);
|
||||
while (nfok >= 2)
|
||||
@ -293,8 +292,8 @@ int Meshing3 :: ApplyRules
|
||||
// not all faces mapped
|
||||
|
||||
ok = 0;
|
||||
locfi = fmapi.Get(nfok);
|
||||
locfr = fmapr.Get(nfok);
|
||||
int locfi = fmapi.Get(nfok);
|
||||
int locfr = fmapr.Get(nfok);
|
||||
|
||||
int actfnp = rule->GetNP(nfok);
|
||||
|
||||
@ -326,28 +325,27 @@ int Meshing3 :: ApplyRules
|
||||
|
||||
|
||||
// reference point already mapped differently ?
|
||||
for (j = 1; j <= actfnp && ok; j++)
|
||||
for (int j = 1; j <= actfnp && ok; j++)
|
||||
{
|
||||
locpi = pmap.Get(rule->GetPointNr (nfok, j));
|
||||
|
||||
if (locpi && locpi != locface->PNumMod(j+locfr))
|
||||
PointIndex locpi = pmap.Get(rule->GetPointNr (nfok, j));
|
||||
if (locpi.IsValid() && locpi != locface->PNumMod(j+locfr))
|
||||
ok = 0;
|
||||
}
|
||||
|
||||
// local point already used or point outside tolerance ?
|
||||
for (j = 1; j <= actfnp && ok; j++)
|
||||
for (int j = 1; j <= actfnp && ok; j++)
|
||||
{
|
||||
refpi = rule->GetPointNr (nfok, j);
|
||||
int refpi = rule->GetPointNr (nfok, j);
|
||||
|
||||
if (pmap.Get(refpi) == 0)
|
||||
if (!pmap.Get(refpi).IsValid())
|
||||
{
|
||||
locpi = locface->PNumMod (j + locfr);
|
||||
PointIndex locpi = locface->PNumMod (j + locfr);
|
||||
|
||||
if (pused.Get(locpi))
|
||||
if (pused[locpi])
|
||||
ok = 0;
|
||||
else
|
||||
{
|
||||
const Point3d & lp = lpoints.Get(locpi);
|
||||
const Point3d & lp = lpoints[locpi];
|
||||
const Point3d & rp = rule->GetPoint(refpi);
|
||||
|
||||
if ( Dist2 (lp, rp) * rule->PointDistFactor(refpi) > minerr)
|
||||
@ -370,16 +368,16 @@ int Meshing3 :: ApplyRules
|
||||
fmapr.Set (nfok, locfr);
|
||||
fused.Set (locfi, 1);
|
||||
|
||||
for (j = 1; j <= rule->GetNP (nfok); j++)
|
||||
for (int j = 1; j <= rule->GetNP (nfok); j++)
|
||||
{
|
||||
locpi = locface->PNumMod(j+locfr);
|
||||
PointIndex locpi = locface->PNumMod(j+locfr);
|
||||
|
||||
if (rule->GetPointNr (nfok, j) <= 3 &&
|
||||
pmap.Get(rule->GetPointNr(nfok, j)) != locpi)
|
||||
(*testout) << "change face1 point, mark1" << endl;
|
||||
|
||||
pmap.Set(rule->GetPointNr (nfok, j), locpi);
|
||||
pused.Elem(locpi)++;
|
||||
pused[locpi]++;
|
||||
}
|
||||
|
||||
nfok++;
|
||||
@ -392,14 +390,15 @@ int Meshing3 :: ApplyRules
|
||||
nfok--;
|
||||
|
||||
fused.Set (fmapi.Get(nfok), 0);
|
||||
for (j = 1; j <= rule->GetNP (nfok); j++)
|
||||
for (int j = 1; j <= rule->GetNP (nfok); j++)
|
||||
{
|
||||
refpi = rule->GetPointNr (nfok, j);
|
||||
pused.Elem(pmap.Get(refpi))--;
|
||||
int refpi = rule->GetPointNr (nfok, j);
|
||||
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();
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -419,12 +418,16 @@ int Meshing3 :: ApplyRules
|
||||
sprintf (problems.Elem(ri), "Faces Ok");
|
||||
}
|
||||
|
||||
npok = 1;
|
||||
incnpok = 1;
|
||||
int npok = 1;
|
||||
int incnpok = 1;
|
||||
|
||||
pfixed.SetSize (pmap.Size());
|
||||
for (i = 1; i <= pmap.Size(); i++)
|
||||
/*
|
||||
for (int i = 1; i <= pmap.Size(); i++)
|
||||
pfixed.Set(i, (pmap.Get(i) != 0) );
|
||||
*/
|
||||
for (int i : pmap.Range())
|
||||
pfixed[i] = pmap[i].IsValid();
|
||||
|
||||
while (npok >= 1)
|
||||
{
|
||||
@ -444,31 +447,31 @@ int Meshing3 :: ApplyRules
|
||||
else
|
||||
|
||||
{
|
||||
locpi = pmap.Elem(npok);
|
||||
PointIndex locpi = pmap.Elem(npok);
|
||||
ok = 0;
|
||||
|
||||
if (locpi)
|
||||
pused.Elem(locpi)--;
|
||||
if (locpi.IsValid())
|
||||
pused[locpi]--;
|
||||
|
||||
while (!ok && locpi < lpoints.Size())
|
||||
while (!ok && locpi < lpoints.Size()-1+PointIndex::BASE)
|
||||
{
|
||||
ok = 1;
|
||||
locpi++;
|
||||
|
||||
if (pused.Get(locpi) ||
|
||||
pnearness.Get(locpi) > rule->GetPNearness(npok))
|
||||
if (pused[locpi] ||
|
||||
pnearness[locpi] > rule->GetPNearness(npok))
|
||||
{
|
||||
ok = 0;
|
||||
}
|
||||
else if (allowpoint.Get(locpi) != 2)
|
||||
else if (allowpoint[locpi] != 2)
|
||||
{
|
||||
ok = 0;
|
||||
if (allowpoint.Get(locpi) == 1)
|
||||
if (allowpoint[locpi] == 1)
|
||||
impossible = 0;
|
||||
}
|
||||
else
|
||||
{
|
||||
const Point3d & lp = lpoints.Get(locpi);
|
||||
const Point3d & lp = lpoints[locpi];
|
||||
const Point3d & rp = rule->GetPoint(npok);
|
||||
|
||||
if ( Dist2 (lp, rp) * rule->PointDistFactor(npok) > minerr)
|
||||
@ -487,7 +490,7 @@ int Meshing3 :: ApplyRules
|
||||
if (npok <= 3)
|
||||
(*testout) << "set face1 point, mark3" << endl;
|
||||
|
||||
pused.Elem(locpi)++;
|
||||
pused[locpi]++;
|
||||
npok++;
|
||||
incnpok = 1;
|
||||
}
|
||||
@ -495,7 +498,8 @@ int Meshing3 :: ApplyRules
|
||||
else
|
||||
|
||||
{
|
||||
pmap.Set (npok, 0);
|
||||
// pmap.Set (npok, 0);
|
||||
pmap.Elem(npok).Invalidate();
|
||||
|
||||
if (npok <= 3)
|
||||
(*testout) << "set face1 point, mark4" << endl;
|
||||
@ -516,8 +520,8 @@ int Meshing3 :: ApplyRules
|
||||
if (loktestmode)
|
||||
{
|
||||
(*testout) << "Mapping found!!: Rule " << rule->Name() << endl;
|
||||
for (i = 1; i <= pmap.Size(); i++)
|
||||
(*testout) << pmap.Get(i) << " ";
|
||||
for (auto pi : pmap)
|
||||
(*testout) << pi << " ";
|
||||
(*testout) << endl;
|
||||
sprintf (problems.Elem(ri), "mapping found");
|
||||
(*testout) << rule->GetNP(1) << " = " << lfaces.Get(1).GetNP() << endl;
|
||||
@ -527,7 +531,7 @@ int Meshing3 :: ApplyRules
|
||||
|
||||
|
||||
// check mapedges:
|
||||
for (i = 1; i <= rule->GetNEd(); i++)
|
||||
for (int i = 1; i <= rule->GetNEd(); i++)
|
||||
{
|
||||
INDEX_2 in2(pmap.Get(rule->GetEdge(i).i1),
|
||||
pmap.Get(rule->GetEdge(i).i2));
|
||||
@ -537,12 +541,12 @@ int Meshing3 :: ApplyRules
|
||||
|
||||
|
||||
// check prism edges:
|
||||
for (i = 1; i <= rule->GetNE(); i++)
|
||||
for (int i = 1; i <= rule->GetNE(); i++)
|
||||
{
|
||||
const Element & el = rule->GetElement (i);
|
||||
if (el.GetType() == PRISM)
|
||||
{
|
||||
for (j = 1; j <= 3; j++)
|
||||
for (int j = 1; j <= 3; j++)
|
||||
{
|
||||
INDEX_2 in2(pmap.Get(el.PNum(j)),
|
||||
pmap.Get(el.PNum(j+3)));
|
||||
@ -554,7 +558,7 @@ int Meshing3 :: ApplyRules
|
||||
{
|
||||
if (loktestmode)
|
||||
(*testout) << "map pyramid, rule = " << rule->Name() << endl;
|
||||
for (j = 1; j <= 2; j++)
|
||||
for (int j = 1; j <= 2; j++)
|
||||
{
|
||||
INDEX_2 in2;
|
||||
if (j == 1)
|
||||
@ -581,7 +585,7 @@ int Meshing3 :: ApplyRules
|
||||
|
||||
|
||||
|
||||
for (i = rule->GetNOldF() + 1; i <= rule->GetNF(); i++)
|
||||
for (int i = rule->GetNOldF() + 1; i <= rule->GetNF(); i++)
|
||||
fmapi.Set(i, 0);
|
||||
|
||||
|
||||
@ -599,9 +603,9 @@ int Meshing3 :: ApplyRules
|
||||
newu.SetSize (3 * (rule->GetNP() - rule->GetNOldP()));
|
||||
allp.SetSize (3 * rule->GetNP());
|
||||
|
||||
for (i = 1; i <= rule->GetNOldP(); i++)
|
||||
for (int i = 1; i <= rule->GetNOldP(); i++)
|
||||
{
|
||||
const Point3d & lp = lpoints.Get(pmap.Get(i));
|
||||
const Point3d & lp = lpoints[pmap.Get(i)];
|
||||
const Point3d & rp = rule->GetPoint(i);
|
||||
oldu (3*i-3) = lp.X()-rp.X();
|
||||
oldu (3*i-2) = lp.Y()-rp.Y();
|
||||
@ -620,7 +624,7 @@ int Meshing3 :: ApplyRules
|
||||
|
||||
// int idiff = 3 * (rule->GetNP()-rule->GetNOldP());
|
||||
int idiff = 3 * rule->GetNOldP();
|
||||
for (i = rule->GetNOldP()+1; i <= rule->GetNP(); i++)
|
||||
for (int i = rule->GetNOldP()+1; i <= rule->GetNP(); i++)
|
||||
{
|
||||
const Point3d & rp = rule->GetPoint(i);
|
||||
allp (3*i-3) = rp.X() + newu(3*i-3 - idiff);
|
||||
@ -644,14 +648,14 @@ int Meshing3 :: ApplyRules
|
||||
{
|
||||
const Array<Point3d> & fz = rule->GetTransFreeZone();
|
||||
(*testout) << "Freezone: " << endl;
|
||||
for (i = 1; i <= fz.Size(); i++)
|
||||
for (int i = 1; i <= fz.Size(); i++)
|
||||
(*testout) << fz.Get(i) << endl;
|
||||
}
|
||||
|
||||
|
||||
// check freezone:
|
||||
|
||||
for (i = 1; i <= lpoints.Size(); i++)
|
||||
for (int i = 1; i <= lpoints.Size(); i++)
|
||||
{
|
||||
if ( !pused.Get(i) )
|
||||
{
|
||||
@ -675,7 +679,7 @@ int Meshing3 :: ApplyRules
|
||||
}
|
||||
}
|
||||
|
||||
for (i = 1; i <= lfaces.Size() && ok; i++)
|
||||
for (int i = 1; i <= lfaces.Size() && ok; i++)
|
||||
{
|
||||
static Array<int> lpi(4);
|
||||
|
||||
@ -692,7 +696,7 @@ int Meshing3 :: ApplyRules
|
||||
for (li = 1; li <= lfacei.GetNP(); li++)
|
||||
{
|
||||
int lpii = 0;
|
||||
int pi = lfacei.PNum(li);
|
||||
PointIndex pi = lfacei.PNum(li);
|
||||
for (lj = 1; lj <= rule->GetNOldP(); lj++)
|
||||
if (pmap.Get(lj) == pi)
|
||||
lpii = lj;
|
||||
@ -704,19 +708,19 @@ int Meshing3 :: ApplyRules
|
||||
{
|
||||
triin = rule->IsTriangleInFreeZone
|
||||
(
|
||||
lpoints.Get(lfacei.PNum(1)),
|
||||
lpoints.Get(lfacei.PNum(2)),
|
||||
lpoints.Get(lfacei.PNum(3)), lpi, 1
|
||||
lpoints[lfacei.PNum(1)],
|
||||
lpoints[lfacei.PNum(2)],
|
||||
lpoints[lfacei.PNum(3)], lpi, 1
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
triin = rule->IsQuadInFreeZone
|
||||
(
|
||||
lpoints.Get(lfacei.PNum(1)),
|
||||
lpoints.Get(lfacei.PNum(2)),
|
||||
lpoints.Get(lfacei.PNum(3)),
|
||||
lpoints.Get(lfacei.PNum(4)),
|
||||
lpoints[lfacei.PNum(1)],
|
||||
lpoints[lfacei.PNum(2)],
|
||||
lpoints[lfacei.PNum(3)],
|
||||
lpoints[lfacei.PNum(4)],
|
||||
lpi, 1
|
||||
);
|
||||
}
|
||||
@ -741,7 +745,7 @@ int Meshing3 :: ApplyRules
|
||||
<< lfaces.Get(i).PNum(3) << " - "
|
||||
<< lfaces.Get(i).PNum(4) << endl;
|
||||
for (int lj = 1; lj <= lfaces.Get(i).GetNP(); lj++)
|
||||
(*testout) << lpoints.Get(lfaces.Get(i).PNum(lj)) << " ";
|
||||
(*testout) << lpoints[lfaces.Get(i).PNum(lj)] << " ";
|
||||
|
||||
(*testout) << endl;
|
||||
|
||||
@ -759,9 +763,9 @@ int Meshing3 :: ApplyRules
|
||||
<< lfacei.PNum(2) << " - "
|
||||
<< lfacei.PNum(3)
|
||||
<< ", or "
|
||||
<< lpoints.Get(lfacei.PNum(1)) << " - "
|
||||
<< lpoints.Get(lfacei.PNum(2)) << " - "
|
||||
<< lpoints.Get(lfacei.PNum(3))
|
||||
<< lpoints[lfacei.PNum(1)] << " - "
|
||||
<< lpoints[lfacei.PNum(2)] << " - "
|
||||
<< lpoints[lfacei.PNum(3)]
|
||||
<< endl;
|
||||
(*testout) << "lpi = " << lpi.Get(1) << ", "
|
||||
<< lpi.Get(2) << ", " << lpi.Get(3) << endl;
|
||||
@ -773,10 +777,10 @@ int Meshing3 :: ApplyRules
|
||||
<< lfacei.PNum(3) << " - "
|
||||
<< lfacei.PNum(4)
|
||||
<< ", or "
|
||||
<< lpoints.Get(lfacei.PNum(1)) << " - "
|
||||
<< lpoints.Get(lfacei.PNum(2)) << " - "
|
||||
<< lpoints.Get(lfacei.PNum(3)) << " - "
|
||||
<< lpoints.Get(lfacei.PNum(4))
|
||||
<< lpoints[lfacei.PNum(1)] << " - "
|
||||
<< lpoints[lfacei.PNum(2)] << " - "
|
||||
<< lpoints[lfacei.PNum(3)] << " - "
|
||||
<< lpoints[lfacei.PNum(4)]
|
||||
<< endl;
|
||||
|
||||
sprintf (problems.Elem(ri), "triangle (%d, %d, %d) in Freezone",
|
||||
@ -786,13 +790,13 @@ int Meshing3 :: ApplyRules
|
||||
}
|
||||
|
||||
hc = 0;
|
||||
for (k = rule->GetNOldF() + 1; k <= rule->GetNF(); k++)
|
||||
for (int k = rule->GetNOldF() + 1; k <= rule->GetNF(); k++)
|
||||
{
|
||||
if (rule->GetPointNr(k, 1) <= rule->GetNOldP() &&
|
||||
rule->GetPointNr(k, 2) <= rule->GetNOldP() &&
|
||||
rule->GetPointNr(k, 3) <= rule->GetNOldP())
|
||||
{
|
||||
for (j = 1; j <= 3; j++)
|
||||
for (int j = 1; j <= 3; j++)
|
||||
if (lfaces.Get(i).PNumMod(j ) == pmap.Get(rule->GetPointNr(k, 1)) &&
|
||||
lfaces.Get(i).PNumMod(j+1) == pmap.Get(rule->GetPointNr(k, 3)) &&
|
||||
lfaces.Get(i).PNumMod(j+2) == pmap.Get(rule->GetPointNr(k, 2)))
|
||||
@ -839,9 +843,9 @@ int Meshing3 :: ApplyRules
|
||||
if (ok)
|
||||
{
|
||||
err = 0;
|
||||
for (i = 1; i <= rule->GetNOldP(); i++)
|
||||
for (int i = 1; i <= rule->GetNOldP(); i++)
|
||||
{
|
||||
hf = rule->CalcPointDist (i, lpoints.Get(pmap.Get(i)));
|
||||
double hf = rule->CalcPointDist (i, lpoints[pmap.Get(i)]);
|
||||
if (hf > err) err = hf;
|
||||
}
|
||||
|
||||
@ -856,40 +860,37 @@ int Meshing3 :: ApplyRules
|
||||
// newu = rule->GetOldUToNewU() * oldu;
|
||||
|
||||
// set new points:
|
||||
int oldnp = rule->GetNOldP();
|
||||
int noldlp = lpoints.Size();
|
||||
int noldlf = lfaces.Size();
|
||||
|
||||
oldnp = rule->GetNOldP();
|
||||
noldlp = lpoints.Size();
|
||||
noldlf = lfaces.Size();
|
||||
|
||||
|
||||
for (i = oldnp + 1; i <= rule->GetNP(); i++)
|
||||
for (int i = oldnp + 1; i <= rule->GetNP(); i++)
|
||||
{
|
||||
np = rule->GetPoint(i);
|
||||
np.X() += newu (3 * (i-oldnp) - 3);
|
||||
np.Y() += newu (3 * (i-oldnp) - 2);
|
||||
np.Z() += newu (3 * (i-oldnp) - 1);
|
||||
|
||||
pmap.Elem(i) = lpoints.Append (np);
|
||||
pmap.Elem(i) = lpoints.Append (np)-1+PointIndex::BASE;
|
||||
}
|
||||
|
||||
// Set new Faces:
|
||||
|
||||
for (i = rule->GetNOldF() + 1; i <= rule->GetNF(); i++)
|
||||
for (int i = rule->GetNOldF() + 1; i <= rule->GetNF(); i++)
|
||||
if (!fmapi.Get(i))
|
||||
{
|
||||
MiniElement2d nface(rule->GetNP(i));
|
||||
for (j = 1; j <= nface.GetNP(); j++)
|
||||
for (int j = 1; j <= nface.GetNP(); j++)
|
||||
nface.PNum(j) = pmap.Get(rule->GetPointNr (i, j));
|
||||
|
||||
lfaces.Append (nface);
|
||||
}
|
||||
|
||||
|
||||
// Delete old Faces:
|
||||
|
||||
for (i = 1; i <= rule->GetNDelF(); i++)
|
||||
for (int i = 1; i <= rule->GetNDelF(); i++)
|
||||
delfaces.Append (fmapi.Get(rule->GetDelFace(i)));
|
||||
for (i = rule->GetNOldF()+1; i <= rule->GetNF(); i++)
|
||||
for (int i = rule->GetNOldF()+1; i <= rule->GetNF(); i++)
|
||||
if (fmapi.Get(i))
|
||||
{
|
||||
delfaces.Append (fmapi.Get(i));
|
||||
@ -898,17 +899,17 @@ int Meshing3 :: ApplyRules
|
||||
|
||||
|
||||
// check orientation
|
||||
for (i = 1; i <= rule->GetNO() && ok; i++)
|
||||
for (int i = 1; i <= rule->GetNO() && ok; i++)
|
||||
{
|
||||
const fourint * fouri;
|
||||
|
||||
fouri = &rule->GetOrientation(i);
|
||||
Vec3d v1 (lpoints.Get(pmap.Get(fouri->i1)),
|
||||
lpoints.Get(pmap.Get(fouri->i2)));
|
||||
Vec3d v2 (lpoints.Get(pmap.Get(fouri->i1)),
|
||||
lpoints.Get(pmap.Get(fouri->i3)));
|
||||
Vec3d v3 (lpoints.Get(pmap.Get(fouri->i1)),
|
||||
lpoints.Get(pmap.Get(fouri->i4)));
|
||||
Vec3d v1 (lpoints[pmap.Get(fouri->i1)],
|
||||
lpoints[pmap.Get(fouri->i2)]);
|
||||
Vec3d v2 (lpoints[pmap.Get(fouri->i1)],
|
||||
lpoints[pmap.Get(fouri->i3)]);
|
||||
Vec3d v3 (lpoints[pmap.Get(fouri->i1)],
|
||||
lpoints[pmap.Get(fouri->i4)]);
|
||||
|
||||
Vec3d n;
|
||||
Cross (v1, v2, n);
|
||||
@ -927,7 +928,7 @@ int Meshing3 :: ApplyRules
|
||||
|
||||
|
||||
// new points in free-zone ?
|
||||
for (i = rule->GetNOldP() + 1; i <= rule->GetNP() && ok; i++)
|
||||
for (int i = rule->GetNOldP() + 1; i <= rule->GetNP() && ok; i++)
|
||||
if (!rule->IsInFreeZone (lpoints.Get(pmap.Get(i))))
|
||||
{
|
||||
if (loktestmode)
|
||||
@ -942,10 +943,10 @@ int Meshing3 :: ApplyRules
|
||||
|
||||
// insert new elements
|
||||
|
||||
for (i = 1; i <= rule->GetNE(); i++)
|
||||
for (int i = 1; i <= rule->GetNE(); i++)
|
||||
{
|
||||
elements.Append (rule->GetElement(i));
|
||||
for (j = 1; j <= elements.Get(i).NP(); j++)
|
||||
for (int j = 1; j <= elements.Get(i).NP(); j++)
|
||||
elements.Elem(i).PNum(j) = pmap.Get(elements.Get(i).PNum(j));
|
||||
}
|
||||
|
||||
@ -953,9 +954,9 @@ int Meshing3 :: ApplyRules
|
||||
// Calculate Element badness
|
||||
|
||||
teterr = 0;
|
||||
for (i = 1; i <= elements.Size(); i++)
|
||||
for (int i = 1; i <= elements.Size(); i++)
|
||||
{
|
||||
hf = CalcElementBadness (lpoints, elements.Get(i));
|
||||
double hf = CalcElementBadness (lpoints, elements.Get(i));
|
||||
if (hf > teterr) teterr = hf;
|
||||
}
|
||||
|
||||
@ -978,29 +979,29 @@ int Meshing3 :: ApplyRules
|
||||
double oldlen = 0;
|
||||
double newlen = 0;
|
||||
|
||||
for (i = 1; i <= rule->GetNDelF(); i++)
|
||||
for (int i = 1; i <= rule->GetNDelF(); i++)
|
||||
{
|
||||
const Element2d & face =
|
||||
rule->GetFace (rule->GetDelFace(i));
|
||||
for (j = 1; j <= 3; j++)
|
||||
for (int j = 1; j <= 3; j++)
|
||||
{
|
||||
const Point3d & p1 =
|
||||
lpoints.Get(pmap.Get(face.PNumMod(j)));
|
||||
lpoints[pmap.Get(face.PNumMod(j))];
|
||||
const Point3d & p2 =
|
||||
lpoints.Get(pmap.Get(face.PNumMod(j+1)));
|
||||
lpoints[pmap.Get(face.PNumMod(j+1))];
|
||||
oldlen += Dist(p1, p2);
|
||||
}
|
||||
}
|
||||
|
||||
for (i = rule->GetNOldF()+1; i <= rule->GetNF(); i++)
|
||||
for (int i = rule->GetNOldF()+1; i <= rule->GetNF(); i++)
|
||||
{
|
||||
const Element2d & face = rule->GetFace (i);
|
||||
for (j = 1; j <= 3; j++)
|
||||
for (int j = 1; j <= 3; j++)
|
||||
{
|
||||
const Point3d & p1 =
|
||||
lpoints.Get(pmap.Get(face.PNumMod(j)));
|
||||
lpoints[pmap.Get(face.PNumMod(j))];
|
||||
const Point3d & p2 =
|
||||
lpoints.Get(pmap.Get(face.PNumMod(j+1)));
|
||||
lpoints[pmap.Get(face.PNumMod(j+1))];
|
||||
newlen += Dist(p1, p2);
|
||||
}
|
||||
}
|
||||
@ -1057,7 +1058,7 @@ int Meshing3 :: ApplyRules
|
||||
|
||||
if (testmode)
|
||||
{
|
||||
for (i = 1; i <= rule->GetNOldP(); i++)
|
||||
for (int i = 1; i <= rule->GetNOldP(); i++)
|
||||
{
|
||||
(*testout) << "P" << i << ": Ref: "
|
||||
<< rule->GetPoint (i) << " is: "
|
||||
@ -1066,19 +1067,19 @@ int Meshing3 :: ApplyRules
|
||||
}
|
||||
|
||||
tempnewpoints.SetSize (0);
|
||||
for (i = noldlp+1; i <= lpoints.Size(); i++)
|
||||
for (int i = noldlp+1; i <= lpoints.Size(); i++)
|
||||
tempnewpoints.Append (lpoints.Get(i));
|
||||
|
||||
tempnewfaces.SetSize (0);
|
||||
for (i = noldlf+1; i <= lfaces.Size(); i++)
|
||||
for (int i = noldlf+1; i <= lfaces.Size(); i++)
|
||||
tempnewfaces.Append (lfaces.Get(i));
|
||||
|
||||
tempdelfaces.SetSize (0);
|
||||
for (i = 1; i <= delfaces.Size(); i++)
|
||||
for (int i = 1; i <= delfaces.Size(); i++)
|
||||
tempdelfaces.Append (delfaces.Get(i));
|
||||
|
||||
tempelements.SetSize (0);
|
||||
for (i = 1; i <= elements.Size(); i++)
|
||||
for (int i = 1; i <= elements.Size(); i++)
|
||||
tempelements.Append (elements.Get(i));
|
||||
}
|
||||
|
||||
@ -1096,15 +1097,13 @@ int Meshing3 :: ApplyRules
|
||||
|
||||
nfok = rule->GetNOldF();
|
||||
|
||||
for (j = 1; j <= rule->GetNP (nfok); j++)
|
||||
for (int j = 1; j <= rule->GetNP (nfok); j++)
|
||||
{
|
||||
refpi = rule->GetPointNr (nfok, j);
|
||||
pused.Elem(pmap.Get(refpi))--;
|
||||
int refpi = rule->GetPointNr (nfok, j);
|
||||
pused[pmap.Get(refpi)]--;
|
||||
|
||||
if (pused.Get(pmap.Get(refpi)) == 0)
|
||||
{
|
||||
pmap.Set(refpi, 0);
|
||||
}
|
||||
if (pused[pmap.Get(refpi)] == 0)
|
||||
pmap.Elem(refpi).Invalidate();
|
||||
}
|
||||
|
||||
}
|
||||
@ -1115,15 +1114,32 @@ int Meshing3 :: ApplyRules
|
||||
|
||||
if (found)
|
||||
{
|
||||
/*
|
||||
for (i = 1; i <= tempnewpoints.Size(); i++)
|
||||
lpoints.Append (tempnewpoints.Get(i));
|
||||
*/
|
||||
for (Point3d p : tempnewpoints)
|
||||
lpoints.Append(p);
|
||||
/*
|
||||
for (i = 1; i <= tempnewfaces.Size(); i++)
|
||||
if (tempnewfaces.Get(i).PNum(1))
|
||||
lfaces.Append (tempnewfaces.Get(i));
|
||||
*/
|
||||
for (int i : tempnewfaces.Range())
|
||||
if (tempnewfaces[i].PNum(1).IsValid())
|
||||
lfaces.Append (tempnewfaces[i]);
|
||||
/*
|
||||
for (i = 1; i <= tempdelfaces.Size(); i++)
|
||||
delfaces.Append (tempdelfaces.Get(i));
|
||||
*/
|
||||
for (int i : tempdelfaces.Range())
|
||||
delfaces.Append (tempdelfaces[i]);
|
||||
/*
|
||||
for (i = 1; i <= tempelements.Size(); i++)
|
||||
elements.Append (tempelements.Get(i));
|
||||
*/
|
||||
for (int i : tempelements.Range())
|
||||
elements.Append (tempelements[i]);
|
||||
}
|
||||
|
||||
retminerr = minerr;
|
||||
|
||||
Loading…
Reference in New Issue
Block a user