ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-24 21:10:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

separate localh trees for different layers

Browse Source 
currenlty used in OCC geometries generated with
shape = netgen.occ.Compound(list_of_shapes, separate_layers=True)
This commit is contained in:
mhochsteger@cerbsim.com 2022-03-10 19:04:44 +01:00
parent eea8054af6
commit 154302605f
11 changed files with 133 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
libsrc/meshing/basegeom.cpp
View File

@ -454,7 +454,7 @@ namespace netgen
mesh.LoadLocalMeshSize(mparam.meshsizefilename);
}
void DivideEdge(GeometryEdge * edge, const MeshingParameters & mparam, const Mesh & mesh, Array<Point<3>> & points, Array<double> & params)
void DivideEdge(GeometryEdge * edge, const MeshingParameters & mparam, const Mesh & mesh, Array<Point<3>> & points, Array<double> & params, int layer)
{
static Timer tdivedgesections("Divide edge sections");
static Timer tdivide("Divide Edges");
@ -470,7 +470,7 @@ namespace netgen
for(auto i : Range(divide_edge_sections))
{
auto pt = edge->GetPoint(double(i+1)/divide_edge_sections);
hvalue[i+1] = hvalue[i] + 1./mesh.GetH(pt) * (pt-old_pt).Length();
hvalue[i+1] = hvalue[i] + 1./mesh.GetH(pt, layer) * (pt-old_pt).Length();
old_pt = pt;
}
int nsubedges = max2(1, int(floor(hvalue[divide_edge_sections]+0.5)));
@ -525,7 +525,7 @@ namespace netgen
std::map<size_t, PointIndex> vert2meshpt;
for(auto & vert : vertices)
{
auto pi = mesh.AddPoint(vert->GetPoint());
auto pi = mesh.AddPoint(vert->GetPoint(), vert->properties.layer);
tree.Insert(mesh[pi], pi);
vert2meshpt[vert->GetHash()] = pi;
mesh[pi].Singularity(vert->properties.hpref);
@ -591,7 +591,7 @@ namespace netgen
}
else
{
DivideEdge(edge, mparam, mesh, edge_points, params);
DivideEdge(edge, mparam, mesh, edge_points, params, edge->properties.layer);
}
}
else
@ -635,7 +635,7 @@ namespace netgen
for(auto i : Range(edge_points))
{
auto pi = mesh.AddPoint(edge_points[i]);
auto pi = mesh.AddPoint(edge_points[i], edge->properties.layer);
tree.Insert(mesh[pi], pi);
pnums[i+1] = pi;
}
@ -724,7 +724,7 @@ namespace netgen
static Timer t("GenerateMesh"); RegionTimer reg(t);
MESHING2_RESULT res = meshing.GenerateMesh(mesh, mparam, mparam.maxh, k+1);
MESHING2_RESULT res = meshing.GenerateMesh(mesh, mparam, mparam.maxh, k+1, face.properties.layer);
for(auto i : Range(noldsurfels, mesh.GetNSE()))
{

5
libsrc/meshing/basegeom.hpp
View File

@ -17,6 +17,7 @@ namespace netgen
optional<Vec<4>> col;
double maxh = 1e99;
double hpref = 0; // number of hp refinement levels (will be multiplied by factor later)
int layer = 1;
optional<bool> quad_dominated;
void Merge(const ShapeProperties & prop2)
{
@ -25,6 +26,7 @@ namespace netgen
maxh = min2(maxh, prop2.maxh);
hpref = max2(hpref, prop2.hpref);
if(!quad_dominated.has_value()) quad_dominated = prop2.quad_dominated;
layer = max(layer, prop2.layer);
}
string GetName() const { return name ? *name : "default"; }
@ -32,7 +34,7 @@ namespace netgen
void DoArchive(Archive& ar)
{
ar & name & col & maxh & hpref;
ar & name & col & maxh & hpref & layer;
}
};
@ -51,6 +53,7 @@ namespace netgen
{
public:
int nr = -1;
int layer = 1;
ShapeProperties properties;
Array<ShapeIdentification> identifications;
GeometryShape * primary;

68
libsrc/meshing/meshclass.cpp
View File

@ -229,7 +229,7 @@ namespace netgen
surfelementht = nullptr;
segmentht = nullptr;
lochfunc = nullptr;
lochfunc = {nullptr};
// mglevels = 1;
elementsearchtree = nullptr;
elementsearchtreets = NextTimeStamp();
@ -3264,7 +3264,7 @@ namespace netgen
void Mesh :: SetLocalH (netgen::Point<3> pmin, netgen::Point<3> pmax, double grading)
void Mesh :: SetLocalH (netgen::Point<3> pmin, netgen::Point<3> pmax, double grading, int layer)
{
using netgen::Point;
Point<3> c = Center (pmin, pmax);
@ -3275,30 +3275,30 @@ namespace netgen
Point<3> pmin2 = c - Vec<3> (d, d, d);
Point<3> pmax2 = c + Vec<3> (d, d, d);
lochfunc = make_unique<LocalH> (pmin2, pmax2, grading, dimension);
SetLocalH(make_unique<LocalH> (pmin2, pmax2, grading, dimension), layer);
}
void Mesh :: RestrictLocalH (const Point3d & p, double hloc)
void Mesh :: RestrictLocalH (const Point3d & p, double hloc, int layer)
{
if(hloc < hmin)
hloc = hmin;
//cout << "restrict h in " << p << " to " << hloc << endl;
if (!lochfunc)
if (!lochfunc[layer-1])
{
PrintWarning("RestrictLocalH called, creating mesh-size tree");
Point3d boxmin, boxmax;
GetBox (boxmin, boxmax);
SetLocalH (boxmin, boxmax, 0.8);
SetLocalH (boxmin, boxmax, 0.8, layer);
}
lochfunc -> SetH (p, hloc);
lochfunc[layer-1] -> SetH (p, hloc);
}
void Mesh :: RestrictLocalHLine (const Point3d & p1,
const Point3d & p2,
double hloc)
double hloc, int layer)
{
if(hloc < hmin)
hloc = hmin;
@ -3311,7 +3311,7 @@ namespace netgen
for (i = 0; i <= steps; i++)
{
Point3d p = p1 + (double(i)/double(steps) * v);
RestrictLocalH (p, hloc);
RestrictLocalH (p, hloc, layer);
}
}
@ -3343,24 +3343,24 @@ namespace netgen
}
double Mesh :: GetH (const Point3d & p) const
double Mesh :: GetH (const Point3d & p, int layer) const
{
double hmin = hglob;
if (lochfunc)
if (lochfunc[layer-1])
{
double hl = lochfunc->GetH (p);
double hl = lochfunc[layer-1]->GetH (p);
if (hl < hglob)
hmin = hl;
}
return hmin;
}
double Mesh :: GetMinH (const Point3d & pmin, const Point3d & pmax)
double Mesh :: GetMinH (const Point3d & pmin, const Point3d & pmax, int layer)
{
double hmin = hglob;
if (lochfunc)
if (lochfunc[layer-1])
{
double hl = lochfunc->GetMinH (pmin, pmax);
double hl = lochfunc[layer-1]->GetMinH (pmin, pmax);
if (hl < hmin)
hmin = hl;
}
@ -3404,16 +3404,16 @@ namespace netgen
void Mesh :: CalcLocalH (double grading)
void Mesh :: CalcLocalH (double grading, int layer)
{
static Timer t("Mesh::CalcLocalH"); RegionTimer reg(t);
if (!lochfunc)
if (!lochfunc[layer-1])
{
Point3d pmin, pmax;
GetBox (pmin, pmax);
// SetLocalH (pmin, pmax, mparam.grading);
SetLocalH (pmin, pmax, grading);
SetLocalH (pmin, pmax, grading, layer);
}
PrintMessage (3,
@ -3459,7 +3459,7 @@ namespace netgen
const Point3d & p1 = points[el.PNum(1)];
const Point3d & p2 = points[el.PNum(2)];
const Point3d & p3 = points[el.PNum(3)];
lochfunc->SetH (Center (p1, p2, p3), hel);
lochfunc[layer-1]->SetH (Center (p1, p2, p3), hel);
}
}
else
@ -3467,12 +3467,12 @@ namespace netgen
{
const Point3d & p1 = points[el.PNum(1)];
const Point3d & p2 = points[el.PNum(2)];
lochfunc->SetH (Center (p1, p2), 2 * Dist (p1, p2));
lochfunc[layer-1]->SetH (Center (p1, p2), 2 * Dist (p1, p2));
}
{
const Point3d & p1 = points[el.PNum(3)];
const Point3d & p2 = points[el.PNum(4)];
lochfunc->SetH (Center (p1, p2), 2 * Dist (p1, p2));
lochfunc[layer-1]->SetH (Center (p1, p2), 2 * Dist (p1, p2));
}
}
}
@ -3490,7 +3490,7 @@ namespace netgen
*/
if (!ident -> UsedSymmetric (seg[0], seg[1]))
{
lochfunc->SetH (Center (p1, p2), Dist (p1, p2));
lochfunc[layer-1]->SetH (Center (p1, p2), Dist (p1, p2));
}
}
/*
@ -3540,17 +3540,17 @@ namespace netgen
}
void Mesh :: CalcLocalHFromPointDistances(double grading)
void Mesh :: CalcLocalHFromPointDistances(double grading, int layer)
{
PrintMessage (3, "Calculating local h from point distances");
if (!lochfunc)
if (!lochfunc[layer-1])
{
Point3d pmin, pmax;
GetBox (pmin, pmax);
// SetLocalH (pmin, pmax, mparam.grading);
SetLocalH (pmin, pmax, grading);
SetLocalH (pmin, pmax, grading, layer);
}
PointIndex i,j;
@ -3575,17 +3575,17 @@ namespace netgen
}
void Mesh :: CalcLocalHFromSurfaceCurvature (double grading, double elperr)
void Mesh :: CalcLocalHFromSurfaceCurvature (double grading, double elperr, int layer)
{
PrintMessage (3, "Calculating local h from surface curvature");
if (!lochfunc)
if (!lochfunc[layer-1])
{
Point3d pmin, pmax;
GetBox (pmin, pmax);
// SetLocalH (pmin, pmax, mparam.grading);
SetLocalH (pmin, pmax, grading);
SetLocalH (pmin, pmax, grading, layer);
}
@ -3839,6 +3839,18 @@ namespace netgen
void Mesh :: SetLocalH(shared_ptr<LocalH> loch, int layer)
{
if(layer>lochfunc.Size())
{
auto pre_size = lochfunc.Size();
lochfunc.SetSize(layer);
for(auto & func : lochfunc.Range(pre_size, layer-1))
func = lochfunc[0];
}
lochfunc[layer-1] = loch;
}
void Mesh :: GetBox (Point3d & pmin, Point3d & pmax, int dom) const
{
if (points.Size() == 0)

31
libsrc/meshing/meshclass.hpp
View File

@ -69,9 +69,9 @@ namespace netgen
/**
Representation of local mesh-size h
Representation of local mesh-size h (one function per mesh layer)
*/
shared_ptr<LocalH> lochfunc;
Array<shared_ptr<LocalH>> lochfunc;
///
double hglob;
///
@ -435,18 +435,18 @@ namespace netgen
*/
DLL_HEADER double AverageH (int surfnr = 0) const;
/// Calculates localh
DLL_HEADER void CalcLocalH (double grading);
DLL_HEADER void CalcLocalH (double grading, int layer=1);
///
DLL_HEADER void SetLocalH (netgen::Point<3> pmin, netgen::Point<3> pmax, double grading);
DLL_HEADER void SetLocalH (netgen::Point<3> pmin, netgen::Point<3> pmax, double grading, int layer=1);
///
DLL_HEADER void RestrictLocalH (const Point3d & p, double hloc);
DLL_HEADER void RestrictLocalH (const Point3d & p, double hloc, int layer=1);
///
DLL_HEADER void RestrictLocalHLine (const Point3d & p1, const Point3d & p2,
double hloc);
double hloc, int layer=1);
/// number of elements per radius
DLL_HEADER void CalcLocalHFromSurfaceCurvature(double grading, double elperr);
DLL_HEADER void CalcLocalHFromSurfaceCurvature(double grading, double elperr, int layer=1);
///
DLL_HEADER void CalcLocalHFromPointDistances(double grading);
DLL_HEADER void CalcLocalHFromPointDistances(double grading, int layer=1);
///
DLL_HEADER void RestrictLocalH (resthtype rht, int nr, double loch);
///
@ -460,19 +460,20 @@ namespace netgen
///
DLL_HEADER void SetMaxHDomain (const NgArray<double> & mhd);
///
DLL_HEADER double GetH (const Point3d & p) const;
DLL_HEADER double GetH (const Point3d & p, int layer=1) const;
DLL_HEADER double GetH (PointIndex pi) const { return GetH(points[pi], points[pi].GetLayer()); }
///
double GetMinH (const Point3d & pmin, const Point3d & pmax);
double GetMinH (const Point3d & pmin, const Point3d & pmax, int layer=1);
///
bool HasLocalHFunction () { return lochfunc != nullptr; }
bool HasLocalHFunction (int layer=1) { return lochfunc[layer-1] != nullptr; }
///
LocalH & LocalHFunction () { return * lochfunc; }
LocalH & LocalHFunction (int layer=1) { return * lochfunc[layer-1]; }
shared_ptr<LocalH> GetLocalH() const { return lochfunc; }
void SetLocalH(shared_ptr<LocalH> loch) { lochfunc = loch; }
shared_ptr<LocalH> GetLocalH(int layer=1) const { return lochfunc[layer-1]; }
void SetLocalH(shared_ptr<LocalH> loch, int layer=1);
///
bool LocalHFunctionGenerated(void) const { return (lochfunc != NULL); }
bool LocalHFunctionGenerated(int layer=1) const { return (lochfunc[layer-1] != NULL); }
/// Find bounding box
DLL_HEADER void GetBox (Point3d & pmin, Point3d & pmax, int dom = -1) const;

