From 154302605f5ae6c10f584d1dd1ea34a1622566d3 Mon Sep 17 00:00:00 2001 From: "mhochsteger@cerbsim.com" Date: Thu, 10 Mar 2022 19:04:44 +0100 Subject: [PATCH 1/2] separate localh trees for different layers currenlty used in OCC geometries generated with shape = netgen.occ.Compound(list_of_shapes, separate_layers=True) --- libsrc/meshing/basegeom.cpp | 12 +++--- libsrc/meshing/basegeom.hpp | 5 ++- libsrc/meshing/meshclass.cpp | 68 +++++++++++++++++++------------- libsrc/meshing/meshclass.hpp | 31 ++++++++------- libsrc/meshing/meshing2.cpp | 4 +- libsrc/meshing/meshing2.hpp | 2 +- libsrc/meshing/meshtype.hpp | 1 + libsrc/meshing/smoothing3.cpp | 20 +++++----- libsrc/occ/occgenmesh.cpp | 56 +++++++++++++++----------- libsrc/occ/occgeom.hpp | 1 + libsrc/occ/python_occ_shapes.cpp | 23 +++++++++-- 11 files changed, 133 insertions(+), 90 deletions(-) diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp index ff994594..63b47bc7 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -454,7 +454,7 @@ namespace netgen mesh.LoadLocalMeshSize(mparam.meshsizefilename); } - void DivideEdge(GeometryEdge * edge, const MeshingParameters & mparam, const Mesh & mesh, Array> & points, Array & params) + void DivideEdge(GeometryEdge * edge, const MeshingParameters & mparam, const Mesh & mesh, Array> & points, Array & 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 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())) { diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index 31eeabba..e0bf18b6 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -17,6 +17,7 @@ namespace netgen optional> col; double maxh = 1e99; double hpref = 0; // number of hp refinement levels (will be multiplied by factor later) + int layer = 1; optional 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 identifications; GeometryShape * primary; diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 729554a5..21822503 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -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 (pmin2, pmax2, grading, dimension); + SetLocalH(make_unique (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 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) diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index 9346ffba..e0b10081 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -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 lochfunc; + Array> 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 & 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 GetLocalH() const { return lochfunc; } - void SetLocalH(shared_ptr loch) { lochfunc = loch; } + shared_ptr GetLocalH(int layer=1) const { return lochfunc[layer-1]; } + void SetLocalH(shared_ptr 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; diff --git a/libsrc/meshing/meshing2.cpp b/libsrc/meshing/meshing2.cpp index c7d60e8d..6700b70c 100644 --- a/libsrc/meshing/meshing2.cpp +++ b/libsrc/meshing/meshing2.cpp @@ -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); diff --git a/libsrc/meshing/meshing2.hpp b/libsrc/meshing/meshing2.hpp index 8eb233a4..a05cd3dd 100644 --- a/libsrc/meshing/meshing2.hpp +++ b/libsrc/meshing/meshing2.hpp @@ -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); diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 9a1ac204..0fb605d0 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -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; diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp index 165f37f1..b30b8156 100644 --- a/libsrc/meshing/smoothing3.cpp +++ b/libsrc/meshing/smoothing3.cpp @@ -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 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 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 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 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 { diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp index 14da5c9a..d36b5d88 100644 --- a/libsrc/occ/occgenmesh.cpp +++ b/libsrc/occ/occgenmesh.cpp @@ -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 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; diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index c36342aa..d266c69d 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -91,6 +91,7 @@ namespace netgen { public: Point<3> p0, p1; + int layer = 1; double Dist (Line l); double Length () { return (p1-p0).Length(); } }; diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index 86a5ebc2..4a4e021d 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -1456,14 +1456,29 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) py::class_ (m, "Solid"); py::class_ (m, "Compound") - .def(py::init([](std::vector shapes) { + .def(py::init([](std::vector 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) ; From 842780feec723113482f16b4c1a1b4a1c173fe5d Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Mon, 14 Mar 2022 14:20:31 +0100 Subject: [PATCH 2/2] if only 1 locahl tree is set use the global one on all layers --- libsrc/meshing/meshclass.cpp | 10 ++++++---- libsrc/meshing/meshclass.hpp | 7 ++++++- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 21822503..855fe169 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -3345,10 +3345,11 @@ namespace netgen double Mesh :: GetH (const Point3d & p, int layer) const { + const auto& lh = GetLocalH(layer); double hmin = hglob; - if (lochfunc[layer-1]) + if (lh) { - double hl = lochfunc[layer-1]->GetH (p); + double hl = lh->GetH (p); if (hl < hglob) hmin = hl; } @@ -3357,10 +3358,11 @@ namespace netgen double Mesh :: GetMinH (const Point3d & pmin, const Point3d & pmax, int layer) { + const auto& lh = GetLocalH(layer); double hmin = hglob; - if (lochfunc[layer-1]) + if (lh) { - double hl = lochfunc[layer-1]->GetMinH (pmin, pmax); + double hl = lh->GetMinH (pmin, pmax); if (hl < hmin) hmin = hl; } diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index e0b10081..01197895 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -469,7 +469,12 @@ namespace netgen /// LocalH & LocalHFunction (int layer=1) { return * lochfunc[layer-1]; } - shared_ptr GetLocalH(int layer=1) const { return lochfunc[layer-1]; } + shared_ptr GetLocalH(int layer=1) const + { + if(lochfunc.Size() == 1) + return lochfunc[0]; + return lochfunc[layer-1]; + } void SetLocalH(shared_ptr loch, int layer=1); ///