netgen/libsrc/meshing/meshfunc.cpp

745 lines
18 KiB
C++
Raw Normal View History

2009-01-13 04:40:13 +05:00
#include <mystdlib.h>
#include "meshing.hpp"
namespace netgen
{
extern const char * tetrules[];
// extern const char * tetrules2[];
extern const char * prismrules2[];
extern const char * pyramidrules[];
extern const char * pyramidrules2[];
extern const char * hexrules[];
2009-01-13 04:40:13 +05:00
2011-01-11 01:18:01 +05:00
// extern double teterrpow;
MESHING3_RESULT MeshVolume (const MeshingParameters & c_mp, Mesh& mesh3d)
2009-01-13 04:40:13 +05:00
{
2019-01-31 22:41:20 +05:00
static Timer t("MeshVolume"); RegionTimer reg(t);
MeshingParameters mp = c_mp; // copy mp to change them here
int oldne;
int meshed;
2009-01-13 04:40:13 +05:00
2019-07-09 13:39:16 +05:00
NgArray<INDEX_2> connectednodes;
2016-02-23 22:22:00 +05:00
if (!mesh3d.HasLocalHFunction()) mesh3d.CalcLocalH(mp.grading);
2009-02-25 21:06:34 +05:00
mesh3d.Compress();
2009-01-13 04:40:13 +05:00
// mesh3d.PrintMemInfo (cout);
2009-01-13 04:40:13 +05:00
if (mp.checkoverlappingboundary)
if (mesh3d.CheckOverlappingBoundary())
throw NgException ("Stop meshing since boundary mesh is overlapping");
2009-01-13 04:40:13 +05:00
int nonconsist = 0;
for (int k = 1; k <= mesh3d.GetNDomains(); k++)
{
if(mp.only3D_domain_nr && mp.only3D_domain_nr !=k)
continue;
PrintMessage (3, "Check subdomain ", k, " / ", mesh3d.GetNDomains());
2009-01-13 04:40:13 +05:00
mesh3d.FindOpenElements(k);
2009-01-13 04:40:13 +05:00
/*
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");
}
2009-01-13 04:40:13 +05:00
double globmaxh = mp.maxh;
2009-01-13 04:40:13 +05:00
for (int k = 1; k <= mesh3d.GetNDomains(); k++)
2010-02-26 22:06:00 +05:00
{
if(mp.only3D_domain_nr && mp.only3D_domain_nr !=k)
continue;
2010-02-26 22:06:00 +05:00
if (multithread.terminate)
break;
2010-02-26 22:06:00 +05:00
PrintMessage (2, "");
PrintMessage (1, "Meshing subdomain ", k, " of ", mesh3d.GetNDomains());
(*testout) << "Meshing subdomain " << k << endl;
mp.maxh = min2 (globmaxh, mesh3d.MaxHDomain(k));
2010-02-26 22:06:00 +05:00
mesh3d.CalcSurfacesOfNode();
mesh3d.FindOpenElements(k);
if (!mesh3d.GetNOpenElements())
continue;
2010-02-26 22:06:00 +05:00
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;
2016-10-11 17:10:36 +05:00
2010-02-26 22:06:00 +05:00
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());
2016-04-20 11:11:59 +05:00
for (int qstep = 0; qstep <= 3; qstep++)
// for (int qstep = 0; qstep <= 0; qstep++) // for hex-filling
2010-02-26 22:06:00 +05:00
{
if (qstep == 0 && !mp.try_hexes) continue;
// cout << "openquads = " << mesh3d.HasOpenQuads() << endl;
2010-02-26 22:06:00 +05:00
if (mesh3d.HasOpenQuads())
{
string rulefile = ngdir;
const char ** rulep = NULL;
switch (qstep)
{
case 0:
2016-04-10 09:36:05 +05:00
rulefile = "/Users/joachim/gitlab/netgen/rules/hexa.rls";
2016-04-20 11:11:59 +05:00
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;
}
2016-04-20 11:11:59 +05:00
// 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)
2013-04-03 02:29:53 +06:00
for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
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);
2014-12-02 18:23:36 +05:00
/*
2013-04-03 02:29:53 +06:00
for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
meshing.AddPoint (mesh3d[pi], pi);
2014-12-02 18:23:36 +05:00
*/
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())
2018-03-09 03:47:58 +05:00
do
{
if (multithread.terminate)
2018-03-09 03:47:58 +05:00
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);
2019-07-09 13:39:16 +05:00
NgArray<int, PointIndex::BASE> glob2loc(mesh3d.GetNP());
glob2loc = -1;
2013-04-03 02:29:53 +06:00
for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
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);
2018-03-09 03:47:58 +05:00
for (ElementIndex ei = oldne; ei < mesh3d.GetNE(); ei++)
mesh3d[ei].SetIndex (k);
2018-03-09 03:47:58 +05:00
mesh3d.CalcSurfacesOfNode();
mesh3d.FindOpenElements(k);
2011-01-11 01:18:01 +05:00
// teterrpow = 2;
if (mesh3d.GetNOpenElements() != 0)
{
meshed = 0;
PrintMessage (5, mesh3d.GetNOpenElements(), " open faces found");
2011-07-25 17:33:19 +06:00
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;
2011-07-25 17:33:19 +06:00
case 'm': mesh3d.ImproveMesh(mp, OPT_REST); break;
}
}
mesh3d.FindOpenElements(k);
PrintMessage (3, "Call remove problem");
RemoveProblem (mesh3d, k);
mesh3d.FindOpenElements(k);
}
else
2018-03-09 03:47:58 +05:00
{
meshed = 1;
PrintMessage (1, "Success !");
2018-03-09 03:47:58 +05:00
}
}
while (!meshed);
PrintMessage (1, mesh3d.GetNP(), " points, ",
mesh3d.GetNE(), " elements");
}
mp.maxh = globmaxh;
MeshQuality3d (mesh3d);
return MESHING3_OK;
2009-01-13 04:40:13 +05:00
}
/*
MESHING3_RESULT MeshVolumeOld (MeshingParameters & mp, Mesh& mesh3d)
{
int i, k, oldne;
int meshed;
int cntsteps;
PlotStatistics3d * pstat;
if (globflags.GetNumFlag("silentflag", 1) <= 2)
pstat = new XPlotStatistics3d;
else
pstat = new TerminalPlotStatistics3d;
cntsteps = 0;
do
{
cntsteps++;
if (cntsteps > mp.maxoutersteps)
{
return MESHING3_OUTERSTEPSEXCEEDED;
}
int noldp = mesh3d.GetNP();
if ( (cntsteps == 1) && globflags.GetDefineFlag ("delaunay"))
{
cntsteps ++;
mesh3d.CalcSurfacesOfNode();
for (k = 1; k <= mesh3d.GetNDomains(); k++)
{
Meshing3 meshing(NULL, pstat);
mesh3d.FindOpenElements(k);
for (i = 1; i <= noldp; i++)
meshing.AddPoint (mesh3d.Point(i), i);
for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
{
if (mesh3d.OpenElement(i).GetIndex() == k)
meshing.AddBoundaryElement (mesh3d.OpenElement(i));
}
oldne = mesh3d.GetNE();
if (globflags.GetDefineFlag ("blockfill"))
{
if (!globflags.GetDefineFlag ("localh"))
meshing.BlockFill
(mesh3d, mp.h * globflags.GetNumFlag ("relblockfillh", 1));
else
meshing.BlockFillLocalH (mesh3d);
}
MeshingParameters mpd;
meshing.Delaunay (mesh3d, mpd);
for (i = oldne + 1; i <= mesh3d.GetNE(); i++)
mesh3d.VolumeElement(i).SetIndex (k);
}
}
noldp = mesh3d.GetNP();
mesh3d.CalcSurfacesOfNode();
mesh3d.FindOpenElements();
for (k = 1; k <= mesh3d.GetNDomains(); k++)
{
Meshing3 meshing(globflags.GetStringFlag ("rules3d", NULL), pstat);
Point3d pmin, pmax;
mesh3d.GetBox (pmin, pmax, k);
rot.SetCenter (Center (pmin, pmax));
for (i = 1; i <= noldp; i++)
meshing.AddPoint (mesh3d.Point(i), i);
for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
{
if (mesh3d.OpenElement(i).GetIndex() == k)
meshing.AddBoundaryElement (mesh3d.OpenElement(i));
}
oldne = mesh3d.GetNE();
if ( (cntsteps == 1) && globflags.GetDefineFlag ("blockfill"))
{
if (!globflags.GetDefineFlag ("localh"))
{
meshing.BlockFill
(mesh3d,
mp.h * globflags.GetNumFlag ("relblockfillh", 1));
}
else
{
meshing.BlockFillLocalH (mesh3d);
}
}
mp.giveuptol = int(globflags.GetNumFlag ("giveuptol", 15));
meshing.GenerateMesh (mesh3d, mp);
for (i = oldne + 1; i <= mesh3d.GetNE(); i++)
mesh3d.VolumeElement(i).SetIndex (k);
}
mesh3d.CalcSurfacesOfNode();
mesh3d.FindOpenElements();
teterrpow = 2;
if (mesh3d.GetNOpenElements() != 0)
{
meshed = 0;
(*mycout) << "Open elements found, old" << endl;
const char * optstr = "mcmcmcmcm";
int j;
for (j = 1; j <= strlen(optstr); j++)
switch (optstr[j-1])
{
case 'c': mesh3d.CombineImprove(); break;
case 'd': mesh3d.SplitImprove(); break;
case 's': mesh3d.SwapImprove(); break;
case 'm': mesh3d.ImproveMesh(2); break;
}
(*mycout) << "Call remove" << endl;
RemoveProblem (mesh3d);
(*mycout) << "Problem removed" << endl;
}
else
meshed = 1;
}
while (!meshed);
MeshQuality3d (mesh3d);
return MESHING3_OK;
}
*/
/*
MESHING3_RESULT MeshMixedVolume(MeshingParameters & mp, Mesh& mesh3d)
{
int i, j;
MESHING3_RESULT res;
Point3d pmin, pmax;
mp.giveuptol = 10;
mp.baseelnp = 4;
mp.starshapeclass = 100;
// TerminalPlotStatistics3d pstat;
Meshing3 meshing1("pyramids.rls");
for (i = 1; i <= mesh3d.GetNP(); i++)
meshing1.AddPoint (mesh3d.Point(i), i);
mesh3d.FindOpenElements();
for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
if (mesh3d.OpenElement(i).GetIndex() == 1)
meshing1.AddBoundaryElement (mesh3d.OpenElement(i));
res = meshing1.GenerateMesh (mesh3d, mp);
mesh3d.GetBox (pmin, pmax);
PrintMessage (1, "Mesh pyramids, res = ", res);
if (res)
exit (1);
for (i = 1; i <= mesh3d.GetNE(); i++)
mesh3d.VolumeElement(i).SetIndex (1);
// do delaunay
mp.baseelnp = 0;
mp.starshapeclass = 5;
Meshing3 meshing2(NULL);
for (i = 1; i <= mesh3d.GetNP(); i++)
meshing2.AddPoint (mesh3d.Point(i), i);
mesh3d.FindOpenElements();
for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
if (mesh3d.OpenElement(i).GetIndex() == 1)
meshing2.AddBoundaryElement (mesh3d.OpenElement(i));
MeshingParameters mpd;
meshing2.Delaunay (mesh3d, mpd);
for (i = 1; i <= mesh3d.GetNE(); i++)
mesh3d.VolumeElement(i).SetIndex (1);
mp.baseelnp = 0;
mp.giveuptol = 10;
for (int trials = 1; trials <= 50; trials++)
{
if (multithread.terminate)
return MESHING3_TERMINATE;
Meshing3 meshing3("tetra.rls");
for (i = 1; i <= mesh3d.GetNP(); i++)
meshing3.AddPoint (mesh3d.Point(i), i);
mesh3d.FindOpenElements();
for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
if (mesh3d.OpenElement(i).GetIndex() == 1)
meshing3.AddBoundaryElement (mesh3d.OpenElement(i));
if (trials > 1)
CheckSurfaceMesh2 (mesh3d);
res = meshing3.GenerateMesh (mesh3d, mp);
for (i = 1; i <= mesh3d.GetNE(); i++)
mesh3d.VolumeElement(i).SetIndex (1);
if (res == 0) break;
for (i = 1; i <= mesh3d.GetNE(); i++)
{
const Element & el = mesh3d.VolumeElement(i);
if (el.GetNP() != 4)
{
for (j = 1; j <= el.GetNP(); j++)
mesh3d.AddLockedPoint (el.PNum(j));
}
}
mesh3d.CalcSurfacesOfNode();
mesh3d.FindOpenElements();
MeshOptimize3d optmesh;
teterrpow = 2;
const char * optstr = "mcmcmcmcm";
for (j = 1; j <= strlen(optstr); j++)
switch (optstr[j-1])
{
case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break;
case 'd': optmesh.SplitImprove(mesh3d); break;
case 's': optmesh.SwapImprove(mesh3d); break;
case 'm': mesh3d.ImproveMesh(); break;
}
RemoveProblem (mesh3d);
}
PrintMessage (1, "Meshing tets, res = ", res);
if (res)
{
mesh3d.FindOpenElements();
2017-12-05 02:32:47 +05:00
PrintSysError (1, "Open elements: ", mesh3d.GetNOpenElements());
2009-01-13 04:40:13 +05:00
exit (1);
}
for (i = 1; i <= mesh3d.GetNE(); i++)
{
const Element & el = mesh3d.VolumeElement(i);
if (el.GetNP() != 4)
{
for (j = 1; j <= el.GetNP(); j++)
mesh3d.AddLockedPoint (el.PNum(j));
}
}
mesh3d.CalcSurfacesOfNode();
mesh3d.FindOpenElements();
MeshOptimize3d optmesh;
teterrpow = 2;
const char * optstr = "mcmcmcmcm";
for (j = 1; j <= strlen(optstr); j++)
switch (optstr[j-1])
{
case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break;
case 'd': optmesh.SplitImprove(mesh3d); break;
case 's': optmesh.SwapImprove(mesh3d); break;
case 'm': mesh3d.ImproveMesh(); break;
}
return MESHING3_OK;
}
*/
MESHING3_RESULT OptimizeVolume (const MeshingParameters & mp,
2009-01-13 04:40:13 +05:00
Mesh & mesh3d)
// const CSGeometry * geometry)
{
2019-01-31 22:41:20 +05:00
static Timer t("OptimizeVolume"); RegionTimer reg(t);
2009-01-13 04:40:13 +05:00
int i;
PrintMessage (1, "Volume Optimization");
/*
if (!mesh3d.PureTetMesh())
return MESHING3_OK;
*/
// (*mycout) << "optstring = " << mp.optimize3d << endl;
/*
const char * optstr = globflags.GetStringFlag ("optimize3d", "cmh");
int optsteps = int (globflags.GetNumFlag ("optsteps3d", 2));
*/
mesh3d.CalcSurfacesOfNode();
for (i = 1; i <= mp.optsteps3d; i++)
{
if (multithread.terminate)
break;
2011-07-25 17:33:19 +06:00
MeshOptimize3d optmesh(mp);
2009-01-13 04:40:13 +05:00
2011-01-11 01:18:01 +05:00
// teterrpow = mp.opterrpow;
2014-08-31 15:14:18 +06:00
// for (size_t j = 1; j <= strlen(mp.optimize3d); j++)
for (size_t j = 1; j <= mp.optimize3d.length(); j++)
2009-01-13 04:40:13 +05:00
{
if (multithread.terminate)
break;
switch (mp.optimize3d[j-1])
{
case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break;
case 'd': optmesh.SplitImprove(mesh3d); break;
case 's': optmesh.SwapImprove(mesh3d); break;
// case 'u': optmesh.SwapImproveSurface(mesh3d); break;
case 't': optmesh.SwapImprove2(mesh3d); break;
#ifdef SOLIDGEOM
case 'm': mesh3d.ImproveMesh(*geometry); break;
case 'M': mesh3d.ImproveMesh(*geometry); break;
#else
2011-07-25 17:33:19 +06:00
case 'm': mesh3d.ImproveMesh(mp); break;
case 'M': mesh3d.ImproveMesh(mp); break;
2009-01-13 04:40:13 +05:00
#endif
2011-07-25 17:33:19 +06:00
case 'j': mesh3d.ImproveMeshJacobian(mp); break;
2009-01-13 04:40:13 +05:00
}
}
// mesh3d.mglevels = 1;
2009-01-13 04:40:13 +05:00
MeshQuality3d (mesh3d);
}
return MESHING3_OK;
}
void RemoveIllegalElements (Mesh & mesh3d)
{
2019-01-31 22:41:20 +05:00
static Timer t("RemoveIllegalElements"); RegionTimer reg(t);
2009-01-13 04:40:13 +05:00
int it = 10;
int nillegal, oldn;
PrintMessage (1, "Remove Illegal Elements");
// return, if non-pure tet-mesh
/*
if (!mesh3d.PureTetMesh())
return;
*/
mesh3d.CalcSurfacesOfNode();
nillegal = mesh3d.MarkIllegalElements();
2011-07-25 17:33:19 +06:00
MeshingParameters dummymp;
MeshOptimize3d optmesh(dummymp);
2009-01-13 04:40:13 +05:00
while (nillegal && (it--) > 0)
{
if (multithread.terminate)
break;
PrintMessage (5, nillegal, " illegal tets");
2016-12-12 00:17:07 +05:00
optmesh.SplitImprove (mesh3d, OPT_LEGAL);
2009-01-13 04:40:13 +05:00
mesh3d.MarkIllegalElements(); // test
optmesh.SwapImprove (mesh3d, OPT_LEGAL);
mesh3d.MarkIllegalElements(); // test
optmesh.SwapImprove2 (mesh3d, OPT_LEGAL);
oldn = nillegal;
nillegal = mesh3d.MarkIllegalElements();
if (oldn != nillegal)
it = 10;
}
PrintMessage (5, nillegal, " illegal tets");
}
}