4
libsrc/meshing/meshing2.cpp
View File

@ -241,7 +241,7 @@ namespace netgen
MESHING2_RESULT Meshing2 :: GenerateMesh (Mesh & mesh, const MeshingParameters & mp, double gh, int facenr)
MESHING2_RESULT Meshing2 :: GenerateMesh (Mesh & mesh, const MeshingParameters & mp, double gh, int facenr, int layer)
{
static Timer timer("surface meshing"); RegionTimer reg(timer);
@ -477,7 +477,7 @@ namespace netgen
double his = Dist (p1, p2);
Point<3> pmid = Center (p1, p2);
double hshould = CalcLocalH (pmid, mesh.GetH (pmid));
double hshould = CalcLocalH (pmid, mesh.GetH (pmid, layer));
if (gh < hshould) hshould = gh;
mesh.RestrictLocalH (pmid, hshould);

2
libsrc/meshing/meshing2.hpp
View File

@ -59,7 +59,7 @@ public:
void LoadRules (const char * filename, bool quad);
///
DLL_HEADER MESHING2_RESULT GenerateMesh (Mesh & mesh, const MeshingParameters & mp, double gh, int facenr);
DLL_HEADER MESHING2_RESULT GenerateMesh (Mesh & mesh, const MeshingParameters & mp, double gh, int facenr, int layer=1);
DLL_HEADER void Delaunay (Mesh & mesh, int domainnr, const MeshingParameters & mp);
DLL_HEADER void BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp);

