mirror of
https://github.com/NGSolve/netgen.git
synced 2025-01-11 21:50:34 +05:00
Merge branch 'parallel_meshing' into test
This commit is contained in:
commit
82a67b27b9
@ -1545,7 +1545,7 @@ namespace netgen
|
|||||||
|
|
||||||
PrintMessage (1, "Delaunay meshing");
|
PrintMessage (1, "Delaunay meshing");
|
||||||
PrintMessage (3, "number of points: ", mesh.GetNP());
|
PrintMessage (3, "number of points: ", mesh.GetNP());
|
||||||
PushStatus ("Delaunay meshing");
|
// PushStatus ("Delaunay meshing");
|
||||||
|
|
||||||
|
|
||||||
NgArray<DelaunayTet> tempels;
|
NgArray<DelaunayTet> tempels;
|
||||||
@ -1675,6 +1675,6 @@ namespace netgen
|
|||||||
mesh.FindOpenElements(domainnr);
|
mesh.FindOpenElements(domainnr);
|
||||||
|
|
||||||
mesh.Compress();
|
mesh.Compress();
|
||||||
PopStatus ();
|
// PopStatus ();
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -130,6 +130,54 @@ namespace netgen
|
|||||||
return lh;
|
return lh;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
unique_ptr<LocalH> LocalH :: Copy( const Box<3> & bbox )
|
||||||
|
{
|
||||||
|
static Timer t("LocalH::Copy with bounding box"); RegionTimer rt(t);
|
||||||
|
auto lh = make_unique<LocalH>(boundingbox, grading, dimension);
|
||||||
|
std::map<GradingBox*, GradingBox*> mapping;
|
||||||
|
lh->boxes.SetAllocSize(boxes.Size());
|
||||||
|
|
||||||
|
for(auto i : boxes.Range())
|
||||||
|
{
|
||||||
|
auto & b = *boxes[i];
|
||||||
|
auto h = b.H2();
|
||||||
|
Vec<3> vh = {h,h,h};
|
||||||
|
Box<3> box( b.PMid() - vh, b.PMid() + vh);
|
||||||
|
if(!box.Intersect(bbox))
|
||||||
|
continue;
|
||||||
|
lh->boxes.Append(new GradingBox());
|
||||||
|
auto & bnew = *lh->boxes.Last();
|
||||||
|
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 & b = *boxes[i];
|
||||||
|
if(mapping.count(&b)==0)
|
||||||
|
continue;
|
||||||
|
|
||||||
|
auto & bnew = *mapping[&b];
|
||||||
|
for(auto k : Range(8))
|
||||||
|
{
|
||||||
|
if(b.childs[k] && mapping.count(b.childs[k]))
|
||||||
|
bnew.childs[k] = mapping[b.childs[k]];
|
||||||
|
}
|
||||||
|
|
||||||
|
if(b.father && mapping.count(b.father))
|
||||||
|
bnew.father = mapping[b.father];
|
||||||
|
}
|
||||||
|
|
||||||
|
lh->root = mapping[root];
|
||||||
|
return lh;
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
void LocalH :: Delete ()
|
void LocalH :: Delete ()
|
||||||
{
|
{
|
||||||
root->DeleteChilds();
|
root->DeleteChilds();
|
||||||
|
@ -99,6 +99,7 @@ namespace netgen
|
|||||||
~LocalH();
|
~LocalH();
|
||||||
///
|
///
|
||||||
unique_ptr<LocalH> Copy();
|
unique_ptr<LocalH> Copy();
|
||||||
|
unique_ptr<LocalH> Copy( const Box<3> & bbox );
|
||||||
///
|
///
|
||||||
void Delete();
|
void Delete();
|
||||||
///
|
///
|
||||||
@ -191,6 +192,8 @@ namespace netgen
|
|||||||
///
|
///
|
||||||
void ConvexifyRec (GradingBox * box);
|
void ConvexifyRec (GradingBox * box);
|
||||||
|
|
||||||
|
unique_ptr<LocalH> CopyRec( const Box<3> & bbox, GradingBox * current );
|
||||||
|
|
||||||
friend ostream & operator<< (ostream & ost, const LocalH & loch);
|
friend ostream & operator<< (ostream & ost, const LocalH & loch);
|
||||||
};
|
};
|
||||||
|
|
||||||
|
@ -21,7 +21,7 @@ namespace netgen
|
|||||||
// maps from local (domain) mesh to global mesh
|
// maps from local (domain) mesh to global mesh
|
||||||
Array<PointIndex, PointIndex> pmap;
|
Array<PointIndex, PointIndex> pmap;
|
||||||
|
|
||||||
Array<INDEX_2> connected_pairs;
|
// Array<INDEX_2> connected_pairs;
|
||||||
|
|
||||||
MeshingParameters mp;
|
MeshingParameters mp;
|
||||||
|
|
||||||
@ -36,9 +36,9 @@ namespace netgen
|
|||||||
Array<MeshingData> ret;
|
Array<MeshingData> ret;
|
||||||
auto num_domains = mesh.GetNDomains();
|
auto num_domains = mesh.GetNDomains();
|
||||||
|
|
||||||
ret.SetSize(num_domains);
|
|
||||||
if(num_domains==1 || mp.only3D_domain_nr)
|
if(num_domains==1 || mp.only3D_domain_nr)
|
||||||
{
|
{
|
||||||
|
ret.SetSize(1);
|
||||||
// no need to divide mesh, just fill in meshing data
|
// no need to divide mesh, just fill in meshing data
|
||||||
ret[0].domain = 1;
|
ret[0].domain = 1;
|
||||||
if(mp.only3D_domain_nr)
|
if(mp.only3D_domain_nr)
|
||||||
@ -48,6 +48,7 @@ namespace netgen
|
|||||||
ret[0].mp = mp;
|
ret[0].mp = mp;
|
||||||
return ret;
|
return ret;
|
||||||
}
|
}
|
||||||
|
ret.SetSize(num_domains);
|
||||||
|
|
||||||
Array<Array<PointIndex, PointIndex>> ipmap;
|
Array<Array<PointIndex, PointIndex>> ipmap;
|
||||||
ipmap.SetSize(num_domains);
|
ipmap.SetSize(num_domains);
|
||||||
@ -55,8 +56,6 @@ namespace netgen
|
|||||||
auto num_points = mesh.GetNP();
|
auto num_points = mesh.GetNP();
|
||||||
auto num_facedescriptors = mesh.GetNFD();
|
auto num_facedescriptors = mesh.GetNFD();
|
||||||
|
|
||||||
auto & identifications = mesh.GetIdentifications();
|
|
||||||
|
|
||||||
for(auto i : Range(ret))
|
for(auto i : Range(ret))
|
||||||
{
|
{
|
||||||
auto & md = ret[i];
|
auto & md = ret[i];
|
||||||
@ -73,6 +72,7 @@ namespace netgen
|
|||||||
ipmap[i].SetSize(num_points);
|
ipmap[i].SetSize(num_points);
|
||||||
ipmap[i] = PointIndex::INVALID;
|
ipmap[i] = PointIndex::INVALID;
|
||||||
m.SetDimension( mesh.GetDimension() );
|
m.SetDimension( mesh.GetDimension() );
|
||||||
|
m.SetGeometry( mesh.GetGeometry() );
|
||||||
|
|
||||||
for(auto i : Range(1, num_facedescriptors+1))
|
for(auto i : Range(1, num_facedescriptors+1))
|
||||||
m.AddFaceDescriptor( mesh.GetFaceDescriptor(i) );
|
m.AddFaceDescriptor( mesh.GetFaceDescriptor(i) );
|
||||||
@ -117,28 +117,51 @@ namespace netgen
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
NgArray<INDEX_2> connectednodes;
|
// add segmetns
|
||||||
for(auto i : Range(ret))
|
for(auto i : Range(ret))
|
||||||
{
|
{
|
||||||
auto & imap = ipmap[i];
|
auto & imap = ipmap[i];
|
||||||
auto & m = *ret[i].mesh;
|
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();
|
||||||
|
|
||||||
|
for(auto i : Range(ret))
|
||||||
|
{
|
||||||
|
auto & m = *ret[i].mesh;
|
||||||
|
auto & imap = ipmap[i];
|
||||||
|
auto nmax = identifications.GetMaxNr ();
|
||||||
|
auto & m_ident = m.GetIdentifications();
|
||||||
|
|
||||||
for (auto & sel : m.SurfaceElements())
|
for (auto & sel : m.SurfaceElements())
|
||||||
for(auto & pi : sel.PNums())
|
for(auto & pi : sel.PNums())
|
||||||
pi = imap[pi];
|
pi = imap[pi];
|
||||||
|
|
||||||
for (int nr = 1; nr <= identifications.GetMaxNr(); nr++)
|
for(auto n : Range(1,nmax+1))
|
||||||
if (identifications.GetType(nr) != Identifications::PERIODIC)
|
{
|
||||||
{
|
NgArray<INDEX_2> pairs;
|
||||||
identifications.GetPairs (nr, connectednodes);
|
identifications.GetPairs(n, pairs);
|
||||||
for (auto pair : connectednodes)
|
|
||||||
{
|
for(auto pair : pairs)
|
||||||
auto pi0 = pair[0];
|
{
|
||||||
auto pi1 = pair[1];
|
auto pi0 = imap[pair[0]];
|
||||||
if(imap[pi0].IsValid() && imap[pi1].IsValid())
|
auto pi1 = imap[pair[1]];
|
||||||
ret[i].connected_pairs.Append({imap[pi0], imap[pi1]});
|
if(!pi0.IsValid() || !pi1.IsValid())
|
||||||
}
|
continue;
|
||||||
}
|
|
||||||
|
if(pi1<pi0)
|
||||||
|
Swap(pi0,pi1);
|
||||||
|
m_ident.Add(pi0, pi1, n);
|
||||||
|
}
|
||||||
|
m_ident.SetType( n, identifications.GetType(n) );
|
||||||
|
}
|
||||||
}
|
}
|
||||||
return ret;
|
return ret;
|
||||||
}
|
}
|
||||||
@ -197,8 +220,16 @@ namespace netgen
|
|||||||
for (PointIndex pi : mesh.Points().Range())
|
for (PointIndex pi : mesh.Points().Range())
|
||||||
meshing.AddPoint (mesh[pi], pi);
|
meshing.AddPoint (mesh[pi], pi);
|
||||||
|
|
||||||
for (auto pair : md.connected_pairs)
|
NgArray<INDEX_2> connectednodes;
|
||||||
meshing.AddConnectedPair (pair);
|
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)
|
||||||
|
meshing.AddConnectedPair (pair);
|
||||||
|
}
|
||||||
|
// for (auto pair : md.connected_pairs)
|
||||||
|
// meshing.AddConnectedPair (pair);
|
||||||
|
|
||||||
for (int i = 1; i <= mesh.GetNOpenElements(); i++)
|
for (int i = 1; i <= mesh.GetNOpenElements(); i++)
|
||||||
{
|
{
|
||||||
@ -229,26 +260,6 @@ namespace netgen
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
void PrepareForBlockFillLocalH(MeshingData & md)
|
|
||||||
{
|
|
||||||
static Timer t("PrepareForBlockFillLocalH"); RegionTimer rt(t);
|
|
||||||
md.meshing = make_unique<Meshing3>(nullptr);
|
|
||||||
|
|
||||||
auto & mesh = *md.mesh;
|
|
||||||
|
|
||||||
mesh.CalcSurfacesOfNode();
|
|
||||||
mesh.FindOpenElements(md.domain);
|
|
||||||
|
|
||||||
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));
|
|
||||||
|
|
||||||
if (mesh.HasLocalHFunction())
|
|
||||||
md.meshing->PrepareBlockFillLocalH(mesh, md.mp);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
void MeshDomain( MeshingData & md)
|
void MeshDomain( MeshingData & md)
|
||||||
{
|
{
|
||||||
@ -256,6 +267,16 @@ namespace netgen
|
|||||||
auto domain = md.domain;
|
auto domain = md.domain;
|
||||||
MeshingParameters & mp = md.mp;
|
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));
|
||||||
|
|
||||||
if (mp.delaunay && mesh.GetNOpenElements())
|
if (mp.delaunay && mesh.GetNOpenElements())
|
||||||
{
|
{
|
||||||
int oldne = mesh.GetNE();
|
int oldne = mesh.GetNE();
|
||||||
@ -378,274 +399,6 @@ namespace netgen
|
|||||||
|
|
||||||
bool res = (mesh.CheckConsistentBoundary() != 0);
|
bool res = (mesh.CheckConsistentBoundary() != 0);
|
||||||
if (res)
|
if (res)
|
||||||
{
|
|
||||||
// mesh.Save("output.vol");
|
|
||||||
PrintError ("Surface mesh not consistent");
|
|
||||||
throw NgException ("Stop meshing since surface mesh not consistent");
|
|
||||||
}
|
|
||||||
}
|
|
||||||
// OptimizeVolume( md.mp, mesh );
|
|
||||||
}
|
|
||||||
|
|
||||||
void MeshDomain(Mesh & mesh3d, const MeshingParameters & c_mp, int k, const Identifications & identifications)
|
|
||||||
{
|
|
||||||
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 <= identifications.GetMaxNr(); nr++)
|
|
||||||
// if (identifications.GetType(nr) != Identifications::PERIODIC)
|
|
||||||
// {
|
|
||||||
// identifications.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");
|
PrintError ("Surface mesh not consistent");
|
||||||
throw NgException ("Stop meshing since surface mesh not consistent");
|
throw NgException ("Stop meshing since surface mesh not consistent");
|
||||||
@ -727,8 +480,7 @@ namespace netgen
|
|||||||
if(mesh3d.GetNDomains()==0)
|
if(mesh3d.GetNDomains()==0)
|
||||||
return MESHING3_OK;
|
return MESHING3_OK;
|
||||||
|
|
||||||
// localh function is built for each domain separately in blockfill ( more efficient )
|
if (!mesh3d.HasLocalHFunction())
|
||||||
if (!mesh3d.HasLocalHFunction() && !mp.blockfill)
|
|
||||||
mesh3d.CalcLocalH(mp.grading);
|
mesh3d.CalcLocalH(mp.grading);
|
||||||
|
|
||||||
auto md = DivideMesh(mesh3d, mp);
|
auto md = DivideMesh(mesh3d, mp);
|
||||||
@ -736,13 +488,6 @@ namespace netgen
|
|||||||
ParallelFor( md.Range(), [&](int i)
|
ParallelFor( md.Range(), [&](int i)
|
||||||
{
|
{
|
||||||
CloseOpenQuads( md[i] );
|
CloseOpenQuads( md[i] );
|
||||||
});
|
|
||||||
|
|
||||||
for(auto & md_ : md)
|
|
||||||
PrepareForBlockFillLocalH(md_);
|
|
||||||
|
|
||||||
ParallelFor( md.Range(), [&](int i)
|
|
||||||
{
|
|
||||||
MeshDomain(md[i]);
|
MeshDomain(md[i]);
|
||||||
});
|
});
|
||||||
|
|
||||||
|
@ -1094,10 +1094,11 @@ static int TestSameSide (const Point3d & p1, const Point3d & p2)
|
|||||||
*/
|
*/
|
||||||
|
|
||||||
|
|
||||||
void Meshing3 :: PrepareBlockFillLocalH (Mesh & mesh,
|
|
||||||
|
void Meshing3 :: BlockFillLocalH (Mesh & mesh,
|
||||||
const MeshingParameters & mp)
|
const MeshingParameters & mp)
|
||||||
{
|
{
|
||||||
static Timer t("Mesing3::PrepareBlockFillLocalH"); RegionTimer reg(t);
|
static Timer t("Mesing3::BlockFillLocalH"); RegionTimer reg(t);
|
||||||
|
|
||||||
double filldist = mp.filldist;
|
double filldist = mp.filldist;
|
||||||
|
|
||||||
@ -1106,8 +1107,11 @@ void Meshing3 :: PrepareBlockFillLocalH (Mesh & mesh,
|
|||||||
PrintMessage (3, "blockfill local h");
|
PrintMessage (3, "blockfill local h");
|
||||||
|
|
||||||
|
|
||||||
|
NgArray<Point<3> > npoints;
|
||||||
|
|
||||||
adfront -> CreateTrees();
|
adfront -> CreateTrees();
|
||||||
|
|
||||||
|
Box<3> bbox ( Box<3>::EMPTY_BOX );
|
||||||
double maxh = 0;
|
double maxh = 0;
|
||||||
|
|
||||||
for (int i = 1; i <= adfront->GetNF(); i++)
|
for (int i = 1; i <= adfront->GetNF(); i++)
|
||||||
@ -1120,15 +1124,28 @@ void Meshing3 :: PrepareBlockFillLocalH (Mesh & mesh,
|
|||||||
|
|
||||||
double hi = Dist (p1, p2);
|
double hi = Dist (p1, p2);
|
||||||
if (hi > maxh) maxh = hi;
|
if (hi > maxh) maxh = hi;
|
||||||
|
|
||||||
|
bbox.Add (p1);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
Point3d mpmin = bbox.PMin();
|
||||||
|
Point3d mpmax = bbox.PMax();
|
||||||
|
Point3d mpc = Center (mpmin, mpmax);
|
||||||
|
double d = max3(mpmax.X()-mpmin.X(),
|
||||||
|
mpmax.Y()-mpmin.Y(),
|
||||||
|
mpmax.Z()-mpmin.Z()) / 2;
|
||||||
|
mpmin = mpc - Vec3d (d, d, d);
|
||||||
|
mpmax = mpc + Vec3d (d, d, d);
|
||||||
|
Box3d meshbox (mpmin, mpmax);
|
||||||
|
|
||||||
|
LocalH loch2 (mpmin, mpmax, 1);
|
||||||
|
|
||||||
if (mp.maxh < maxh) maxh = mp.maxh;
|
if (mp.maxh < maxh) maxh = mp.maxh;
|
||||||
|
|
||||||
// auto loch_ptr = mesh.LocalHFunction().Copy();
|
auto loch_ptr = mesh.LocalHFunction().Copy(bbox);
|
||||||
// auto & loch = *loch_ptr;
|
auto & loch = *loch_ptr;
|
||||||
auto & loch = mesh.LocalHFunction();
|
|
||||||
|
|
||||||
bool changed;
|
bool changed;
|
||||||
static Timer t1("loop1");
|
static Timer t1("loop1");
|
||||||
@ -1174,60 +1191,6 @@ void Meshing3 :: PrepareBlockFillLocalH (Mesh & mesh,
|
|||||||
while (changed);
|
while (changed);
|
||||||
t1.Stop();
|
t1.Stop();
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
void Meshing3 :: BlockFillLocalH (Mesh & mesh,
|
|
||||||
const MeshingParameters & mp)
|
|
||||||
{
|
|
||||||
static Timer t("Mesing3::BlockFillLocalH"); RegionTimer reg(t);
|
|
||||||
|
|
||||||
if (!mesh.HasLocalHFunction())
|
|
||||||
{
|
|
||||||
mesh.CalcLocalH(mp.grading);
|
|
||||||
PrepareBlockFillLocalH(mesh, mp);
|
|
||||||
}
|
|
||||||
|
|
||||||
double filldist = mp.filldist;
|
|
||||||
|
|
||||||
// (*testout) << "blockfill local h" << endl;
|
|
||||||
// (*testout) << "rel filldist = " << filldist << endl;
|
|
||||||
PrintMessage (3, "blockfill local h");
|
|
||||||
|
|
||||||
Box<3> bbox ( Box<3>::EMPTY_BOX );
|
|
||||||
double maxh = 0;
|
|
||||||
|
|
||||||
for (int i = 1; i <= adfront->GetNF(); i++)
|
|
||||||
{
|
|
||||||
const MiniElement2d & el = adfront->GetFace(i);
|
|
||||||
for (int j = 1; j <= 3; j++)
|
|
||||||
{
|
|
||||||
const Point3d & p1 = adfront->GetPoint (el.PNumMod(j));
|
|
||||||
const Point3d & p2 = adfront->GetPoint (el.PNumMod(j+1));
|
|
||||||
|
|
||||||
double hi = Dist (p1, p2);
|
|
||||||
if (hi > maxh) maxh = hi;
|
|
||||||
|
|
||||||
bbox.Add (p1);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
Point3d mpmin = bbox.PMin();
|
|
||||||
Point3d mpmax = bbox.PMax();
|
|
||||||
Point3d mpc = Center (mpmin, mpmax);
|
|
||||||
double d = max3(mpmax.X()-mpmin.X(),
|
|
||||||
mpmax.Y()-mpmin.Y(),
|
|
||||||
mpmax.Z()-mpmin.Z()) / 2;
|
|
||||||
mpmin = mpc - Vec3d (d, d, d);
|
|
||||||
mpmax = mpc + Vec3d (d, d, d);
|
|
||||||
Box3d meshbox (mpmin, mpmax);
|
|
||||||
|
|
||||||
LocalH loch2 (mpmin, mpmax, 1);
|
|
||||||
|
|
||||||
if (mp.maxh < maxh) maxh = mp.maxh;
|
|
||||||
|
|
||||||
if (debugparam.slowchecks)
|
if (debugparam.slowchecks)
|
||||||
(*testout) << "Blockfill with points: " << endl;
|
(*testout) << "Blockfill with points: " << endl;
|
||||||
for (int i = 1; i <= npoints.Size(); i++)
|
for (int i = 1; i <= npoints.Size(); i++)
|
||||||
|
@ -28,7 +28,6 @@ class Meshing3
|
|||||||
NgArray<char*> problems;
|
NgArray<char*> problems;
|
||||||
/// tolerance criterion
|
/// tolerance criterion
|
||||||
double tolfak;
|
double tolfak;
|
||||||
NgArray<Point<3> > npoints;
|
|
||||||
public:
|
public:
|
||||||
///
|
///
|
||||||
Meshing3 (const string & rulefilename);
|
Meshing3 (const string & rulefilename);
|
||||||
@ -64,7 +63,6 @@ public:
|
|||||||
///
|
///
|
||||||
void BlockFill (Mesh & mesh, double gh);
|
void BlockFill (Mesh & mesh, double gh);
|
||||||
///
|
///
|
||||||
void PrepareBlockFillLocalH (Mesh & mesh, const MeshingParameters & mp);
|
|
||||||
void BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp);
|
void BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp);
|
||||||
|
|
||||||
/// uses points of adfront, and puts new elements into mesh
|
/// uses points of adfront, and puts new elements into mesh
|
||||||
|
@ -134,7 +134,6 @@ void ResetStatus()
|
|||||||
|
|
||||||
void PushStatus(const MyStr& s)
|
void PushStatus(const MyStr& s)
|
||||||
{
|
{
|
||||||
return;
|
|
||||||
msgstatus_stack.Append(new MyStr (s));
|
msgstatus_stack.Append(new MyStr (s));
|
||||||
SetStatMsg(s);
|
SetStatMsg(s);
|
||||||
threadpercent_stack.Append(0);
|
threadpercent_stack.Append(0);
|
||||||
@ -142,7 +141,6 @@ void PushStatus(const MyStr& s)
|
|||||||
|
|
||||||
void PushStatusF(const MyStr& s)
|
void PushStatusF(const MyStr& s)
|
||||||
{
|
{
|
||||||
return;
|
|
||||||
msgstatus_stack.Append(new MyStr (s));
|
msgstatus_stack.Append(new MyStr (s));
|
||||||
SetStatMsg(s);
|
SetStatMsg(s);
|
||||||
threadpercent_stack.Append(0);
|
threadpercent_stack.Append(0);
|
||||||
@ -151,7 +149,6 @@ void PushStatusF(const MyStr& s)
|
|||||||
|
|
||||||
void PopStatus()
|
void PopStatus()
|
||||||
{
|
{
|
||||||
return;
|
|
||||||
if (msgstatus_stack.Size())
|
if (msgstatus_stack.Size())
|
||||||
{
|
{
|
||||||
if (msgstatus_stack.Size() > 1)
|
if (msgstatus_stack.Size() > 1)
|
||||||
|
@ -1059,18 +1059,18 @@
|
|||||||
},
|
},
|
||||||
{
|
{
|
||||||
"angles_tet": [
|
"angles_tet": [
|
||||||
25.826,
|
25.773,
|
||||||
139.09
|
140.19
|
||||||
],
|
],
|
||||||
"angles_trig": [
|
"angles_trig": [
|
||||||
25.414,
|
25.416,
|
||||||
117.11
|
117.5
|
||||||
],
|
],
|
||||||
"ne1d": 0,
|
"ne1d": 0,
|
||||||
"ne2d": 1616,
|
"ne2d": 1616,
|
||||||
"ne3d": 5547,
|
"ne3d": 5533,
|
||||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 2, 8, 53, 131, 290, 453, 688, 906, 1032, 1025, 755, 204]",
|
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 1, 10, 48, 135, 266, 450, 664, 886, 1061, 1008, 779, 225]",
|
||||||
"total_badness": 7050.1376548
|
"total_badness": 7010.7623078
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
"angles_tet": [
|
"angles_tet": [
|
||||||
@ -3119,18 +3119,18 @@
|
|||||||
},
|
},
|
||||||
{
|
{
|
||||||
"angles_tet": [
|
"angles_tet": [
|
||||||
21.326,
|
21.239,
|
||||||
146.04
|
146.04
|
||||||
],
|
],
|
||||||
"angles_trig": [
|
"angles_trig": [
|
||||||
23.262,
|
23.611,
|
||||||
121.87
|
121.98
|
||||||
],
|
],
|
||||||
"ne1d": 0,
|
"ne1d": 0,
|
||||||
"ne2d": 5890,
|
"ne2d": 5890,
|
||||||
"ne3d": 25271,
|
"ne3d": 25307,
|
||||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 1, 10, 41, 108, 355, 777, 1656, 2789, 4110, 5146, 5263, 3910, 1105]",
|
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 1, 10, 42, 120, 371, 783, 1649, 2855, 4238, 5103, 5217, 3821, 1097]",
|
||||||
"total_badness": 31387.456922
|
"total_badness": 31484.35982
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
"angles_tet": [
|
"angles_tet": [
|
||||||
@ -3145,7 +3145,7 @@
|
|||||||
"ne2d": 16290,
|
"ne2d": 16290,
|
||||||
"ne3d": 174928,
|
"ne3d": 174928,
|
||||||
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 3, 71, 256, 868, 2763, 7163, 15703, 26959, 37225, 41474, 32186, 10257]",
|
"quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 3, 71, 256, 868, 2763, 7163, 15703, 26959, 37225, 41474, 32186, 10257]",
|
||||||
"total_badness": 211389.4278
|
"total_badness": 211389.42727
|
||||||
}
|
}
|
||||||
],
|
],
|
||||||
"trafo.geo": [
|
"trafo.geo": [
|
||||||
|
Loading…
Reference in New Issue
Block a user