netgen/libsrc/meshing/meshfunc.cpp

765 lines
22 KiB
C++
Raw Normal View History

#include <set>
2009-01-13 04:40:13 +05:00
#include <mystdlib.h>
#include "meshing.hpp"
2021-06-16 17:04:06 +05:00
#include "debugging.hpp"
2009-01-13 04:40:13 +05:00
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
2021-06-16 17:04:06 +05:00
struct MeshingData
{
int domain;
// mesh for one domain (contains all adjacent surface elements)
2021-06-22 14:16:28 +05:00
unique_ptr<Mesh> mesh;
2021-06-16 17:04:06 +05:00
// maps from local (domain) mesh to global mesh
2021-06-16 17:04:06 +05:00
Array<PointIndex, PointIndex> pmap;
2021-07-16 14:30:11 +05:00
// Array<INDEX_2> connected_pairs;
2021-06-16 17:04:06 +05:00
MeshingParameters mp;
unique_ptr<Meshing3> meshing;
};
// extract surface meshes belonging to individual domains
Array<MeshingData> DivideMesh(Mesh & mesh, const MeshingParameters & mp)
2021-06-16 17:04:06 +05:00
{
static Timer timer("DivideMesh"); RegionTimer rt(timer);
Array<MeshingData> ret;
auto num_domains = mesh.GetNDomains();
if(num_domains==1 || mp.only3D_domain_nr)
{
ret.SetSize(1);
// no need to divide mesh, just fill in meshing data
ret[0].domain = 1;
if(mp.only3D_domain_nr)
ret[0].domain = mp.only3D_domain_nr;
ret[0].mesh.reset(&mesh); // careful, this unique_ptr must not delete &mesh! (it will be released in MergeMeshes after meshing)
ret[0].mp = mp;
return ret;
}
ret.SetSize(num_domains);
2021-06-16 17:04:06 +05:00
Array<Array<PointIndex, PointIndex>> ipmap;
ipmap.SetSize(num_domains);
2023-08-05 15:01:01 +05:00
// auto dim = mesh.GetDimension();
2021-06-16 17:04:06 +05:00
auto num_points = mesh.GetNP();
auto num_facedescriptors = mesh.GetNFD();
for(auto i : Range(ret))
{
auto & md = ret[i];
md.domain = i+1;
md.mp = mp;
md.mp.maxh = min2 (mp.maxh, mesh.MaxHDomain(md.domain));
2021-06-22 14:16:28 +05:00
ret[i].mesh = make_unique<Mesh>();
auto & m = *ret[i].mesh;
2021-06-16 17:04:06 +05:00
m.SetLocalH(mesh.GetLocalH());
2021-06-16 17:04:06 +05:00
ipmap[i].SetSize(num_points);
ipmap[i] = PointIndex::INVALID;
m.SetDimension( mesh.GetDimension() );
2021-07-16 14:30:11 +05:00
m.SetGeometry( mesh.GetGeometry() );
2021-06-16 17:04:06 +05:00
for(auto i : Range(1, num_facedescriptors+1))
m.AddFaceDescriptor( mesh.GetFaceDescriptor(i) );
}
// mark used points for each domain, add surface elements (with wrong point numbers) to domain mesh
for(const auto & sel : mesh.SurfaceElements())
{
const auto & fd = mesh.GetFaceDescriptor(sel.GetIndex());
int dom_in = fd.DomainIn();
int dom_out = fd.DomainOut();
for( auto dom : {dom_in, dom_out} )
{
if(dom==0)
continue;
2021-06-22 14:16:28 +05:00
auto & sels = ret[dom-1].mesh->SurfaceElements();
2021-06-16 17:04:06 +05:00
for(auto pi : sel.PNums())
ipmap[dom-1][pi] = 1;
sels.Append(sel);
}
}
// mark used points for already existing volume elements, add them (with wrong point numbers) to domain mesh
for(const auto & el : mesh.VolumeElements())
{
auto dom = el.GetIndex();
auto & els = ret[dom-1].mesh->VolumeElements();
for(auto pi : el.PNums())
ipmap[dom-1][pi] = 1;
els.Append(el);
}
2021-06-22 14:35:44 +05:00
// mark locked/fixed points for each domain TODO: domain bounding box to add only relevant points?
for(auto pi : mesh.LockedPoints())
for(auto i : Range(ret))
ipmap[i][pi] = 2;
2021-06-22 14:35:44 +05:00
2021-06-16 17:04:06 +05:00
// add used points to domain mesh, build point mapping
for(auto i : Range(ret))
{
2021-06-22 14:16:28 +05:00
auto & m = *ret[i].mesh;
2021-06-16 17:04:06 +05:00
auto & pmap = ret[i].pmap;
for(auto pi : Range(ipmap[i]))
if(ipmap[i][pi])
{
const auto& mp = mesh[pi];
auto pi_new = m.AddPoint( mp, mp.GetLayer(), mp.Type() );
if(ipmap[i][pi] == 2)
mesh.AddLockedPoint(pi_new);
2021-06-16 17:04:06 +05:00
ipmap[i][pi] = pi_new;
pmap.Append( pi );
}
}
// add segments
for(auto i : Range(ret))
{
auto & imap = ipmap[i];
auto & m = *ret[i].mesh;
for(auto seg : mesh.LineSegments())
if(imap[seg[0]].IsValid() && imap[seg[1]].IsValid())
{
seg[0] = imap[seg[0]];
seg[1] = imap[seg[1]];
m.AddSegment(seg);
}
}
auto & identifications = mesh.GetIdentifications();
2021-06-16 17:04:06 +05:00
for(auto i : Range(ret))
{
2021-06-22 14:16:28 +05:00
auto & m = *ret[i].mesh;
2021-07-16 14:30:11 +05:00
auto & imap = ipmap[i];
auto nmax = identifications.GetMaxNr ();
auto & m_ident = m.GetIdentifications();
2021-06-16 17:04:06 +05:00
for (auto & sel : m.SurfaceElements())
for(auto & pi : sel.PNums())
pi = imap[pi];
for (auto & el : m.VolumeElements())
for(auto & pi : el.PNums())
pi = imap[pi];
2021-07-16 14:30:11 +05:00
for(auto n : Range(1,nmax+1))
{
NgArray<INDEX_2> pairs;
identifications.GetPairs(n, pairs);
for(auto pair : pairs)
{
auto pi0 = imap[pair[0]];
auto pi1 = imap[pair[1]];
if(!pi0.IsValid() || !pi1.IsValid())
continue;
m_ident.Add(pi0, pi1, n);
2021-07-16 14:30:11 +05:00
}
m_ident.SetType( n, identifications.GetType(n) );
2021-07-16 14:30:11 +05:00
}
2021-06-16 17:04:06 +05:00
}
return ret;
}
// Add between identified surface elements (only consider closesurface identifications)
void FillCloseSurface( MeshingData & md)
{
static Timer timer("FillCloseSurface"); RegionTimer rtimer(timer);
auto & mesh = *md.mesh;
auto & identifications = mesh.GetIdentifications();
auto nmax = identifications.GetMaxNr();
bool have_closesurfaces = false;
for(auto i : Range(1,nmax+1))
if(identifications.GetType(i) == Identifications::CLOSESURFACES)
have_closesurfaces = true;
if(!have_closesurfaces)
return;
NgArray<int, PointIndex::BASE> map;
std::set<std::tuple<int,int,int>> hex_faces;
for(auto identnr : Range(1,nmax+1))
{
if(identifications.GetType(identnr) != Identifications::CLOSESURFACES)
continue;
identifications.GetMap(identnr, map);
mesh.FindOpenElements(md.domain);
for(auto & sel : mesh.OpenElements())
{
// For quads: check if this open element is already closed by a hex
2022-09-07 04:43:32 +05:00
// this happens when we have identifications in two directions
if(sel.GetNP() == 4)
{
Element2d face = sel;
face.NormalizeNumbering();
if(hex_faces.count({face[0], face[1], face[2]}))
continue;
}
bool is_mapped = true;
for(auto pi : sel.PNums())
if(!PointIndex(map[pi]).IsValid())
is_mapped = false;
if(!is_mapped)
continue;
// insert prism/hex
auto np = sel.GetNP();
Element el(2*np);
std::set<int> pis;
for(auto i : Range(np))
{
el[i] = sel[i];
el[i+np] = map[sel[i]];
pis.insert(sel[i]);
pis.insert(map[sel[i]]);
}
// degenerate element (mapped element onto itself, might happen for surface elements connecting two identified faces)
if(pis.size() < 2*np)
continue;
// check if new element is inside current domain
auto p0 = mesh[sel[0]];
Vec<3> n = -Cross(mesh[sel[1]] - p0, mesh[sel[2]] - p0 );
if(n*(mesh[el[np]]-p0) < 0.0)
continue;
el.SetIndex(md.domain);
mesh.AddVolumeElement(el);
if(el.NP()==8)
{
// remember all adjacent faces of the new hex (to skip corresponding openelements accordingly)
for(auto facei : Range(1,7))
{
Element2d face;
el.GetFace(facei, face);
face.NormalizeNumbering();
hex_faces.insert({face[0], face[1], face[2]});
}
}
}
}
}
2021-06-16 17:04:06 +05:00
void CloseOpenQuads( MeshingData & md)
{
static Timer t("CloseOpenQuads"); RegionTimer rt(t);
2021-06-22 14:16:28 +05:00
auto & mesh = *md.mesh;
2021-06-16 17:04:06 +05:00
auto domain = md.domain;
MeshingParameters & mp = md.mp;
int oldne;
if (multithread.terminate)
return;
mesh.CalcSurfacesOfNode();
2021-06-16 18:35:50 +05:00
mesh.FindOpenElements(domain);
2021-06-16 17:04:06 +05:00
if (!mesh.GetNOpenElements())
return;
for (int qstep = 0; qstep <= 3; qstep++)
{
if (qstep == 0 && !mp.try_hexes) continue;
if (mesh.HasOpenQuads())
{
string rulefile = ngdir;
const char ** rulep = NULL;
switch (qstep)
{
case 0:
rulep = hexrules;
break;
case 1:
rulep = prismrules2;
break;
case 2: // connect pyramid to triangle
rulep = pyramidrules2;
break;
case 3: // connect to vis-a-vis point
rulep = pyramidrules;
break;
}
Meshing3 meshing(rulep);
MeshingParameters mpquad = mp;
mpquad.giveuptol = mp.giveuptolopenquads;
2021-06-16 17:04:06 +05:00
mpquad.baseelnp = 4;
mpquad.starshapeclass = 1000;
mpquad.check_impossible = qstep == 1; // for prisms only (air domain in trafo)
for (PointIndex pi : mesh.Points().Range())
meshing.AddPoint (mesh[pi], pi);
2021-07-16 14:30:11 +05:00
NgArray<INDEX_2> connectednodes;
for (int nr = 1; nr <= mesh.GetIdentifications().GetMaxNr(); nr++)
if (mesh.GetIdentifications().GetType(nr) != Identifications::PERIODIC)
{
mesh.GetIdentifications().GetPairs (nr, connectednodes);
for (auto pair : connectednodes)
{
2021-07-16 14:30:11 +05:00
meshing.AddConnectedPair (pair);
meshing.AddConnectedPair ({pair[1], pair[0]});
}
2021-07-16 14:30:11 +05:00
}
2021-06-16 17:04:06 +05:00
for (int i = 1; i <= mesh.GetNOpenElements(); i++)
{
Element2d hel = mesh.OpenElement(i);
meshing.AddBoundaryElement (hel);
}
oldne = mesh.GetNE();
meshing.GenerateMesh (mesh, mpquad);
for (int i = oldne + 1; i <= mesh.GetNE(); i++)
mesh.VolumeElement(i).SetIndex (domain);
(*testout)
<< "mesh has " << mesh.GetNE() << " prism/pyramid elements" << endl;
2021-06-16 18:35:50 +05:00
mesh.FindOpenElements(domain);
2021-06-16 17:04:06 +05:00
}
}
if (mesh.HasOpenQuads())
{
if(debugparam.write_mesh_on_error)
md.mesh->Save("open_quads_"+ToString(md.domain)+".vol.gz");
2021-06-16 17:04:06 +05:00
PrintSysError ("mesh has still open quads");
throw NgException ("Stop meshing since too many attempts");
// return MESHING3_GIVEUP;
}
}
void MeshDomain( MeshingData & md)
{
2021-06-22 14:16:28 +05:00
auto & mesh = *md.mesh;
auto domain = md.domain;
MeshingParameters & mp = md.mp;
mesh.CalcSurfacesOfNode();
mesh.FindOpenElements(md.domain);
md.meshing = make_unique<Meshing3>(nullptr);
for (PointIndex pi : mesh.Points().Range())
md.meshing->AddPoint (mesh[pi], pi);
for (int i = 1; i <= mesh.GetNOpenElements(); i++)
md.meshing->AddBoundaryElement (mesh.OpenElement(i));
2021-06-16 17:04:06 +05:00
if (mp.delaunay && mesh.GetNOpenElements())
{
int oldne = mesh.GetNE();
md.meshing->Delaunay (mesh, domain, mp);
for (int i = oldne + 1; i <= mesh.GetNE(); i++)
mesh.VolumeElement(i).SetIndex (domain);
PrintMessage (3, mesh.GetNP(), " points, ",
mesh.GetNE(), " elements");
mesh.FindOpenElements(domain);
2021-06-16 17:04:06 +05:00
}
Box<3> domain_bbox( Box<3>::EMPTY_BOX );
for (auto & sel : mesh.SurfaceElements())
{
if (sel.IsDeleted() ) continue;
for (auto pi : sel.PNums())
domain_bbox.Add (mesh[pi]);
}
domain_bbox.Increase (0.01 * domain_bbox.Diam());
2021-06-16 17:04:06 +05:00
int cntsteps = 0;
int meshed;
if (mesh.GetNOpenElements())
do
{
if (multithread.terminate)
break;
mesh.FindOpenElements(domain);
2021-06-16 17:04:06 +05:00
PrintMessage (5, mesh.GetNOpenElements(), " open faces");
2021-09-21 17:39:35 +05:00
// GetOpenElements( mesh, domain )->Save("open_"+ToString(cntsteps)+".vol");
2021-06-16 17:04:06 +05:00
cntsteps++;
if (cntsteps > mp.maxoutersteps)
{
if(debugparam.write_mesh_on_error)
md.mesh->Save("meshing_error_domain_"+ToString(md.domain)+".vol.gz");
throw NgException ("Stop meshing since too many attempts in domain " + ToString(md.domain));
}
2021-06-16 17:04:06 +05:00
PrintMessage (1, "start tetmeshing");
Meshing3 meshing(tetrules);
Array<PointIndex, PointIndex> glob2loc(mesh.GetNP());
glob2loc = PointIndex::INVALID;
for (PointIndex pi : mesh.Points().Range())
if (domain_bbox.IsIn (mesh[pi]))
glob2loc[pi] = meshing.AddPoint (mesh[pi], pi);
2021-06-16 17:04:06 +05:00
for (auto sel : mesh.OpenElements() )
{
for(auto & pi : sel.PNums())
pi = glob2loc[pi];
meshing.AddBoundaryElement (sel);
}
int oldne = mesh.GetNE();
mp.giveuptol = 15 + 10 * cntsteps;
mp.sloppy = 5;
meshing.GenerateMesh (mesh, mp);
for (ElementIndex ei = oldne; ei < mesh.GetNE(); ei++)
mesh[ei].SetIndex (domain);
mesh.CalcSurfacesOfNode();
mesh.FindOpenElements(domain);
2021-06-16 17:04:06 +05:00
// teterrpow = 2;
if (mesh.GetNOpenElements() != 0)
{
meshed = 0;
PrintMessage (5, mesh.GetNOpenElements(), " open faces found");
2024-02-23 21:41:15 +05:00
MeshOptimize3d optmesh(mesh, mp, OPT_REST);
2021-06-16 17:04:06 +05:00
const char * optstr = "mcmstmcmstmcmstmcm";
for (size_t j = 1; j <= strlen(optstr); j++)
{
mesh.FindOpenElements();
2021-06-16 17:04:06 +05:00
mesh.CalcSurfacesOfNode();
mesh.FreeOpenElementsEnvironment(2);
mesh.CalcSurfacesOfNode();
switch (optstr[j-1])
{
2024-02-23 21:41:15 +05:00
case 'c': optmesh.CombineImprove(); break;
case 'd': optmesh.SplitImprove(); break;
case 's': optmesh.SwapImprove(); break;
case 't': optmesh.SwapImprove2(); break;
case 'm': optmesh.ImproveMesh(); break;
2021-06-16 17:04:06 +05:00
}
}
mesh.FindOpenElements(domain);
2021-06-16 17:04:06 +05:00
PrintMessage (3, "Call remove problem");
RemoveProblem (mesh, domain);
mesh.FindOpenElements(domain);
2021-06-16 17:04:06 +05:00
}
else
{
meshed = 1;
PrintMessage (1, "Success !");
}
}
while (!meshed);
{
PrintMessage (3, "Check subdomain ", domain, " / ", mesh.GetNDomains());
mesh.FindOpenElements(domain);
2021-06-16 17:04:06 +05:00
bool res = (mesh.CheckConsistentBoundary() != 0);
if (res)
{
if(debugparam.write_mesh_on_error)
md.mesh->Save("inconsistent_surface_domain_"+ToString(md.domain)+".vol.gz");
PrintError ("Surface mesh not consistent");
throw NgException ("Stop meshing since surface mesh not consistent");
}
}
RemoveIllegalElements (mesh);
}
2021-06-16 17:04:06 +05:00
void MergeMeshes( Mesh & mesh, Array<MeshingData> & md )
{
// todo: optimize: count elements, alloc all memory, copy vol elements in parallel
static Timer t("MergeMeshes"); RegionTimer rt(t);
if(md.Size()==1)
{
// assume that mesh was never divided, no need to do anything
if(&mesh != md[0].mesh.get())
throw Exception("Illegal Mesh pointer in MeshingData");
md[0].mesh.release();
return;
}
mesh.VolumeElements().DeleteAll();
2021-06-16 17:04:06 +05:00
for(auto & m_ : md)
{
auto first_new_pi = m_.pmap.Range().Next();
2021-06-22 14:16:28 +05:00
auto & m = *m_.mesh;
2021-06-16 17:04:06 +05:00
Array<PointIndex, PointIndex> pmap(m.Points().Size());
for(auto pi : Range(PointIndex(PointIndex::BASE), first_new_pi))
pmap[pi] = m_.pmap[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]];
el.SetIndex(m_.domain);
mesh.AddVolumeElement(el);
}
}
}
void MergeMeshes( Mesh & mesh, FlatArray<Mesh> meshes, PointIndex first_new_pi )
{
// todo: optimize: count elements, alloc all memory, copy vol elements in parallel
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;
2009-01-13 04:40:13 +05:00
for (auto pi : Range(first_new_pi, m.Points().Range().Next()))
pmap[pi] = mesh.AddPoint(m[pi]);
2009-01-13 04:40:13 +05:00
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 & mp, Mesh& mesh3d)
2021-06-16 17:04:06 +05:00
{
static Timer t("MeshVolume"); RegionTimer reg(t);
mesh3d.Compress();
if(mesh3d.GetNDomains()==0)
return MESHING3_OK;
2021-07-16 21:59:18 +05:00
if (!mesh3d.HasLocalHFunction())
mesh3d.CalcLocalH(mp.grading);
2021-06-16 17:04:06 +05:00
auto md = DivideMesh(mesh3d, mp);
try
{
2021-06-16 17:04:06 +05:00
ParallelFor( md.Range(), [&](int i)
{
if (mp.checkoverlappingboundary)
if (md[i].mesh->CheckOverlappingBoundary())
{
if(debugparam.write_mesh_on_error)
md[i].mesh->Save("overlapping_mesh_domain_"+ToString(md[i].domain)+".vol.gz");
throw NgException ("Stop meshing since boundary mesh is overlapping");
}
if(md[i].mesh->GetGeometry()->GetGeomType() == Mesh::GEOM_OCC)
FillCloseSurface( md[i] );
2021-06-16 17:04:06 +05:00
CloseOpenQuads( md[i] );
MeshDomain(md[i]);
}, md.Size());
}
catch(...)
{
MergeMeshes(mesh3d, md);
return MESHING3_GIVEUP;
}
2021-06-16 17:04:06 +05:00
MergeMeshes(mesh3d, md);
MeshQuality3d (mesh3d);
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);
2023-09-05 20:13:29 +05:00
#ifndef EMSCRIPTEN
RegionTaskManager rtm(mp.parallel_meshing ? mp.nthreads : 0);
2023-09-05 20:13:29 +05:00
#endif // EMSCRIPTEN
const char* savetask = multithread.task;
multithread.task = "Optimize Volume";
2019-01-31 22:41:20 +05:00
2023-08-05 15:01:01 +05:00
// int i;
2009-01-13 04:40:13 +05:00
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();
MeshOptimize3d optmesh(mesh3d, mp);
// optimize only bad elements first
optmesh.SetMinBadness(1000.);
bool do_split = mp.optimize3d.find('d') != string::npos;
bool do_swap = mp.optimize3d.find('s') != string::npos;
bool do_swap2 = mp.optimize3d.find('t') != string::npos;
for(auto i : Range(mp.optsteps3d))
{
auto [total_badness, max_badness, bad_els] = optmesh.UpdateBadness();
if(bad_els==0) break;
if(do_split) optmesh.SplitImprove();
if(do_swap) optmesh.SwapImprove();
if(do_swap2) optmesh.SwapImprove2();
}
// Now optimize all elements
optmesh.SetMinBadness(0);
for (auto i : Range(mp.optsteps3d))
2009-01-13 04:40:13 +05:00
{
if (multithread.terminate)
break;
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 (auto j : Range(mp.optimize3d.size()))
2009-01-13 04:40:13 +05:00
{
multithread.percent = 100.* (double(j)/mp.optimize3d.size() + i)/mp.optsteps3d;
2009-01-13 04:40:13 +05:00
if (multithread.terminate)
break;
switch (mp.optimize3d[j])
2009-01-13 04:40:13 +05:00
{
2024-02-23 21:41:15 +05:00
case 'c':
optmesh.SetGoal(OPT_REST);
optmesh.CombineImprove();
optmesh.SetGoal(OPT_QUALITY);
break;
case 'd': optmesh.SplitImprove(); break;
case 'D': optmesh.SplitImprove2(); break;
case 's': optmesh.SwapImprove(); break;
2009-01-13 04:40:13 +05:00
// case 'u': optmesh.SwapImproveSurface(mesh3d); break;
2024-02-23 21:41:15 +05:00
case 't': optmesh.SwapImprove2(); break;
2009-01-13 04:40:13 +05:00
#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);
}
multithread.task = savetask;
2009-01-13 04:40:13 +05:00
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;
// return, if non-pure tet-mesh
/*
if (!mesh3d.PureTetMesh())
return;
*/
mesh3d.CalcSurfacesOfNode();
nillegal = mesh3d.MarkIllegalElements();
if(nillegal)
PrintMessage (1, "Remove Illegal Elements");
2009-01-13 04:40:13 +05:00
2011-07-25 17:33:19 +06:00
MeshingParameters dummymp;
2024-02-23 21:41:15 +05:00
MeshOptimize3d optmesh(mesh3d, dummymp, OPT_LEGAL);
2009-01-13 04:40:13 +05:00
while (nillegal && (it--) > 0)
{
if (multithread.terminate)
break;
PrintMessage (5, nillegal, " illegal tets");
2024-02-23 21:41:15 +05:00
optmesh.SplitImprove ();
2009-01-13 04:40:13 +05:00
mesh3d.MarkIllegalElements(); // test
2024-02-23 21:41:15 +05:00
optmesh.SwapImprove ();
2009-01-13 04:40:13 +05:00
mesh3d.MarkIllegalElements(); // test
2024-02-23 21:41:15 +05:00
optmesh.SwapImprove2 ();
2009-01-13 04:40:13 +05:00
oldn = nillegal;
nillegal = mesh3d.MarkIllegalElements();
if (oldn != nillegal)
it = 10;
}
PrintMessage (5, nillegal, " illegal tets");
}
}