1
libsrc/meshing/meshtype.hpp
View File

@ -1350,6 +1350,7 @@ namespace netgen
public:
Point<3> pnt;
double h;
int layer = 1;
MeshSizePoint (Point<3> _pnt, double _h) : pnt(_pnt), h(_h) { ; }
MeshSizePoint () = default;
MeshSizePoint (const MeshSizePoint &) = default;

20
libsrc/meshing/smoothing3.cpp
View File

@ -1243,8 +1243,8 @@ void Mesh :: ImproveMesh (const CSG eometry & geometry, OPTIMIZEGOAL goal)
// if (uselocalh)
{
double lh = GetH (points.Get(i));
pf->SetLocalH (GetH (points.Get(i)));
double lh = GetH(points.Get(i), points.Get(i).GetLayer());
pf->SetLocalH (lh);
par.typx = lh / 10;
// (*testout) << "lh(" << points.Get(i) << ") = " << lh << "\n";
}
@ -1366,10 +1366,10 @@ void Mesh :: ImproveMeshSequential (const MeshingParameters & mp, OPTIMIZEGOAL g
NgArray<double, PointIndex::BASE> pointh (points.Size());
if(lochfunc)
if(HasLocalHFunction())
{
for (PointIndex pi : points.Range())
pointh[pi] = GetH(points[pi]);
pointh[pi] = GetH(pi);
}
else
{
@ -1505,13 +1505,13 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
Array<double, PointIndex> pointh (points.Size());
if(lochfunc)
if(HasLocalHFunction())
{
RegionTimer rt(tloch);
ParallelForRange(points.Range(), [&] (auto myrange)
{
for(auto pi : myrange)
pointh[pi] = GetH(points[pi]);
pointh[pi] = GetH(pi);
});
}
else
@ -1640,11 +1640,11 @@ void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp,
NgArray<double, PointIndex::BASE, PointIndex> pointh (points.Size());
if(lochfunc)
if(HasLocalHFunction())
{
// for(i = 1; i<=points.Size(); i++)
for (PointIndex pi : points.Range())
pointh[pi] = GetH(points[pi]);
pointh[pi] = GetH(pi);
}
else
{
@ -1797,11 +1797,11 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp,
NgArray<double, PointIndex::BASE> pointh (points.Size());
if(lochfunc)
if(HasLocalHFunction())
{
// for(i=1; i<=points.Size(); i++)
for (PointIndex pi : points.Range())
pointh[pi] = GetH(points[pi]);
pointh[pi] = GetH(pi);
}
else
{

56
libsrc/occ/occgenmesh.cpp
View File

@ -78,8 +78,7 @@ namespace netgen
void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
BRepLProp_SLProps * prop, BRepLProp_SLProps * prop2, Mesh & mesh, int depth, double h,
const MeshingParameters & mparam)
BRepLProp_SLProps * prop, BRepLProp_SLProps * prop2, Mesh & mesh, int depth, double h, int layer, const MeshingParameters & mparam)
{
int ls = -1;
@ -190,20 +189,20 @@ namespace netgen
if(ls == 0)
{
pm.SetX(0.5*(par1.X()+par2.X())); pm.SetY(0.5*(par1.Y()+par2.Y()));
RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, layer, mparam);
RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, layer, mparam);
}
else if(ls == 1)
{
pm.SetX(0.5*(par0.X()+par2.X())); pm.SetY(0.5*(par0.Y()+par2.Y()));
RestrictHTriangle(pm, par1, par2, prop, prop2, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par1, par2, prop, prop2, mesh, depth+1, h, layer, mparam);
RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, layer, mparam);
}
else if(ls == 2)
{
pm.SetX(0.5*(par0.X()+par1.X())); pm.SetY(0.5*(par0.Y()+par1.Y()));
RestrictHTriangle(pm, par1, par2, prop, prop2, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, mparam);
RestrictHTriangle(pm, par1, par2, prop, prop2, mesh, depth+1, h, layer, mparam);
RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, layer, mparam);
}
}
@ -215,16 +214,16 @@ namespace netgen
prop->SetParameters (parmid.X(), parmid.Y());
pnt = prop->Value();
p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
mesh.RestrictLocalH (p3d, h);
mesh.RestrictLocalH (p3d, h, layer);
p3d = Point3d(pnt0.X(), pnt0.Y(), pnt0.Z());
mesh.RestrictLocalH (p3d, h);
mesh.RestrictLocalH (p3d, h, layer);
p3d = Point3d(pnt1.X(), pnt1.Y(), pnt1.Z());
mesh.RestrictLocalH (p3d, h);
mesh.RestrictLocalH (p3d, h, layer);
p3d = Point3d(pnt2.X(), pnt2.Y(), pnt2.Z());
mesh.RestrictLocalH (p3d, h);
mesh.RestrictLocalH (p3d, h, layer);
//(*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl;
@ -375,6 +374,7 @@ namespace netgen
//double maxh = mparam.maxh;
// int noldpoints = mesh->GetNP();
int noldsurfel = mesh.GetNSE();
int layer = OCCGeometry::global_shape_properties[TopoDS::Face(geom.fmap(k)).TShape()].layer;
static Timer tsurfprop("surfprop");
tsurfprop.Start();
@ -391,7 +391,7 @@ namespace netgen
try {
static Timer t("GenerateMesh"); RegionTimer reg(t);
res = meshing.GenerateMesh (mesh, mparam_without_overlap, maxh, k);
res = meshing.GenerateMesh (mesh, mparam_without_overlap, maxh, k, layer);
}
catch (SingularMatrixException)
@ -437,24 +437,28 @@ namespace netgen
NgArray<double> maxhdom;
maxhdom.SetSize (geom.NrSolids());
maxhdom = mparam.maxh;
int maxlayer = 1;
int dom = 0;
for (TopExp_Explorer e(geom.GetShape(), TopAbs_SOLID); e.More(); e.Next(), dom++)
{
maxhdom[dom] = min2(maxhdom[dom], OCCGeometry::global_shape_properties[e.Current().TShape()].maxh);
maxlayer = max2(maxlayer, OCCGeometry::global_shape_properties[e.Current().TShape()].layer);
}
mesh.SetMaxHDomain (maxhdom);
Box<3> bb = geom.GetBoundingBox();
bb.Increase (bb.Diam()/10);
mesh.SetLocalH (bb.PMin(), bb.PMax(), 0.5);
if (mparam.uselocalh)
{
const char * savetask = multithread.task;
multithread.percent = 0;
mesh.SetLocalH (bb.PMin(), bb.PMax(), mparam.grading);
for(auto layer : Range(1, maxlayer+1))
mesh.SetLocalH (bb.PMin(), bb.PMax(), mparam.grading, layer);
int nedges = geom.emap.Extent();
@ -482,6 +486,7 @@ namespace netgen
for (int i = 1; i <= nedges && !multithread.terminate; i++)
{
TopoDS_Edge e = TopoDS::Edge (geom.emap(i));
int layer = OCCGeometry::global_shape_properties[e.TShape()].layer;
multithread.percent = 100 * (i-1)/double(nedges);
if (BRep_Tool::Degenerated(e)) continue;
@ -540,7 +545,7 @@ namespace netgen
for (int j = 0; j <= maxj; j++)
{
gp_Pnt pnt = c->Value (s0+double(j)/maxj*(s1-s0));
mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh);
mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh, layer);
}
}
@ -555,6 +560,7 @@ namespace netgen
double maxcur = 0;
multithread.percent = 100 * (i-1)/double(nedges);
TopoDS_Edge edge = TopoDS::Edge (geom.emap(i));
int layer = OCCGeometry::global_shape_properties[edge.TShape()].layer;
if (BRep_Tool::Degenerated(edge)) continue;
double s0, s1;
Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1);
@ -575,7 +581,7 @@ namespace netgen
gp_Pnt pnt = c->Value (s);
mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), ComputeH (fabs(curvature), mparam));
mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), ComputeH (fabs(curvature), mparam), layer);
}
}
@ -589,6 +595,7 @@ namespace netgen
{
multithread.percent = 100 * (i-1)/double(nfaces);
TopoDS_Face face = TopoDS::Face(geom.fmap(i));
int layer = OCCGeometry::global_shape_properties[face.TShape()].layer;
TopLoc_Location loc;
Handle(Geom_Surface) surf = BRep_Tool::Surface (face);
Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc);
@ -628,7 +635,7 @@ namespace netgen
//maxside = max (maxside, p[1].Distance(p[2]));
//cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush;
RestrictHTriangle (par[0], par[1], par[2], &prop, &prop2, mesh, 0, 0, mparam);
RestrictHTriangle (par[0], par[1], par[2], &prop, &prop2, mesh, 0, 0, layer, mparam);
//cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush;
}
}
@ -654,6 +661,7 @@ namespace netgen
for (int i = 1; i <= nedges && !multithread.terminate; i++)
{
TopoDS_Edge edge = TopoDS::Edge (geom.emap(i));
int layer = OCCGeometry::global_shape_properties[edge.TShape()].layer;
if (BRep_Tool::Degenerated(edge)) continue;
double s0, s1;
@ -678,6 +686,7 @@ namespace netgen
gp_Pnt p1 = c->Value (s);
lines[nlines].p0 = Point<3> (p0.X(), p0.Y(), p0.Z());
lines[nlines].p1 = Point<3> (p1.X(), p1.Y(), p1.Z());
lines[nlines].layer = layer;
Box3d box;
box.SetPoint (Point3d(lines[nlines].p0));
@ -719,8 +728,8 @@ namespace netgen
Box3d box;
box.SetPoint (Point3d(line.p0));
box.AddPoint (Point3d(line.p1));
double maxhline = max (mesh.GetH(box.PMin()),
mesh.GetH(box.PMax()));
double maxhline = max (mesh.GetH(box.PMin(), line.layer),
mesh.GetH(box.PMax(), line.layer));
box.Increase(maxhline);
double mindist = 1e99;
@ -731,6 +740,7 @@ namespace netgen
{
int num = linenums[j]-1;
if (i == num) continue;
if (line.layer != lines[num].layer) continue;
if( is_identified_edge(edgenumber[i], edgenumber[num]) ) continue;
if ((line.p0-lines[num].p0).Length2() < 1e-15) continue;
if ((line.p0-lines[num].p1).Length2() < 1e-15) continue;
@ -749,12 +759,12 @@ namespace netgen
mindist = 1e-3 * bb.Diam();
}
mesh.RestrictLocalHLine(line.p0, line.p1, mindist);
mesh.RestrictLocalHLine(line.p0, line.p1, mindist, line.layer);
}
}
for (auto mspnt : mparam.meshsize_points)
mesh.RestrictLocalH(mspnt.pnt, mspnt.h);
mesh.RestrictLocalH(mspnt.pnt, mspnt.h, mspnt.layer);
multithread.task = savetask;

