mirror of
https://github.com/NGSolve/netgen.git
synced 2024-12-24 21:10:33 +05:00
Mesh 3d domains in parallel
To get consistent results, copy the LocalH tree in BlockFillLocalH
This commit is contained in:
Matthias Hochsteger
parent
5e3505b897
commit
97623db219
@ -91,6 +91,45 @@ namespace netgen
|
||||
delete root;
|
||||
}
|
||||
|
||||
unique_ptr<LocalH> LocalH :: Copy ()
|
||||
{
|
||||
static Timer t("LocalH::Copy"); RegionTimer rt(t);
|
||||
auto lh = make_unique<LocalH>(boundingbox, grading, dimension);
|
||||
std::map<GradingBox*, GradingBox*> mapping;
|
||||
lh->boxes.SetSize(boxes.Size());
|
||||
|
||||
for(auto i : boxes.Range())
|
||||
{
|
||||
lh->boxes[i] = new GradingBox();
|
||||
auto & bnew = *lh->boxes[i];
|
||||
auto & b = *boxes[i];
|
||||
bnew.xmid[0] = b.xmid[0];
|
||||
bnew.xmid[1] = b.xmid[1];
|
||||
bnew.xmid[2] = b.xmid[2];
|
||||
bnew.h2 = b.h2;
|
||||
bnew.hopt = b.hopt;
|
||||
bnew.flags = b.flags;
|
||||
mapping[&b] = &bnew;
|
||||
}
|
||||
|
||||
for(auto i : boxes.Range())
|
||||
{
|
||||
auto & bnew = *lh->boxes[i];
|
||||
auto & b = *boxes[i];
|
||||
for(auto k : Range(8))
|
||||
{
|
||||
if(b.childs[k])
|
||||
bnew.childs[k] = mapping[b.childs[k]];
|
||||
}
|
||||
|
||||
if(b.father)
|
||||
bnew.father = mapping[b.father];
|
||||
}
|
||||
|
||||
lh->root = mapping[root];
|
||||
return lh;
|
||||
}
|
||||
|
||||
void LocalH :: Delete ()
|
||||
{
|
||||
root->DeleteChilds();
|
||||
|
||||
@ -98,6 +98,8 @@ namespace netgen
|
||||
|
||||
~LocalH();
|
||||
///
|
||||
unique_ptr<LocalH> Copy();
|
||||
///
|
||||
void Delete();
|
||||
///
|
||||
void DoArchive(Archive& ar);
|
||||
|
||||
@ -10,18 +10,299 @@ namespace netgen
|
||||
extern const char * pyramidrules2[];
|
||||
extern const char * hexrules[];
|
||||
|
||||
void MeshDomain(Mesh & mesh3d, const MeshingParameters & c_mp, int k)
|
||||
{
|
||||
MeshingParameters mp = c_mp; // copy mp to change them here
|
||||
NgArray<INDEX_2> connectednodes;
|
||||
|
||||
int oldne;
|
||||
int meshed;
|
||||
if(mp.only3D_domain_nr && mp.only3D_domain_nr !=k)
|
||||
return;
|
||||
if (multithread.terminate)
|
||||
return;
|
||||
|
||||
PrintMessage (2, "");
|
||||
PrintMessage (1, "Meshing subdomain ", k, " of ", mesh3d.GetNDomains());
|
||||
(*testout) << "Meshing subdomain " << k << endl;
|
||||
|
||||
mp.maxh = min2 (mp.maxh, mesh3d.MaxHDomain(k));
|
||||
|
||||
mesh3d.CalcSurfacesOfNode();
|
||||
mesh3d.FindOpenElements(k);
|
||||
|
||||
if (!mesh3d.GetNOpenElements())
|
||||
return;
|
||||
|
||||
|
||||
|
||||
Box<3> domain_bbox( Box<3>::EMPTY_BOX );
|
||||
|
||||
for (SurfaceElementIndex sei = 0; sei < mesh3d.GetNSE(); sei++)
|
||||
{
|
||||
const Element2d & el = mesh3d[sei];
|
||||
if (el.IsDeleted() ) continue;
|
||||
|
||||
if (mesh3d.GetFaceDescriptor(el.GetIndex()).DomainIn() == k ||
|
||||
mesh3d.GetFaceDescriptor(el.GetIndex()).DomainOut() == k)
|
||||
|
||||
for (int j = 0; j < el.GetNP(); j++)
|
||||
domain_bbox.Add (mesh3d[el[j]]);
|
||||
}
|
||||
domain_bbox.Increase (0.01 * domain_bbox.Diam());
|
||||
|
||||
|
||||
for (int qstep = 0; qstep <= 3; qstep++)
|
||||
// for (int qstep = 0; qstep <= 0; qstep++) // for hex-filling
|
||||
{
|
||||
if (qstep == 0 && !mp.try_hexes) continue;
|
||||
|
||||
// cout << "openquads = " << mesh3d.HasOpenQuads() << endl;
|
||||
if (mesh3d.HasOpenQuads())
|
||||
{
|
||||
string rulefile = ngdir;
|
||||
|
||||
const char ** rulep = NULL;
|
||||
switch (qstep)
|
||||
{
|
||||
case 0:
|
||||
rulefile = "/Users/joachim/gitlab/netgen/rules/hexa.rls";
|
||||
rulep = hexrules;
|
||||
break;
|
||||
case 1:
|
||||
rulefile += "/rules/prisms2.rls";
|
||||
rulep = prismrules2;
|
||||
break;
|
||||
case 2: // connect pyramid to triangle
|
||||
rulefile += "/rules/pyramids2.rls";
|
||||
rulep = pyramidrules2;
|
||||
break;
|
||||
case 3: // connect to vis-a-vis point
|
||||
rulefile += "/rules/pyramids.rls";
|
||||
rulep = pyramidrules;
|
||||
break;
|
||||
}
|
||||
|
||||
// Meshing3 meshing(rulefile);
|
||||
Meshing3 meshing(rulep);
|
||||
|
||||
MeshingParameters mpquad = mp;
|
||||
|
||||
mpquad.giveuptol = 15;
|
||||
mpquad.baseelnp = 4;
|
||||
mpquad.starshapeclass = 1000;
|
||||
mpquad.check_impossible = qstep == 1; // for prisms only (air domain in trafo)
|
||||
|
||||
|
||||
// for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
|
||||
for (PointIndex pi : mesh3d.Points().Range())
|
||||
meshing.AddPoint (mesh3d[pi], pi);
|
||||
|
||||
/*
|
||||
mesh3d.GetIdentifications().GetPairs (0, connectednodes);
|
||||
for (int i = 1; i <= connectednodes.Size(); i++)
|
||||
meshing.AddConnectedPair (connectednodes.Get(i));
|
||||
*/
|
||||
for (int nr = 1; nr <= mesh3d.GetIdentifications().GetMaxNr(); nr++)
|
||||
if (mesh3d.GetIdentifications().GetType(nr) != Identifications::PERIODIC)
|
||||
{
|
||||
mesh3d.GetIdentifications().GetPairs (nr, connectednodes);
|
||||
for (auto pair : connectednodes)
|
||||
meshing.AddConnectedPair (pair);
|
||||
}
|
||||
|
||||
for (int i = 1; i <= mesh3d.GetNOpenElements(); i++)
|
||||
{
|
||||
Element2d hel = mesh3d.OpenElement(i);
|
||||
meshing.AddBoundaryElement (hel);
|
||||
}
|
||||
|
||||
oldne = mesh3d.GetNE();
|
||||
|
||||
meshing.GenerateMesh (mesh3d, mpquad);
|
||||
|
||||
for (int i = oldne + 1; i <= mesh3d.GetNE(); i++)
|
||||
mesh3d.VolumeElement(i).SetIndex (k);
|
||||
|
||||
(*testout)
|
||||
<< "mesh has " << mesh3d.GetNE() << " prism/pyramid elements" << endl;
|
||||
|
||||
mesh3d.FindOpenElements(k);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if (mesh3d.HasOpenQuads())
|
||||
{
|
||||
PrintSysError ("mesh has still open quads");
|
||||
throw NgException ("Stop meshing since too many attempts");
|
||||
// return MESHING3_GIVEUP;
|
||||
}
|
||||
|
||||
|
||||
if (mp.delaunay && mesh3d.GetNOpenElements())
|
||||
{
|
||||
Meshing3 meshing((const char**)NULL);
|
||||
|
||||
mesh3d.FindOpenElements(k);
|
||||
|
||||
/*
|
||||
for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
|
||||
meshing.AddPoint (mesh3d[pi], pi);
|
||||
*/
|
||||
for (PointIndex pi : mesh3d.Points().Range())
|
||||
meshing.AddPoint (mesh3d[pi], pi);
|
||||
|
||||
for (int i = 1; i <= mesh3d.GetNOpenElements(); i++)
|
||||
meshing.AddBoundaryElement (mesh3d.OpenElement(i));
|
||||
|
||||
oldne = mesh3d.GetNE();
|
||||
|
||||
meshing.Delaunay (mesh3d, k, mp);
|
||||
|
||||
for (int i = oldne + 1; i <= mesh3d.GetNE(); i++)
|
||||
mesh3d.VolumeElement(i).SetIndex (k);
|
||||
|
||||
PrintMessage (3, mesh3d.GetNP(), " points, ",
|
||||
mesh3d.GetNE(), " elements");
|
||||
}
|
||||
|
||||
|
||||
int cntsteps = 0;
|
||||
if (mesh3d.GetNOpenElements())
|
||||
do
|
||||
{
|
||||
if (multithread.terminate)
|
||||
break;
|
||||
|
||||
mesh3d.FindOpenElements(k);
|
||||
PrintMessage (5, mesh3d.GetNOpenElements(), " open faces");
|
||||
cntsteps++;
|
||||
|
||||
if (cntsteps > mp.maxoutersteps)
|
||||
throw NgException ("Stop meshing since too many attempts");
|
||||
|
||||
string rulefile = ngdir + "/tetra.rls";
|
||||
PrintMessage (1, "start tetmeshing");
|
||||
|
||||
// Meshing3 meshing(rulefile);
|
||||
Meshing3 meshing(tetrules);
|
||||
|
||||
NgArray<int, PointIndex::BASE> glob2loc(mesh3d.GetNP());
|
||||
glob2loc = -1;
|
||||
|
||||
// for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
|
||||
for (PointIndex pi : mesh3d.Points().Range())
|
||||
if (domain_bbox.IsIn (mesh3d[pi]))
|
||||
glob2loc[pi] =
|
||||
meshing.AddPoint (mesh3d[pi], pi);
|
||||
|
||||
for (int i = 1; i <= mesh3d.GetNOpenElements(); i++)
|
||||
{
|
||||
Element2d hel = mesh3d.OpenElement(i);
|
||||
for (int j = 0; j < hel.GetNP(); j++)
|
||||
hel[j] = glob2loc[hel[j]];
|
||||
meshing.AddBoundaryElement (hel);
|
||||
// meshing.AddBoundaryElement (mesh3d.OpenElement(i));
|
||||
}
|
||||
|
||||
oldne = mesh3d.GetNE();
|
||||
|
||||
mp.giveuptol = 15 + 10 * cntsteps;
|
||||
mp.sloppy = 5;
|
||||
meshing.GenerateMesh (mesh3d, mp);
|
||||
|
||||
for (ElementIndex ei = oldne; ei < mesh3d.GetNE(); ei++)
|
||||
mesh3d[ei].SetIndex (k);
|
||||
|
||||
|
||||
mesh3d.CalcSurfacesOfNode();
|
||||
mesh3d.FindOpenElements(k);
|
||||
|
||||
// teterrpow = 2;
|
||||
if (mesh3d.GetNOpenElements() != 0)
|
||||
{
|
||||
meshed = 0;
|
||||
PrintMessage (5, mesh3d.GetNOpenElements(), " open faces found");
|
||||
|
||||
MeshOptimize3d optmesh(mp);
|
||||
|
||||
const char * optstr = "mcmstmcmstmcmstmcm";
|
||||
for (size_t j = 1; j <= strlen(optstr); j++)
|
||||
{
|
||||
mesh3d.CalcSurfacesOfNode();
|
||||
mesh3d.FreeOpenElementsEnvironment(2);
|
||||
mesh3d.CalcSurfacesOfNode();
|
||||
|
||||
switch (optstr[j-1])
|
||||
{
|
||||
case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break;
|
||||
case 'd': optmesh.SplitImprove(mesh3d, OPT_REST); break;
|
||||
case 's': optmesh.SwapImprove(mesh3d, OPT_REST); break;
|
||||
case 't': optmesh.SwapImprove2(mesh3d, OPT_REST); break;
|
||||
case 'm': mesh3d.ImproveMesh(mp, OPT_REST); break;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
mesh3d.FindOpenElements(k);
|
||||
PrintMessage (3, "Call remove problem");
|
||||
RemoveProblem (mesh3d, k);
|
||||
mesh3d.FindOpenElements(k);
|
||||
}
|
||||
else
|
||||
{
|
||||
meshed = 1;
|
||||
PrintMessage (1, "Success !");
|
||||
}
|
||||
}
|
||||
while (!meshed);
|
||||
|
||||
PrintMessage (1, mesh3d.GetNP(), " points, ",
|
||||
mesh3d.GetNE(), " elements");
|
||||
{
|
||||
if(mp.only3D_domain_nr && mp.only3D_domain_nr !=k)
|
||||
return;
|
||||
PrintMessage (3, "Check subdomain ", k, " / ", mesh3d.GetNDomains());
|
||||
|
||||
mesh3d.FindOpenElements(k);
|
||||
|
||||
bool res = (mesh3d.CheckConsistentBoundary() != 0);
|
||||
if (res)
|
||||
{
|
||||
PrintError ("Surface mesh not consistent");
|
||||
throw NgException ("Stop meshing since surface mesh not consistent");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void MergeMeshes( Mesh & mesh, FlatArray<Mesh> meshes, PointIndex first_new_pi )
|
||||
{
|
||||
static Timer t("MergeMeshes"); RegionTimer rt(t);
|
||||
for(auto & m : meshes)
|
||||
{
|
||||
Array<PointIndex, PointIndex> pmap(m.Points().Size());
|
||||
for(auto pi : Range(PointIndex(PointIndex::BASE), first_new_pi))
|
||||
pmap[pi] = pi;
|
||||
|
||||
for (auto pi : Range(first_new_pi, m.Points().Range().Next()))
|
||||
pmap[pi] = mesh.AddPoint(m[pi]);
|
||||
|
||||
|
||||
for ( auto el : m.VolumeElements() )
|
||||
{
|
||||
for (auto i : Range(el.GetNP()))
|
||||
el[i] = pmap[el[i]];
|
||||
mesh.AddVolumeElement(el);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// extern double teterrpow;
|
||||
MESHING3_RESULT MeshVolume (const MeshingParameters & c_mp, Mesh& mesh3d)
|
||||
MESHING3_RESULT MeshVolume (const MeshingParameters & mp, Mesh& mesh3d)
|
||||
{
|
||||
static Timer t("MeshVolume"); RegionTimer reg(t);
|
||||
|
||||
MeshingParameters mp = c_mp; // copy mp to change them here
|
||||
int oldne;
|
||||
int meshed;
|
||||
|
||||
NgArray<INDEX_2> connectednodes;
|
||||
|
||||
if (!mesh3d.HasLocalHFunction()) mesh3d.CalcLocalH(mp.grading);
|
||||
|
||||
mesh3d.Compress();
|
||||
@ -32,288 +313,32 @@ namespace netgen
|
||||
if (mesh3d.CheckOverlappingBoundary())
|
||||
throw NgException ("Stop meshing since boundary mesh is overlapping");
|
||||
|
||||
int nonconsist = 0;
|
||||
for (int k = 1; k <= mesh3d.GetNDomains(); k++)
|
||||
|
||||
if(task_manager)
|
||||
{
|
||||
if(mp.only3D_domain_nr && mp.only3D_domain_nr !=k)
|
||||
continue;
|
||||
PrintMessage (3, "Check subdomain ", k, " / ", mesh3d.GetNDomains());
|
||||
Array<Mesh> meshes(mesh3d.GetNDomains()-1);
|
||||
auto first_new_pi = mesh3d.Points().Range().Next();
|
||||
|
||||
mesh3d.FindOpenElements(k);
|
||||
|
||||
/*
|
||||
bool res = mesh3d.CheckOverlappingBoundary();
|
||||
if (res)
|
||||
{
|
||||
PrintError ("Surface is overlapping !!");
|
||||
nonconsist = 1;
|
||||
}
|
||||
*/
|
||||
|
||||
bool res = (mesh3d.CheckConsistentBoundary() != 0);
|
||||
if (res)
|
||||
{
|
||||
PrintError ("Surface mesh not consistent");
|
||||
nonconsist = 1;
|
||||
}
|
||||
}
|
||||
|
||||
if (nonconsist)
|
||||
{
|
||||
PrintError ("Stop meshing since surface mesh not consistent");
|
||||
throw NgException ("Stop meshing since surface mesh not consistent");
|
||||
}
|
||||
|
||||
double globmaxh = mp.maxh;
|
||||
|
||||
for (int k = 1; k <= mesh3d.GetNDomains(); k++)
|
||||
for(auto & m : meshes)
|
||||
{
|
||||
if(mp.only3D_domain_nr && mp.only3D_domain_nr !=k)
|
||||
continue;
|
||||
if (multithread.terminate)
|
||||
break;
|
||||
|
||||
PrintMessage (2, "");
|
||||
PrintMessage (1, "Meshing subdomain ", k, " of ", mesh3d.GetNDomains());
|
||||
(*testout) << "Meshing subdomain " << k << endl;
|
||||
|
||||
mp.maxh = min2 (globmaxh, mesh3d.MaxHDomain(k));
|
||||
|
||||
mesh3d.CalcSurfacesOfNode();
|
||||
mesh3d.FindOpenElements(k);
|
||||
|
||||
if (!mesh3d.GetNOpenElements())
|
||||
continue;
|
||||
|
||||
|
||||
|
||||
Box<3> domain_bbox( Box<3>::EMPTY_BOX );
|
||||
|
||||
for (SurfaceElementIndex sei = 0; sei < mesh3d.GetNSE(); sei++)
|
||||
{
|
||||
const Element2d & el = mesh3d[sei];
|
||||
if (el.IsDeleted() ) continue;
|
||||
|
||||
if (mesh3d.GetFaceDescriptor(el.GetIndex()).DomainIn() == k ||
|
||||
mesh3d.GetFaceDescriptor(el.GetIndex()).DomainOut() == k)
|
||||
|
||||
for (int j = 0; j < el.GetNP(); j++)
|
||||
domain_bbox.Add (mesh3d[el[j]]);
|
||||
}
|
||||
domain_bbox.Increase (0.01 * domain_bbox.Diam());
|
||||
|
||||
|
||||
for (int qstep = 0; qstep <= 3; qstep++)
|
||||
// for (int qstep = 0; qstep <= 0; qstep++) // for hex-filling
|
||||
{
|
||||
if (qstep == 0 && !mp.try_hexes) continue;
|
||||
|
||||
// cout << "openquads = " << mesh3d.HasOpenQuads() << endl;
|
||||
if (mesh3d.HasOpenQuads())
|
||||
{
|
||||
string rulefile = ngdir;
|
||||
|
||||
const char ** rulep = NULL;
|
||||
switch (qstep)
|
||||
{
|
||||
case 0:
|
||||
rulefile = "/Users/joachim/gitlab/netgen/rules/hexa.rls";
|
||||
rulep = hexrules;
|
||||
break;
|
||||
case 1:
|
||||
rulefile += "/rules/prisms2.rls";
|
||||
rulep = prismrules2;
|
||||
break;
|
||||
case 2: // connect pyramid to triangle
|
||||
rulefile += "/rules/pyramids2.rls";
|
||||
rulep = pyramidrules2;
|
||||
break;
|
||||
case 3: // connect to vis-a-vis point
|
||||
rulefile += "/rules/pyramids.rls";
|
||||
rulep = pyramidrules;
|
||||
break;
|
||||
}
|
||||
|
||||
// Meshing3 meshing(rulefile);
|
||||
Meshing3 meshing(rulep);
|
||||
|
||||
MeshingParameters mpquad = mp;
|
||||
|
||||
mpquad.giveuptol = 15;
|
||||
mpquad.baseelnp = 4;
|
||||
mpquad.starshapeclass = 1000;
|
||||
mpquad.check_impossible = qstep == 1; // for prisms only (air domain in trafo)
|
||||
|
||||
|
||||
// for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
|
||||
for (PointIndex pi : mesh3d.Points().Range())
|
||||
meshing.AddPoint (mesh3d[pi], pi);
|
||||
|
||||
/*
|
||||
mesh3d.GetIdentifications().GetPairs (0, connectednodes);
|
||||
for (int i = 1; i <= connectednodes.Size(); i++)
|
||||
meshing.AddConnectedPair (connectednodes.Get(i));
|
||||
*/
|
||||
for (int nr = 1; nr <= mesh3d.GetIdentifications().GetMaxNr(); nr++)
|
||||
if (mesh3d.GetIdentifications().GetType(nr) != Identifications::PERIODIC)
|
||||
{
|
||||
mesh3d.GetIdentifications().GetPairs (nr, connectednodes);
|
||||
for (auto pair : connectednodes)
|
||||
meshing.AddConnectedPair (pair);
|
||||
}
|
||||
|
||||
for (int i = 1; i <= mesh3d.GetNOpenElements(); i++)
|
||||
{
|
||||
Element2d hel = mesh3d.OpenElement(i);
|
||||
meshing.AddBoundaryElement (hel);
|
||||
}
|
||||
|
||||
oldne = mesh3d.GetNE();
|
||||
|
||||
meshing.GenerateMesh (mesh3d, mpquad);
|
||||
|
||||
for (int i = oldne + 1; i <= mesh3d.GetNE(); i++)
|
||||
mesh3d.VolumeElement(i).SetIndex (k);
|
||||
|
||||
(*testout)
|
||||
<< "mesh has " << mesh3d.GetNE() << " prism/pyramid elements" << endl;
|
||||
|
||||
mesh3d.FindOpenElements(k);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if (mesh3d.HasOpenQuads())
|
||||
{
|
||||
PrintSysError ("mesh has still open quads");
|
||||
throw NgException ("Stop meshing since too many attempts");
|
||||
// return MESHING3_GIVEUP;
|
||||
}
|
||||
|
||||
|
||||
if (mp.delaunay && mesh3d.GetNOpenElements())
|
||||
{
|
||||
Meshing3 meshing((const char**)NULL);
|
||||
|
||||
mesh3d.FindOpenElements(k);
|
||||
|
||||
/*
|
||||
for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
|
||||
meshing.AddPoint (mesh3d[pi], pi);
|
||||
*/
|
||||
for (PointIndex pi : mesh3d.Points().Range())
|
||||
meshing.AddPoint (mesh3d[pi], pi);
|
||||
|
||||
for (int i = 1; i <= mesh3d.GetNOpenElements(); i++)
|
||||
meshing.AddBoundaryElement (mesh3d.OpenElement(i));
|
||||
|
||||
oldne = mesh3d.GetNE();
|
||||
|
||||
meshing.Delaunay (mesh3d, k, mp);
|
||||
|
||||
for (int i = oldne + 1; i <= mesh3d.GetNE(); i++)
|
||||
mesh3d.VolumeElement(i).SetIndex (k);
|
||||
|
||||
PrintMessage (3, mesh3d.GetNP(), " points, ",
|
||||
mesh3d.GetNE(), " elements");
|
||||
}
|
||||
|
||||
|
||||
int cntsteps = 0;
|
||||
if (mesh3d.GetNOpenElements())
|
||||
do
|
||||
{
|
||||
if (multithread.terminate)
|
||||
break;
|
||||
|
||||
mesh3d.FindOpenElements(k);
|
||||
PrintMessage (5, mesh3d.GetNOpenElements(), " open faces");
|
||||
cntsteps++;
|
||||
|
||||
if (cntsteps > mp.maxoutersteps)
|
||||
throw NgException ("Stop meshing since too many attempts");
|
||||
|
||||
string rulefile = ngdir + "/tetra.rls";
|
||||
PrintMessage (1, "start tetmeshing");
|
||||
|
||||
// Meshing3 meshing(rulefile);
|
||||
Meshing3 meshing(tetrules);
|
||||
|
||||
NgArray<int, PointIndex::BASE> glob2loc(mesh3d.GetNP());
|
||||
glob2loc = -1;
|
||||
|
||||
// for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
|
||||
for (PointIndex pi : mesh3d.Points().Range())
|
||||
if (domain_bbox.IsIn (mesh3d[pi]))
|
||||
glob2loc[pi] =
|
||||
meshing.AddPoint (mesh3d[pi], pi);
|
||||
|
||||
for (int i = 1; i <= mesh3d.GetNOpenElements(); i++)
|
||||
{
|
||||
Element2d hel = mesh3d.OpenElement(i);
|
||||
for (int j = 0; j < hel.GetNP(); j++)
|
||||
hel[j] = glob2loc[hel[j]];
|
||||
meshing.AddBoundaryElement (hel);
|
||||
// meshing.AddBoundaryElement (mesh3d.OpenElement(i));
|
||||
}
|
||||
|
||||
oldne = mesh3d.GetNE();
|
||||
|
||||
mp.giveuptol = 15 + 10 * cntsteps;
|
||||
mp.sloppy = 5;
|
||||
meshing.GenerateMesh (mesh3d, mp);
|
||||
|
||||
for (ElementIndex ei = oldne; ei < mesh3d.GetNE(); ei++)
|
||||
mesh3d[ei].SetIndex (k);
|
||||
|
||||
|
||||
mesh3d.CalcSurfacesOfNode();
|
||||
mesh3d.FindOpenElements(k);
|
||||
|
||||
// teterrpow = 2;
|
||||
if (mesh3d.GetNOpenElements() != 0)
|
||||
{
|
||||
meshed = 0;
|
||||
PrintMessage (5, mesh3d.GetNOpenElements(), " open faces found");
|
||||
|
||||
MeshOptimize3d optmesh(mp);
|
||||
|
||||
const char * optstr = "mcmstmcmstmcmstmcm";
|
||||
for (size_t j = 1; j <= strlen(optstr); j++)
|
||||
{
|
||||
mesh3d.CalcSurfacesOfNode();
|
||||
mesh3d.FreeOpenElementsEnvironment(2);
|
||||
mesh3d.CalcSurfacesOfNode();
|
||||
|
||||
switch (optstr[j-1])
|
||||
{
|
||||
case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break;
|
||||
case 'd': optmesh.SplitImprove(mesh3d, OPT_REST); break;
|
||||
case 's': optmesh.SwapImprove(mesh3d, OPT_REST); break;
|
||||
case 't': optmesh.SwapImprove2(mesh3d, OPT_REST); break;
|
||||
case 'm': mesh3d.ImproveMesh(mp, OPT_REST); break;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
mesh3d.FindOpenElements(k);
|
||||
PrintMessage (3, "Call remove problem");
|
||||
RemoveProblem (mesh3d, k);
|
||||
mesh3d.FindOpenElements(k);
|
||||
}
|
||||
else
|
||||
{
|
||||
meshed = 1;
|
||||
PrintMessage (1, "Success !");
|
||||
}
|
||||
}
|
||||
while (!meshed);
|
||||
|
||||
PrintMessage (1, mesh3d.GetNP(), " points, ",
|
||||
mesh3d.GetNE(), " elements");
|
||||
m = mesh3d;
|
||||
m.SetLocalH(mesh3d.GetLocalH());
|
||||
}
|
||||
|
||||
mp.maxh = globmaxh;
|
||||
ParallelFor(Range(1, mesh3d.GetNDomains()+1), [&](int k)
|
||||
{
|
||||
if(k==1)
|
||||
MeshDomain(mesh3d, mp, k);
|
||||
else
|
||||
MeshDomain(meshes[k-2], mp, k);
|
||||
});
|
||||
MergeMeshes(mesh3d, meshes, first_new_pi);
|
||||
}
|
||||
else
|
||||
for (int k = 1; k <= mesh3d.GetNDomains(); k++)
|
||||
MeshDomain(mesh3d, mp, k);
|
||||
|
||||
|
||||
|
||||
MeshQuality3d (mesh3d);
|
||||
|
||||
|
||||
@ -1144,10 +1144,13 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
|
||||
|
||||
if (mp.maxh < maxh) maxh = mp.maxh;
|
||||
|
||||
auto loch_ptr = mesh.LocalHFunction().Copy();
|
||||
auto & loch = *loch_ptr;
|
||||
|
||||
bool changed;
|
||||
do
|
||||
{
|
||||
mesh.LocalHFunction().ClearFlags();
|
||||
loch.ClearFlags();
|
||||
|
||||
for (int i = 1; i <= adfront->GetNF(); i++)
|
||||
{
|
||||
@ -1161,21 +1164,21 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
|
||||
double filld = filldist * bbox.Diam();
|
||||
bbox.Increase (filld);
|
||||
|
||||
mesh.LocalHFunction().CutBoundary (bbox); // .PMin(), bbox.PMax());
|
||||
loch.CutBoundary (bbox); // .PMin(), bbox.PMax());
|
||||
}
|
||||
|
||||
// locadfront = adfront;
|
||||
mesh.LocalHFunction().FindInnerBoxes (adfront, NULL);
|
||||
loch.FindInnerBoxes (adfront, NULL);
|
||||
|
||||
npoints.SetSize(0);
|
||||
mesh.LocalHFunction().GetInnerPoints (npoints);
|
||||
loch.GetInnerPoints (npoints);
|
||||
|
||||
changed = false;
|
||||
for (int i = 1; i <= npoints.Size(); i++)
|
||||
{
|
||||
if (mesh.LocalHFunction().GetH(npoints.Get(i)) > 1.5 * maxh)
|
||||
if (loch.GetH(npoints.Get(i)) > 1.5 * maxh)
|
||||
{
|
||||
mesh.LocalHFunction().SetH (npoints.Get(i), maxh);
|
||||
loch.SetH (npoints.Get(i), maxh);
|
||||
changed = true;
|
||||
}
|
||||
}
|
||||
|
||||
@ -82,13 +82,13 @@
|
||||
],
|
||||
"angles_trig": [
|
||||
24.858,
|
||||
104.73
|
||||
107.08
|
||||
],
|
||||
"ne1d": 181,
|
||||
"ne2d": 313,
|
||||
"ne3d": 506,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 6, 13, 25, 47, 64, 85, 94, 100, 51, 17]",
|
||||
"total_badness": 650.35553401
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 4, 13, 28, 45, 66, 88, 88, 95, 57, 18]",
|
||||
"total_badness": 650.55380342
|
||||
}
|
||||
],
|
||||
"boxcyl.geo": [
|
||||
@ -1840,9 +1840,9 @@
|
||||
],
|
||||
"ne1d": 174,
|
||||
"ne2d": 1190,
|
||||
"ne3d": 5113,
|
||||
"quality_histogram": "[0, 0, 15, 128, 183, 43, 57, 114, 124, 184, 294, 387, 492, 590, 635, 626, 528, 418, 234, 61]",
|
||||
"total_badness": 8982.334633
|
||||
"ne3d": 5073,
|
||||
"quality_histogram": "[0, 0, 15, 128, 183, 44, 57, 112, 124, 171, 322, 373, 483, 622, 623, 558, 526, 429, 240, 63]",
|
||||
"total_badness": 8929.8790628
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
@ -1865,14 +1865,14 @@
|
||||
169.26
|
||||
],
|
||||
"angles_trig": [
|
||||
8.6573,
|
||||
161.61
|
||||
8.6562,
|
||||
161.62
|
||||
],
|
||||
"ne1d": 132,
|
||||
"ne2d": 812,
|
||||
"ne3d": 2514,
|
||||
"quality_histogram": "[0, 0, 7, 63, 69, 154, 134, 117, 148, 199, 247, 286, 257, 214, 177, 162, 120, 92, 45, 23]",
|
||||
"total_badness": 5239.4475113
|
||||
"ne3d": 2498,
|
||||
"quality_histogram": "[0, 0, 7, 63, 69, 153, 138, 115, 153, 201, 253, 273, 267, 201, 164, 161, 120, 95, 44, 21]",
|
||||
"total_badness": 5225.1671437
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
@ -1885,39 +1885,39 @@
|
||||
],
|
||||
"ne1d": 174,
|
||||
"ne2d": 1190,
|
||||
"ne3d": 5067,
|
||||
"quality_histogram": "[0, 0, 14, 117, 182, 41, 52, 114, 110, 155, 268, 313, 495, 590, 645, 645, 535, 485, 238, 68]",
|
||||
"total_badness": 8750.9467413
|
||||
"ne3d": 4995,
|
||||
"quality_histogram": "[0, 0, 14, 117, 182, 40, 52, 109, 106, 145, 275, 310, 458, 591, 622, 604, 541, 495, 267, 67]",
|
||||
"total_badness": 8621.0579835
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
12.76,
|
||||
147.23
|
||||
14.607,
|
||||
144.88
|
||||
],
|
||||
"angles_trig": [
|
||||
15.824,
|
||||
16.508,
|
||||
143.02
|
||||
],
|
||||
"ne1d": 248,
|
||||
"ne2d": 2320,
|
||||
"ne3d": 16409,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 4, 20, 55, 119, 192, 288, 630, 958, 1489, 2106, 2552, 2926, 2613, 1913, 544]",
|
||||
"total_badness": 21635.781365
|
||||
"ne3d": 16414,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 4, 22, 61, 107, 178, 334, 607, 1066, 1532, 2200, 2498, 2897, 2655, 1737, 516]",
|
||||
"total_badness": 21736.358053
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
18.203,
|
||||
17.436,
|
||||
144.68
|
||||
],
|
||||
"angles_trig": [
|
||||
17.821,
|
||||
130.51
|
||||
127.08
|
||||
],
|
||||
"ne1d": 418,
|
||||
"ne2d": 5958,
|
||||
"ne3d": 102414,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 1, 4, 6, 45, 124, 368, 1028, 2576, 5782, 10622, 16119, 21223, 22353, 16809, 5354]",
|
||||
"total_badness": 125921.99427
|
||||
"ne3d": 100894,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 1, 5, 7, 34, 114, 369, 997, 2522, 5452, 10512, 15927, 20701, 21949, 16958, 5346]",
|
||||
"total_badness": 123900.96306
|
||||
}
|
||||
],
|
||||
"ortho.geo": [
|
||||
@ -2092,8 +2092,8 @@
|
||||
"period.geo": [
|
||||
{
|
||||
"angles_tet": [
|
||||
13.348,
|
||||
152.73
|
||||
14.261,
|
||||
145.23
|
||||
],
|
||||
"angles_trig": [
|
||||
15.314,
|
||||
@ -2101,9 +2101,9 @@
|
||||
],
|
||||
"ne1d": 344,
|
||||
"ne2d": 1118,
|
||||
"ne3d": 3303,
|
||||
"quality_histogram": "[0, 0, 0, 0, 1, 3, 19, 28, 60, 79, 167, 274, 340, 411, 474, 453, 427, 342, 172, 53]",
|
||||
"total_badness": 4787.3615998
|
||||
"ne3d": 3285,
|
||||
"quality_histogram": "[0, 0, 0, 0, 1, 2, 16, 26, 56, 87, 157, 264, 331, 405, 486, 454, 436, 340, 173, 51]",
|
||||
"total_badness": 4744.2301558
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
@ -2133,12 +2133,12 @@
|
||||
"ne2d": 566,
|
||||
"ne3d": 1302,
|
||||
"quality_histogram": "[0, 0, 0, 1, 4, 17, 29, 39, 64, 86, 116, 123, 148, 143, 127, 134, 119, 83, 58, 11]",
|
||||
"total_badness": 2151.3920824
|
||||
"total_badness": 2151.3920827
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
15.291,
|
||||
144.92
|
||||
15.577,
|
||||
142.87
|
||||
],
|
||||
"angles_trig": [
|
||||
15.314,
|
||||
@ -2146,39 +2146,39 @@
|
||||
],
|
||||
"ne1d": 344,
|
||||
"ne2d": 1118,
|
||||
"ne3d": 3261,
|
||||
"quality_histogram": "[0, 0, 0, 0, 1, 2, 14, 23, 49, 62, 156, 225, 323, 419, 476, 458, 421, 385, 185, 62]",
|
||||
"total_badness": 4648.2224319
|
||||
"ne3d": 3240,
|
||||
"quality_histogram": "[0, 0, 0, 0, 1, 2, 15, 22, 54, 65, 154, 226, 315, 404, 494, 459, 406, 383, 179, 61]",
|
||||
"total_badness": 4628.23351
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
21.808,
|
||||
142.31
|
||||
21.896,
|
||||
143.75
|
||||
],
|
||||
"angles_trig": [
|
||||
23.022,
|
||||
128.64
|
||||
23.103,
|
||||
123.18
|
||||
],
|
||||
"ne1d": 480,
|
||||
"ne2d": 2248,
|
||||
"ne3d": 11618,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 8, 29, 90, 215, 539, 874, 1483, 1964, 2225, 2201, 1553, 437]",
|
||||
"total_badness": 14695.782197
|
||||
"ne3d": 11666,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 1, 4, 25, 100, 251, 514, 914, 1518, 1994, 2234, 2110, 1565, 436]",
|
||||
"total_badness": 14784.15449
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
22.578,
|
||||
141.63
|
||||
21.617,
|
||||
142.69
|
||||
],
|
||||
"angles_trig": [
|
||||
22.146,
|
||||
120.55
|
||||
122.89
|
||||
],
|
||||
"ne1d": 820,
|
||||
"ne2d": 6206,
|
||||
"ne3d": 68494,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 8, 53, 159, 556, 1496, 3478, 6612, 10914, 14587, 15239, 11731, 3661]",
|
||||
"total_badness": 83666.140475
|
||||
"ne3d": 68535,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 1, 5, 42, 157, 567, 1545, 3554, 6751, 10948, 14379, 15271, 11665, 3650]",
|
||||
"total_badness": 83784.962076
|
||||
}
|
||||
],
|
||||
"plane.stl": [
|
||||
@ -2760,24 +2760,24 @@
|
||||
],
|
||||
"ne1d": 74,
|
||||
"ne2d": 412,
|
||||
"ne3d": 1681,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 5, 5, 17, 18, 30, 48, 94, 133, 194, 249, 242, 240, 204, 159, 43]",
|
||||
"total_badness": 2334.8383469
|
||||
"ne3d": 1688,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 5, 5, 17, 18, 31, 48, 98, 133, 190, 250, 259, 245, 198, 142, 49]",
|
||||
"total_badness": 2349.0686154
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
25.029,
|
||||
138.94
|
||||
24.909,
|
||||
140.9
|
||||
],
|
||||
"angles_trig": [
|
||||
22.069,
|
||||
127.5
|
||||
21.973,
|
||||
127.7
|
||||
],
|
||||
"ne1d": 122,
|
||||
"ne2d": 1076,
|
||||
"ne3d": 14037,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 3, 24, 84, 179, 418, 822, 1431, 2256, 2852, 2929, 2328, 711]",
|
||||
"total_badness": 17344.477334
|
||||
"ne3d": 14075,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 3, 28, 68, 192, 457, 844, 1466, 2272, 2839, 2856, 2353, 697]",
|
||||
"total_badness": 17420.850186
|
||||
}
|
||||
],
|
||||
"square.in2d": [
|
||||
@ -3156,13 +3156,13 @@
|
||||
],
|
||||
"angles_trig": [
|
||||
14.916,
|
||||
130.79
|
||||
132.02
|
||||
],
|
||||
"ne1d": 690,
|
||||
"ne2d": 1670,
|
||||
"ne3d": 5168,
|
||||
"quality_histogram": "[0, 0, 1, 0, 0, 8, 30, 37, 106, 198, 284, 368, 447, 559, 691, 708, 596, 539, 460, 136]",
|
||||
"total_badness": 7457.2380943
|
||||
"ne3d": 5169,
|
||||
"quality_histogram": "[0, 0, 1, 0, 1, 7, 31, 35, 107, 196, 275, 376, 451, 556, 687, 707, 599, 545, 457, 138]",
|
||||
"total_badness": 7455.986288
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
@ -3191,8 +3191,8 @@
|
||||
"ne1d": 512,
|
||||
"ne2d": 866,
|
||||
"ne3d": 2373,
|
||||
"quality_histogram": "[0, 0, 0, 5, 9, 17, 44, 69, 124, 143, 189, 210, 313, 383, 339, 233, 138, 87, 46, 24]",
|
||||
"total_badness": 3943.0636847
|
||||
"quality_histogram": "[0, 0, 0, 5, 9, 17, 44, 69, 124, 143, 189, 211, 313, 381, 339, 234, 138, 87, 46, 24]",
|
||||
"total_badness": 3943.062114
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
@ -3201,13 +3201,13 @@
|
||||
],
|
||||
"angles_trig": [
|
||||
14.916,
|
||||
130.79
|
||||
132.02
|
||||
],
|
||||
"ne1d": 690,
|
||||
"ne2d": 1670,
|
||||
"ne3d": 5111,
|
||||
"quality_histogram": "[0, 0, 1, 0, 0, 3, 21, 38, 107, 192, 269, 348, 432, 542, 664, 707, 611, 564, 471, 141]",
|
||||
"total_badness": 7317.6330225
|
||||
"ne3d": 5117,
|
||||
"quality_histogram": "[0, 0, 1, 0, 0, 3, 22, 36, 106, 193, 265, 351, 428, 546, 668, 707, 609, 567, 474, 141]",
|
||||
"total_badness": 7321.564248
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
@ -3220,24 +3220,24 @@
|
||||
],
|
||||
"ne1d": 1050,
|
||||
"ne2d": 3784,
|
||||
"ne3d": 17749,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 3, 14, 35, 69, 180, 557, 1402, 2124, 2373, 2627, 2696, 2716, 2280, 673]",
|
||||
"total_badness": 23160.553922
|
||||
"ne3d": 17750,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 3, 14, 35, 68, 188, 555, 1411, 2124, 2368, 2623, 2700, 2707, 2279, 675]",
|
||||
"total_badness": 23168.68937
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
15.34,
|
||||
149.42
|
||||
14.338,
|
||||
149.44
|
||||
],
|
||||
"angles_trig": [
|
||||
20.032,
|
||||
19.234,
|
||||
129.78
|
||||
],
|
||||
"ne1d": 1722,
|
||||
"ne2d": 10022,
|
||||
"ne3d": 85138,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 3, 49, 1423, 714, 399, 655, 1222, 2404, 5670, 9020, 13424, 16402, 16956, 12744, 4053]",
|
||||
"total_badness": 108990.18852
|
||||
"ne3d": 84843,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 2, 47, 1403, 696, 390, 665, 1234, 2433, 5392, 8824, 13169, 16695, 17000, 12791, 4102]",
|
||||
"total_badness": 108464.06506
|
||||
}
|
||||
],
|
||||
"twobricks.geo": [
|
||||
@ -3318,18 +3318,18 @@
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
28.752,
|
||||
130.07
|
||||
28.509,
|
||||
131.79
|
||||
],
|
||||
"angles_trig": [
|
||||
25.5,
|
||||
109.19
|
||||
27.418,
|
||||
108.8
|
||||
],
|
||||
"ne1d": 186,
|
||||
"ne2d": 334,
|
||||
"ne3d": 596,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 21, 35, 49, 93, 105, 118, 99, 54, 17]",
|
||||
"total_badness": 770.34262233
|
||||
"ne3d": 589,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 20, 31, 52, 90, 109, 114, 101, 59, 8]",
|
||||
"total_badness": 760.70063244
|
||||
}
|
||||
],
|
||||
"twocubes.geo": [
|
||||
@ -3410,18 +3410,18 @@
|
||||
},
|
||||
{
|
||||
"angles_tet": [
|
||||
28.752,
|
||||
130.07
|
||||
28.509,
|
||||
131.79
|
||||
],
|
||||
"angles_trig": [
|
||||
25.5,
|
||||
109.19
|
||||
27.418,
|
||||
108.8
|
||||
],
|
||||
"ne1d": 186,
|
||||
"ne2d": 334,
|
||||
"ne3d": 596,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 21, 35, 49, 93, 105, 118, 99, 54, 17]",
|
||||
"total_badness": 770.34262233
|
||||
"ne3d": 589,
|
||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 20, 31, 52, 90, 109, 114, 101, 59, 8]",
|
||||
"total_badness": 760.70063244
|
||||
}
|
||||
],
|
||||
"twocyl.geo": [
|
||||
|
||||
Loading…
Reference in New Issue
Block a user