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

Merge branch 'separate_localh' into 'master'

Browse Source 
separate localh trees for different layers

See merge request jschoeberl/netgen!490
This commit is contained in:
Joachim Schöberl 2022-03-14 13:48:47 +00:00
commit 7f8e0f1276
11 changed files with 140 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); 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 tdivedgesections("Divide edge sections");
static Timer tdivide("Divide Edges"); static Timer tdivide("Divide Edges");
@ -470,7 +470,7 @@ namespace netgen
for(auto i : Range(divide_edge_sections)) for(auto i : Range(divide_edge_sections))
{ {
auto pt = edge->GetPoint(double(i+1)/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; old_pt = pt;
} }
int nsubedges = max2(1, int(floor(hvalue[divide_edge_sections]+0.5))); int nsubedges = max2(1, int(floor(hvalue[divide_edge_sections]+0.5)));
@ -525,7 +525,7 @@ namespace netgen
std::map<size_t, PointIndex> vert2meshpt; std::map<size_t, PointIndex> vert2meshpt;
for(auto & vert : vertices) for(auto & vert : vertices)
{ {
auto pi = mesh.AddPoint(vert->GetPoint()); auto pi = mesh.AddPoint(vert->GetPoint(), vert->properties.layer);
tree.Insert(mesh[pi], pi); tree.Insert(mesh[pi], pi);
vert2meshpt[vert->GetHash()] = pi; vert2meshpt[vert->GetHash()] = pi;
mesh[pi].Singularity(vert->properties.hpref); mesh[pi].Singularity(vert->properties.hpref);
@ -591,7 +591,7 @@ namespace netgen
} }
else else
{ {
DivideEdge(edge, mparam, mesh, edge_points, params); DivideEdge(edge, mparam, mesh, edge_points, params, edge->properties.layer);
} }
} }
else else
@ -635,7 +635,7 @@ namespace netgen
for(auto i : Range(edge_points)) 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); tree.Insert(mesh[pi], pi);
pnums[i+1] = pi; pnums[i+1] = pi;
} }
@ -724,7 +724,7 @@ namespace netgen
static Timer t("GenerateMesh"); RegionTimer reg(t); 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())) for(auto i : Range(noldsurfels, mesh.GetNSE()))
{ {

5
libsrc/meshing/basegeom.hpp
View File

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

70
libsrc/meshing/meshclass.cpp
View File

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

36
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; double hglob;
/// ///
@ -435,18 +435,18 @@ namespace netgen
*/ */
DLL_HEADER double AverageH (int surfnr = 0) const; DLL_HEADER double AverageH (int surfnr = 0) const;
/// Calculates localh /// 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, DLL_HEADER void RestrictLocalHLine (const Point3d & p1, const Point3d & p2,
double hloc); double hloc, int layer=1);
/// number of elements per radius /// 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); DLL_HEADER void RestrictLocalH (resthtype rht, int nr, double loch);
/// ///
@ -460,19 +460,25 @@ namespace netgen
/// ///
DLL_HEADER void SetMaxHDomain (const NgArray<double> & mhd); 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; } shared_ptr<LocalH> GetLocalH(int layer=1) const
void SetLocalH(shared_ptr<LocalH> loch) { lochfunc = loch; } {
if(lochfunc.Size() == 1)
return lochfunc[0];
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 /// Find bounding box
DLL_HEADER void GetBox (Point3d & pmin, Point3d & pmax, int dom = -1) const; 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); static Timer timer("surface meshing"); RegionTimer reg(timer);
@ -477,7 +477,7 @@ namespace netgen
double his = Dist (p1, p2); double his = Dist (p1, p2);
Point<3> pmid = Center (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; if (gh < hshould) hshould = gh;
mesh.RestrictLocalH (pmid, hshould); mesh.RestrictLocalH (pmid, hshould);

2
libsrc/meshing/meshing2.hpp
View File

@ -59,7 +59,7 @@ public:
void LoadRules (const char * filename, bool quad); 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 Delaunay (Mesh & mesh, int domainnr, const MeshingParameters & mp);
DLL_HEADER void BlockFillLocalH (Mesh & mesh, 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: public:
Point<3> pnt; Point<3> pnt;
double h; double h;
int layer = 1;
MeshSizePoint (Point<3> _pnt, double _h) : pnt(_pnt), h(_h) { ; } MeshSizePoint (Point<3> _pnt, double _h) : pnt(_pnt), h(_h) { ; }
MeshSizePoint () = default; MeshSizePoint () = default;
MeshSizePoint (const 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) // if (uselocalh)
{ {
double lh = GetH (points.Get(i)); double lh = GetH(points.Get(i), points.Get(i).GetLayer());
pf->SetLocalH (GetH (points.Get(i))); pf->SetLocalH (lh);
par.typx = lh / 10; par.typx = lh / 10;
// (*testout) << "lh(" << points.Get(i) << ") = " << lh << "\n"; // (*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()); NgArray<double, PointIndex::BASE> pointh (points.Size());
if(lochfunc) if(HasLocalHFunction())
{ {
for (PointIndex pi : points.Range()) for (PointIndex pi : points.Range())
pointh[pi] = GetH(points[pi]); pointh[pi] = GetH(pi);
} }
else else
{ {
@ -1505,13 +1505,13 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
Array<double, PointIndex> pointh (points.Size()); Array<double, PointIndex> pointh (points.Size());
if(lochfunc) if(HasLocalHFunction())
{ {
RegionTimer rt(tloch); RegionTimer rt(tloch);
ParallelForRange(points.Range(), [&] (auto myrange) ParallelForRange(points.Range(), [&] (auto myrange)
{ {
for(auto pi : myrange) for(auto pi : myrange)
pointh[pi] = GetH(points[pi]); pointh[pi] = GetH(pi);
}); });
} }
else else
@ -1640,11 +1640,11 @@ void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp,
NgArray<double, PointIndex::BASE, PointIndex> pointh (points.Size()); NgArray<double, PointIndex::BASE, PointIndex> pointh (points.Size());
if(lochfunc) if(HasLocalHFunction())
{ {
// for(i = 1; i<=points.Size(); i++) // for(i = 1; i<=points.Size(); i++)
for (PointIndex pi : points.Range()) for (PointIndex pi : points.Range())
pointh[pi] = GetH(points[pi]); pointh[pi] = GetH(pi);
} }
else else
{ {
@ -1797,11 +1797,11 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp,
NgArray<double, PointIndex::BASE> pointh (points.Size()); NgArray<double, PointIndex::BASE> pointh (points.Size());
if(lochfunc) if(HasLocalHFunction())
{ {
// for(i=1; i<=points.Size(); i++) // for(i=1; i<=points.Size(); i++)
for (PointIndex pi : points.Range()) for (PointIndex pi : points.Range())
pointh[pi] = GetH(points[pi]); pointh[pi] = GetH(pi);
} }
else 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, void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
BRepLProp_SLProps * prop, BRepLProp_SLProps * prop2, Mesh & mesh, int depth, double h, BRepLProp_SLProps * prop, BRepLProp_SLProps * prop2, Mesh & mesh, int depth, double h, int layer, const MeshingParameters & mparam)
const MeshingParameters & mparam)
{ {
int ls = -1; int ls = -1;
@ -190,20 +189,20 @@ namespace netgen
if(ls == 0) if(ls == 0)
{ {
pm.SetX(0.5*(par1.X()+par2.X())); pm.SetY(0.5*(par1.Y()+par2.Y())); 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, par2, par0, prop, prop2, mesh, depth+1, h, layer, mparam);
RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, mparam); RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, layer, mparam);
} }
else if(ls == 1) else if(ls == 1)
{ {
pm.SetX(0.5*(par0.X()+par2.X())); pm.SetY(0.5*(par0.Y()+par2.Y())); 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, par1, par2, prop, prop2, mesh, depth+1, h, layer, mparam);
RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, mparam); RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, layer, mparam);
} }
else if(ls == 2) else if(ls == 2)
{ {
pm.SetX(0.5*(par0.X()+par1.X())); pm.SetY(0.5*(par0.Y()+par1.Y())); 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, par1, par2, prop, prop2, mesh, depth+1, h, layer, mparam);
RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, 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()); prop->SetParameters (parmid.X(), parmid.Y());
pnt = prop->Value(); pnt = prop->Value();
p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z()); 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()); 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()); 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()); 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; //(*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl;
@ -375,6 +374,7 @@ namespace netgen
//double maxh = mparam.maxh; //double maxh = mparam.maxh;
// int noldpoints = mesh->GetNP(); // int noldpoints = mesh->GetNP();
int noldsurfel = mesh.GetNSE(); int noldsurfel = mesh.GetNSE();
int layer = OCCGeometry::global_shape_properties[TopoDS::Face(geom.fmap(k)).TShape()].layer;
static Timer tsurfprop("surfprop"); static Timer tsurfprop("surfprop");
tsurfprop.Start(); tsurfprop.Start();
@ -391,7 +391,7 @@ namespace netgen
try { try {
static Timer t("GenerateMesh"); RegionTimer reg(t); 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) catch (SingularMatrixException)
@ -437,24 +437,28 @@ namespace netgen
NgArray<double> maxhdom; NgArray<double> maxhdom;
maxhdom.SetSize (geom.NrSolids()); maxhdom.SetSize (geom.NrSolids());
maxhdom = mparam.maxh; maxhdom = mparam.maxh;
int maxlayer = 1;
int dom = 0; int dom = 0;
for (TopExp_Explorer e(geom.GetShape(), TopAbs_SOLID); e.More(); e.Next(), dom++) 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); 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); mesh.SetMaxHDomain (maxhdom);
Box<3> bb = geom.GetBoundingBox(); Box<3> bb = geom.GetBoundingBox();
bb.Increase (bb.Diam()/10); bb.Increase (bb.Diam()/10);
mesh.SetLocalH (bb.PMin(), bb.PMax(), 0.5);
if (mparam.uselocalh) if (mparam.uselocalh)
{ {
const char * savetask = multithread.task; const char * savetask = multithread.task;
multithread.percent = 0; 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(); int nedges = geom.emap.Extent();
@ -482,6 +486,7 @@ namespace netgen
for (int i = 1; i <= nedges && !multithread.terminate; i++) for (int i = 1; i <= nedges && !multithread.terminate; i++)
{ {
TopoDS_Edge e = TopoDS::Edge (geom.emap(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); multithread.percent = 100 * (i-1)/double(nedges);
if (BRep_Tool::Degenerated(e)) continue; if (BRep_Tool::Degenerated(e)) continue;
@ -540,7 +545,7 @@ namespace netgen
for (int j = 0; j <= maxj; j++) for (int j = 0; j <= maxj; j++)
{ {
gp_Pnt pnt = c->Value (s0+double(j)/maxj*(s1-s0)); 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; double maxcur = 0;
multithread.percent = 100 * (i-1)/double(nedges); multithread.percent = 100 * (i-1)/double(nedges);
TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); TopoDS_Edge edge = TopoDS::Edge (geom.emap(i));
int layer = OCCGeometry::global_shape_properties[edge.TShape()].layer;
if (BRep_Tool::Degenerated(edge)) continue; if (BRep_Tool::Degenerated(edge)) continue;
double s0, s1; double s0, s1;
Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1);
@ -575,7 +581,7 @@ namespace netgen
gp_Pnt pnt = c->Value (s); 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); multithread.percent = 100 * (i-1)/double(nfaces);
TopoDS_Face face = TopoDS::Face(geom.fmap(i)); TopoDS_Face face = TopoDS::Face(geom.fmap(i));
int layer = OCCGeometry::global_shape_properties[face.TShape()].layer;
TopLoc_Location loc; TopLoc_Location loc;
Handle(Geom_Surface) surf = BRep_Tool::Surface (face); Handle(Geom_Surface) surf = BRep_Tool::Surface (face);
Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc); Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc);
@ -628,7 +635,7 @@ namespace netgen
//maxside = max (maxside, p[1].Distance(p[2])); //maxside = max (maxside, p[1].Distance(p[2]));
//cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush; //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; //cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush;
} }
} }
@ -654,6 +661,7 @@ namespace netgen
for (int i = 1; i <= nedges && !multithread.terminate; i++) for (int i = 1; i <= nedges && !multithread.terminate; i++)
{ {
TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); TopoDS_Edge edge = TopoDS::Edge (geom.emap(i));
int layer = OCCGeometry::global_shape_properties[edge.TShape()].layer;
if (BRep_Tool::Degenerated(edge)) continue; if (BRep_Tool::Degenerated(edge)) continue;
double s0, s1; double s0, s1;
@ -678,6 +686,7 @@ namespace netgen
gp_Pnt p1 = c->Value (s); gp_Pnt p1 = c->Value (s);
lines[nlines].p0 = Point<3> (p0.X(), p0.Y(), p0.Z()); lines[nlines].p0 = Point<3> (p0.X(), p0.Y(), p0.Z());
lines[nlines].p1 = Point<3> (p1.X(), p1.Y(), p1.Z()); lines[nlines].p1 = Point<3> (p1.X(), p1.Y(), p1.Z());
lines[nlines].layer = layer;
Box3d box; Box3d box;
box.SetPoint (Point3d(lines[nlines].p0)); box.SetPoint (Point3d(lines[nlines].p0));
@ -719,8 +728,8 @@ namespace netgen
Box3d box; Box3d box;
box.SetPoint (Point3d(line.p0)); box.SetPoint (Point3d(line.p0));
box.AddPoint (Point3d(line.p1)); box.AddPoint (Point3d(line.p1));
double maxhline = max (mesh.GetH(box.PMin()), double maxhline = max (mesh.GetH(box.PMin(), line.layer),
mesh.GetH(box.PMax())); mesh.GetH(box.PMax(), line.layer));
box.Increase(maxhline); box.Increase(maxhline);
double mindist = 1e99; double mindist = 1e99;
@ -731,6 +740,7 @@ namespace netgen
{ {
int num = linenums[j]-1; int num = linenums[j]-1;
if (i == num) continue; if (i == num) continue;
if (line.layer != lines[num].layer) continue;
if( is_identified_edge(edgenumber[i], edgenumber[num]) ) 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].p0).Length2() < 1e-15) continue;
if ((line.p0-lines[num].p1).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(); 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) for (auto mspnt : mparam.meshsize_points)
mesh.RestrictLocalH(mspnt.pnt, mspnt.h); mesh.RestrictLocalH(mspnt.pnt, mspnt.h, mspnt.layer);
multithread.task = savetask; multithread.task = savetask;

1
libsrc/occ/occgeom.hpp
View File

@ -91,6 +91,7 @@ namespace netgen
{ {
public: public:
Point<3> p0, p1; Point<3> p0, p1;
int layer = 1;
double Dist (Line l); double Dist (Line l);
double Length () { return (p1-p0).Length(); } 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_Solid, TopoDS_Shape> (m, "Solid");
py::class_<TopoDS_Compound, TopoDS_Shape> (m, "Compound") 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; BRep_Builder builder;
TopoDS_Compound comp; TopoDS_Compound comp;
builder.MakeCompound(comp); builder.MakeCompound(comp);
for(auto& s : shapes) for(auto i : Range(shapes.size()))
builder.Add(comp, s); {
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; return comp;
})) }), py::arg("shapes"), py::arg("separate_layers")=false)
; ;