1
libsrc/occ/occgeom.hpp
View File

@ -91,6 +91,7 @@ namespace netgen
{
public:
Point<3> p0, p1;
int layer = 1;
double Dist (Line l);
double Length () { return (p1-p0).Length(); }
};

23
libsrc/occ/python_occ_shapes.cpp
View File

@ -1456,14 +1456,29 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
py::class_<TopoDS_Solid, TopoDS_Shape> (m, "Solid");
py::class_<TopoDS_Compound, TopoDS_Shape> (m, "Compound")
.def(py::init([](std::vector<TopoDS_Shape> shapes) {
.def(py::init([](std::vector<TopoDS_Shape> shapes, bool separate_layers) {
BRep_Builder builder;
TopoDS_Compound comp;
builder.MakeCompound(comp);
for(auto& s : shapes)
builder.Add(comp, s);
for(auto i : Range(shapes.size()))
{
builder.Add(comp, shapes[i]);
if(separate_layers)
{
auto & props = OCCGeometry::global_shape_properties;
for(auto & s : GetSolids(shapes[i]))
props[s.TShape()].layer = i+1;
for(auto & s : GetFaces(shapes[i]))
props[s.TShape()].layer = i+1;
for(auto & s : GetEdges(shapes[i]))
props[s.TShape()].layer = i+1;
for(auto & s : GetVertices(shapes[i]))
props[s.TShape()].layer = i+1;
}
}
return comp;
}))
}), py::arg("shapes"), py::arg("separate_layers")=false)
;