diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp index bcc85153..38c793e9 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -1245,8 +1245,7 @@ namespace netgen void NetgenGeometry :: FinalizeMesh(Mesh& mesh) const { - if(solids.Size()) - for (int i = 0; i < mesh.GetNDomains(); i++) + for (int i = 0; i < std::min(solids.Size(), (size_t)mesh.GetNDomains()); i++) if (auto name = solids[i]->properties.name) mesh.SetMaterial (i+1, *name); diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index bd3cee11..1ac762a5 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -267,7 +267,7 @@ namespace netgen for(SurfaceElementIndex sei : mesh.SurfaceElements().Range()) { auto facei = mesh[sei].GetIndex(); - if(facei < nfd_old && !params.surfid.Contains(facei)) + if(facei < nfd_old && !par_surfid.Contains(facei)) continue; auto sel = mesh[sei]; @@ -322,7 +322,7 @@ namespace netgen for(auto pi : sel.PNums()) box.Add(mesh[pi]); // also add moved points to bounding box - if(params.surfid.Contains(sel.GetIndex())) + if(par_surfid.Contains(sel.GetIndex())) for(auto pi : sel.PNums()) box.Add(mesh[pi]+limits[pi]*height*growthvectors[pi]); tree.Insert(box, sei); @@ -346,7 +346,7 @@ namespace netgen return false; auto face = GetMappedFace(sei, -2); double lam_ = 999; - bool is_bl_sel = params.surfid.Contains(sel.GetIndex()); + bool is_bl_sel = par_surfid.Contains(sel.GetIndex()); if (step == 0) { @@ -521,7 +521,7 @@ namespace netgen for(SurfaceElementIndex sei : myrange) { auto facei = mesh[sei].GetIndex(); - if(facei < nfd_old && !params.surfid.Contains(facei)) + if(facei < nfd_old && !par_surfid.Contains(facei)) { for(auto pi : mesh[sei].PNums()) if(mesh[pi].Type() == SURFACEPOINT) @@ -593,12 +593,13 @@ namespace netgen { static Timer timer("BoundaryLayerTool::ctor"); RegionTimer regt(timer); + ProcessParameters(); //for(auto & seg : mesh.LineSegments()) //seg.edgenr = seg.epgeominfo[1].edgenr; height = 0.0; - for (auto h : params.heights) + for (auto h : par_heights) height += h; max_edge_nr = -1; @@ -611,7 +612,7 @@ namespace netgen new_mat_nrs.SetSize(mesh.FaceDescriptors().Size() + 1); new_mat_nrs = -1; - for(auto [bcname, matname] : params.new_mat) + for(auto [bcname, matname] : par_new_mat) { mesh.SetMaterial(++ndom, matname); regex pattern(bcname); @@ -623,7 +624,6 @@ namespace netgen } } - domains = params.domains; if(!params.outside) domains.Invert(); @@ -660,7 +660,7 @@ namespace netgen { const auto& fd = mesh.GetFaceDescriptor(i); string name = fd.GetBCName(); - if(params.surfid.Contains(i)) + if(par_surfid.Contains(i)) { if(auto isIn = domains.Test(fd.DomainIn()); isIn != domains.Test(fd.DomainOut())) { @@ -685,7 +685,7 @@ namespace netgen } } - for(auto si : params.surfid) + for(auto si : par_surfid) if(surfacefacs[si] == 0.0) throw Exception("Surface " + to_string(si) + " is not a boundary of the domain to be grown into!"); } @@ -747,7 +747,7 @@ namespace netgen { const auto & sel = mesh[sei]; auto facei = sel.GetIndex(); - if(!params.surfid.Contains(facei)) + if(!par_surfid.Contains(facei)) continue; auto n = surfacefacs[sel.GetIndex()] * getNormal(sel); @@ -943,7 +943,7 @@ namespace netgen else continue; - if(!params.project_boundaries.Contains(sel.GetIndex())) + if(!par_project_boundaries.Contains(sel.GetIndex())) continue; auto& g = growthvectors[pi]; auto ng = n * g; @@ -962,7 +962,7 @@ namespace netgen { int count = 0; for(const auto& seg2 : segments) - if(((seg[0] == seg2[0] && seg[1] == seg2[1]) || (seg[0] == seg2[1] && seg[1] == seg2[0])) && params.surfid.Contains(seg2.si)) + if(((seg[0] == seg2[0] && seg[1] == seg2[1]) || (seg[0] == seg2[1] && seg[1] == seg2[0])) && par_surfid.Contains(seg2.si)) count++; if(count == 1) { @@ -1099,9 +1099,9 @@ namespace netgen if (growthvectors[pi].Length2() != 0) { Point<3> p = mesh[pi]; - for(auto i : Range(params.heights)) + for(auto i : Range(par_heights)) { - p += params.heights[i] * growthvectors[pi]; + p += par_heights[i] * growthvectors[pi]; mapto[pi].Append(mesh.AddPoint(p)); } } @@ -1152,7 +1152,7 @@ namespace netgen s0.si = segj.si; new_segments.Append(s0); - for(auto i : Range(params.heights)) + for(auto i : Range(par_heights)) { Element2d sel(QUAD); p3 = mapto[pp2][i]; @@ -1220,7 +1220,7 @@ namespace netgen { Array points(sel.PNums()); if(surfacefacs[sel.GetIndex()] > 0) Swap(points[0], points[2]); - for(auto j : Range(params.heights)) + for(auto j : Range(par_heights)) { auto eltype = points.Size() == 3 ? PRISM : HEX; Element el(eltype); @@ -1328,7 +1328,7 @@ namespace netgen PointIndex p4 = p1; PointIndex p5 = p2; PointIndex p6 = p3; - for(auto i : Range(params.heights)) + for(auto i : Range(par_heights)) { Element nel(PRISM); nel[0] = p4; nel[1] = p5; nel[2] = p6; @@ -1344,7 +1344,7 @@ namespace netgen { PointIndex p1 = moved[0]; PointIndex p2 = moved[1]; - for(auto i : Range(params.heights)) + for(auto i : Range(par_heights)) { PointIndex p3 = mapto[moved[1]][i]; PointIndex p4 = mapto[moved[0]][i]; @@ -1366,7 +1366,7 @@ namespace netgen if(moved.Size() == 1 && fixed.Size() == 1) { PointIndex p1 = moved[0]; - for(auto i : Range(params.heights)) + for(auto i : Range(par_heights)) { Element nel = el; PointIndex p2 = mapto[moved[0]][i]; @@ -1390,7 +1390,7 @@ namespace netgen throw Exception("This case is not implemented yet! Fixed size = " + ToString(fixed.Size())); PointIndex p1 = moved[0]; PointIndex p2 = moved[1]; - for(auto i : Range(params.heights)) + for(auto i : Range(par_heights)) { PointIndex p3 = mapto[moved[1]][i]; PointIndex p4 = mapto[moved[0]][i]; @@ -1517,6 +1517,94 @@ namespace netgen } } + void BoundaryLayerTool :: ProcessParameters() + { + if(int* bc = get_if(¶ms.boundary); bc) + { + for (int i = 1; i <= mesh.GetNFD(); i++) + if(mesh.GetFaceDescriptor(i).BCProperty() == *bc) + par_surfid.Append(i); + } + else if(string* s = get_if(¶ms.boundary); s) + { + regex pattern(*s); + BitArray boundaries(mesh.GetNFD()+1); + boundaries.Clear(); + for(int i = 1; i<=mesh.GetNFD(); i++) + { + auto& fd = mesh.GetFaceDescriptor(i); + if(regex_match(fd.GetBCName(), pattern)) + { + boundaries.SetBit(i); + auto dom_pattern = get_if(¶ms.domain); + // only add if adjacent to domain + if(dom_pattern) + { + regex pattern(*dom_pattern); + bool mat1_match = fd.DomainIn() > 0 && regex_match(mesh.GetMaterial(fd.DomainIn()), pattern); + bool mat2_match = fd.DomainOut() > 0 && regex_match(mesh.GetMaterial(fd.DomainOut()), pattern); + // if boundary is inner or outer remove from list + if(mat1_match == mat2_match) + boundaries.Clear(i); + // if((fd.DomainIn() > 0 && regex_match(mesh.GetMaterial(fd.DomainIn()), pattern)) || (fd.DomainOut() > 0 && regex_match(self.GetMaterial(fd.DomainOut()), pattern))) + // boundaries.Clear(i); + // par_surfid.Append(i); + } + // else + // par_surfid.Append(i); + } + } + for(int i = 1; i<=mesh.GetNFD(); i++) + if(boundaries.Test(i)) + par_surfid.Append(i); + } + else + { + auto & surfids = *get_if>(¶ms.boundary); + for(auto id : surfids) + par_surfid.Append(id); + } + if(string* mat = get_if(¶ms.new_material); mat) + par_new_mat = { { ".*", *mat } }; + else + par_new_mat = *get_if>(¶ms.new_material); + + if(params.project_boundaries.has_value()) + { + regex pattern(*params.project_boundaries); + for(int i = 1; i<=mesh.GetNFD(); i++) + if(regex_match(mesh.GetFaceDescriptor(i).GetBCName(), pattern)) + par_project_boundaries.Append(i); + } + + if(double* height = get_if(¶ms.thickness); height) + { + par_heights.Append(*height); + } + else + { + auto & heights = *get_if>(¶ms.thickness); + for(auto val : heights) + par_heights.Append(val); + } + + int nr_domains = mesh.GetNDomains(); + domains.SetSize(nr_domains + 1); // one based + domains.Clear(); + if(string* pdomain = get_if(¶ms.domain); pdomain) + { + regex pattern(*pdomain); + for(auto i : Range(1, nr_domains+1)) + if(regex_match(mesh.GetMaterial(i), pattern)) + domains.SetBit(i); + } + else + { + auto idomain = *get_if(¶ms.domain); + domains.SetBit(idomain); + } + } + void BoundaryLayerTool :: Perform() { CreateNewFaceDescriptors(); diff --git a/libsrc/meshing/boundarylayer.hpp b/libsrc/meshing/boundarylayer.hpp index c879e7bb..264819bd 100644 --- a/libsrc/meshing/boundarylayer.hpp +++ b/libsrc/meshing/boundarylayer.hpp @@ -10,23 +10,6 @@ DLL_HEADER extern void InsertVirtualBoundaryLayer (Mesh & mesh); /// Create a typical prismatic boundary layer on the given /// surfaces -class BoundaryLayerParameters -{ -public: - // parameters by Philippose .. - Array surfid; - Array heights; - map new_mat; - BitArray domains; - bool outside = false; // set the boundary layer on the outside - bool grow_edges = false; - bool limit_growth_vectors = true; - double limit_safety = 0.3; // alloow only 30% of the growth vector length - bool sides_keep_surfaceindex = false; - bool keep_surfaceindex = false; - Array project_boundaries; -}; - DLL_HEADER void GenerateBoundaryLayer (Mesh & mesh, const BoundaryLayerParameters & blp); @@ -36,6 +19,7 @@ class BoundaryLayerTool { public: BoundaryLayerTool(Mesh & mesh_, const BoundaryLayerParameters & params_); + void ProcessParameters(); void Perform(); protected: @@ -53,6 +37,12 @@ class BoundaryLayerTool int np, nseg, nse, ne; double height; + // These parameters are derived from given BoundaryLayerParameters and the Mesh + Array par_heights; + Array par_surfid; + map par_new_mat; + Array par_project_boundaries; + bool have_single_segments; Array segments, new_segments; diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index f187abfa..c551301b 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -592,8 +592,8 @@ namespace netgen static Timer t("MeshVolume"); RegionTimer reg(t); mesh3d.Compress(); - - + for (auto bl : mp.boundary_layers) + GenerateBoundaryLayer(mesh3d, bl); if(mesh3d.GetNDomains()==0) return MESHING3_OK; diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index 5d15dd1b..a14d9ca5 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -2905,6 +2905,69 @@ namespace netgen ost << "table: " << endl << idpoints_table << endl; } + ostream & operator<< (ostream & ost, const BoundaryLayerParameters & mp) + { + ost << "BoundaryLayerParameters" << endl; + ost << " boundary: "; + switch(mp.boundary.index()) + { + case 0: + ost << std::get<0>(mp.boundary); + break; + case 1: + ost << std::get<1>(mp.boundary); + break; + case 2: + ost << "["; + for (auto val : std::get<2>(mp.boundary)) + ost << val << " "; + ost << "]"; + break; + } + ost << "\n thickness: "; + switch(mp.thickness.index()) + { + case 0: + ost << std::get<0>(mp.thickness); + break; + case 1: + ost << "["; + for (auto val : std::get<1>(mp.thickness)) + ost << val << " "; + ost << "]"; + break; + } + ost <<"\n new_material: "; + switch(mp.new_material.index()) + { + case 0: + ost << std::get<0>(mp.new_material); + break; + case 1: + for (const auto & [key, value] : std::get<1>(mp.new_material)) + ost << key << " -> " << value << ", "; + break; + } + ost << "\n domain: "; + switch(mp.domain.index()) + { + case 0: + ost << std::get<0>(mp.domain); + break; + case 1: + ost << std::get<1>(mp.domain); + break; + } + ost << "\n outside: " << mp.outside; + ost << "\n project_boundaries: " << (mp.project_boundaries ? ToString(*mp.project_boundaries) : "nullopt"); + ost << "\n grow_edges: " << mp.grow_edges; + ost << "\n limit_growth_vectors: " << mp.limit_growth_vectors; + ost << "\n sides_keep_surfaceindex: " << mp.sides_keep_surfaceindex; + ost << "\n keep_surfaceindex: " << mp.keep_surfaceindex; + ost << "\n limit_safety: " << mp.limit_safety; + ost << endl; + return ost; + } MeshingParameters :: MeshingParameters () { diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 67928a91..ec189ad6 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -8,6 +8,8 @@ /* Date: 01. Okt. 95 */ /**************************************************************************/ +#include + #include #include #include @@ -1272,6 +1274,24 @@ namespace netgen }; + struct BoundaryLayerParameters + { + std::variant> boundary; + std::variant> thickness; + std::variant> new_material; + std::variant domain; + bool outside; + std::optional project_boundaries; + bool grow_edges; + bool limit_growth_vectors; + bool sides_keep_surfaceindex; + bool keep_surfaceindex; + + double limit_safety = 0.3; // alloow only 30% of the growth vector length + }; + + + ostream & operator<< (ostream & ost, const BoundaryLayerParameters & mp); class DLL_HEADER MeshingParameters { @@ -1397,6 +1417,8 @@ namespace netgen int nthreads = 4; Flags geometrySpecificParameters; + + Array boundary_layers; /// MeshingParameters (); /// diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 3404fb16..6937c905 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -1469,8 +1469,8 @@ py::arg("point_tolerance") = -1.) .def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard()) .def ("BoundaryLayer2", GenerateBoundaryLayer2, py::arg("domain"), py::arg("thicknesses"), py::arg("make_new_domain")=true, py::arg("boundaries")=Array{}) - .def ("BoundaryLayer", [](Mesh & self, variant boundary, - variant thickness, + .def ("BoundaryLayer", [](Mesh & self, variant> boundary, + variant> thickness, variant> material, variant domain, bool outside, optional project_boundaries, @@ -1479,127 +1479,22 @@ py::arg("point_tolerance") = -1.) bool keep_surfaceindex) { BoundaryLayerParameters blp; - BitArray boundaries(self.GetNFD()+1); - boundaries.Clear(); - if(int* bc = get_if(&boundary); bc) - { - for (int i = 1; i <= self.GetNFD(); i++) - if(self.GetFaceDescriptor(i).BCProperty() == *bc) - boundaries.SetBit(i); - } - else - { - regex pattern(*get_if(&boundary)); - for(int i = 1; i<=self.GetNFD(); i++) - { - auto& fd = self.GetFaceDescriptor(i); - if(regex_match(fd.GetBCName(), pattern)) - { - boundaries.SetBit(i); - auto dom_pattern = get_if(&domain); - // only add if adjacent to domain - if(dom_pattern) - { - regex pattern(*dom_pattern); - bool mat1_match = fd.DomainIn() > 0 && regex_match(self.GetMaterial(fd.DomainIn()), pattern); - bool mat2_match = fd.DomainOut() > 0 && regex_match(self.GetMaterial(fd.DomainOut()), pattern); - // if boundary is inner or outer remove from list - if(mat1_match == mat2_match) - boundaries.Clear(i); - // if((fd.DomainIn() > 0 && regex_match(self.GetMaterial(fd.DomainIn()), pattern)) || (fd.DomainOut() > 0 && regex_match(self.GetMaterial(fd.DomainOut()), pattern))) - // boundaries.Clear(i); - // blp.surfid.Append(i); - } - // else - // blp.surfid.Append(i); - } - } - } - for(int i = 1; i<=self.GetNFD(); i++) - if(boundaries.Test(i)) - blp.surfid.Append(i); - if(string* mat = get_if(&material); mat) - blp.new_mat = { { ".*", *mat } }; - else - blp.new_mat = *get_if>(&material); - - if(project_boundaries.has_value()) - { - regex pattern(*project_boundaries); - for(int i = 1; i<=self.GetNFD(); i++) - if(regex_match(self.GetFaceDescriptor(i).GetBCName(), pattern)) - blp.project_boundaries.Append(i); - } - - if(double* pthickness = get_if(&thickness); pthickness) - { - blp.heights.Append(*pthickness); - } - else - { - auto thicknesses = *get_if(&thickness); - for(auto val : thicknesses) - blp.heights.Append(val.cast()); - } - - int nr_domains = self.GetNDomains(); - blp.domains.SetSize(nr_domains + 1); // one based - blp.domains.Clear(); - if(string* pdomain = get_if(&domain); pdomain) - { - regex pattern(*pdomain); - for(auto i : Range(1, nr_domains+1)) - if(regex_match(self.GetMaterial(i), pattern)) - blp.domains.SetBit(i); - } - else - { - auto idomain = *get_if(&domain); - blp.domains.SetBit(idomain); - } - + blp.boundary = boundary; + blp.thickness = thickness; + blp.new_material = material; + blp.domain = domain; blp.outside = outside; + blp.project_boundaries = project_boundaries; blp.grow_edges = grow_edges; blp.limit_growth_vectors = limit_growth_vectors; blp.sides_keep_surfaceindex = sides_keep_surfaceindex; blp.keep_surfaceindex = keep_surfaceindex; - GenerateBoundaryLayer (self, blp); self.UpdateTopology(); }, py::arg("boundary"), py::arg("thickness"), py::arg("material"), py::arg("domains") = ".*", py::arg("outside") = false, py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true, py::arg("limit_growth_vectors") = true, py::arg("sides_keep_surfaceindex")=false, - py::arg("keep_surfaceindex")=false, - R"delimiter( -Add boundary layer to mesh. - -Parameters ----------- - -boundary : string or int - Boundary name or number. - -thickness : float or List[float] - Thickness of boundary layer(s). - -material : str or List[str] - Material name of boundary layer(s). - -domain : str or int - Regexp for domain boundarylayer is going into. - -outside : bool = False - If true add the layer on the outside - -grow_edges : bool = False - Grow boundary layer over edges. - -project_boundaries : Optional[str] = None - Project boundarylayer to these boundaries if they meet them. Set - to boundaries that meet boundarylayer at a non-orthogonal edge and - layer-ending should be projected to that boundary. - -)delimiter") + py::arg("keep_surfaceindex")=false, "Add boundary layer to mesh. see help(BoundaryLayerParameters) for details.") .def_static ("EnableTableClass", [] (string name, bool set) { @@ -1825,6 +1720,95 @@ project_boundaries : Optional[str] = None m.attr("debugparam") = py::cast(&debugparam); + py::class_(m, "BoundaryLayerParameters") + .def(py::init([]( + std::variant> boundary, + std::variant> thickness, + std::variant> new_material, + std::variant domain, + bool outside, + std::optional project_boundaries, + bool grow_edges, + bool limit_growth_vectors, + bool sides_keep_surfaceindex, + bool keep_surfaceindex, + double limit_safety) + { + BoundaryLayerParameters blp; + blp.boundary = boundary; + blp.thickness = thickness; + blp.new_material = new_material; + blp.domain = domain; + blp.outside = outside; + blp.project_boundaries = project_boundaries; + blp.grow_edges = grow_edges; + blp.limit_growth_vectors = limit_growth_vectors; + blp.sides_keep_surfaceindex = sides_keep_surfaceindex; + blp.keep_surfaceindex = keep_surfaceindex; + blp.limit_safety = limit_safety; + return blp; + }), + py::arg("boundary"), py::arg("thickness"), py::arg("new_material"), + py::arg("domain") = ".*", py::arg("outside") = false, + py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true, + py::arg("limit_growth_vectors") = true, py::arg("sides_keep_surfaceindex")=false, + py::arg("keep_surfaceindex")=false, py::arg("limit_safety")=0.3, + R"delimiter( +Add boundary layer to mesh. + +Parameters +---------- + +boundary : string or int + Boundary name or number. + +thickness : float or List[float] + Thickness of boundary layer(s). + +material : str or List[str] + Material name of boundary layer(s). + +domain : str or int + Regexp for domain boundarylayer is going into. + +outside : bool = False + If true add the layer on the outside + +grow_edges : bool = False + Grow boundary layer over edges. + +project_boundaries : Optional[str] = None + Project boundarylayer to these boundaries if they meet them. Set + to boundaries that meet boundarylayer at a non-orthogonal edge and + layer-ending should be projected to that boundary. + +)delimiter") + .def(py::init([]( const py::dict & d ) { + try { + // Call other constructor with named arguments by unpacking the dictionary + py::object cls = py::type::of(); + return cls(**d).cast(); + } + catch (py::error_already_set & e) { + cerr << "Error creating BoundaryLayerParameters from dict:" << endl; + cerr << e.what() << endl; + throw; + } + })) + .def_readwrite("boundary", &BoundaryLayerParameters::boundary) + .def_readwrite("thickness", &BoundaryLayerParameters::thickness) + .def_readwrite("new_material", &BoundaryLayerParameters::new_material) + .def_readwrite("domain", &BoundaryLayerParameters::domain) + .def_readwrite("outside", &BoundaryLayerParameters::outside) + .def_readwrite("project_boundaries", &BoundaryLayerParameters::project_boundaries) + .def_readwrite("grow_edges", &BoundaryLayerParameters::grow_edges) + .def_readwrite("limit_growth_vectors", &BoundaryLayerParameters::limit_growth_vectors) + .def_readwrite("sides_keep_surfaceindex", &BoundaryLayerParameters::sides_keep_surfaceindex) + .def_readwrite("keep_surfaceindex", &BoundaryLayerParameters::keep_surfaceindex) + .def_readwrite("limit_safety", &BoundaryLayerParameters::limit_safety) + ; + py::implicitly_convertible(); + #ifdef NG_CGNS m.def("ReadCGNSFile", &ReadCGNSFile, py::arg("filename"), py::arg("base")=1, "Read mesh and solution vectors from CGNS file"); m.def("WriteCGNSFile", &WriteCGNSFile, py::arg("mesh"), py::arg("filename"), py::arg("names"), py::arg("values"), py::arg("locations"), diff --git a/libsrc/meshing/python_mesh.hpp b/libsrc/meshing/python_mesh.hpp index 4a5db879..d0186a1f 100644 --- a/libsrc/meshing/python_mesh.hpp +++ b/libsrc/meshing/python_mesh.hpp @@ -185,6 +185,14 @@ inline void CreateMPfromKwargs(MeshingParameters& mp, py::kwargs kwargs, bool th mp.nthreads = py::cast(kwargs.attr("pop")("nthreads")); if(kwargs.contains("closeedgefac")) mp.closeedgefac = py::cast>(kwargs.attr("pop")("closeedgefac")); + + if(kwargs.contains("boundary_layers")) + { + auto layers = py::list(kwargs.attr("pop")("boundary_layers")); + for(auto layer : layers) + mp.boundary_layers.Append(py::cast(layer)); + } + if(kwargs.size()) { if(throw_if_not_all_parsed) diff --git a/ng/ngpkg.cpp b/ng/ngpkg.cpp index e7120710..8b6a245b 100644 --- a/ng/ngpkg.cpp +++ b/ng/ngpkg.cpp @@ -1108,7 +1108,7 @@ namespace netgen // Use an array to support creation of boundary // layers for multiple surfaces in the future... - Array surfid; + std::vector surfid; int surfinp = 0; int prismlayers = 1; double hfirst = 0.01; @@ -1119,13 +1119,13 @@ namespace netgen { cout << "Enter Surface ID (-1 to end list): "; cin >> surfinp; - if(surfinp >= 0) surfid.Append(surfinp); + if(surfinp >= 0) surfid.push_back(surfinp); } - cout << "Number of surfaces entered = " << surfid.Size() << endl; + cout << "Number of surfaces entered = " << surfid.size() << endl; cout << "Selected surfaces are:" << endl; - for(auto i : Range(surfid)) + for(auto i : Range(surfid.size())) cout << "Surface " << i << ": " << surfid[i] << endl; cout << endl << "Enter number of prism layers: "; @@ -1141,15 +1141,17 @@ namespace netgen if(growthfactor <= 0.0) growthfactor = 0.5; BoundaryLayerParameters blp; - blp.surfid = surfid; + blp.boundary = surfid; + std::vector thickness; for(auto i : Range(prismlayers)) { auto layer = i+1; if(growthfactor == 1) - blp.heights.Append(layer * hfirst); + thickness.push_back(layer * hfirst); else - blp.heights.Append(hfirst * (pow(growthfactor, (layer+1))-1)/(growthfactor-1)); + thickness.push_back(hfirst * (pow(growthfactor, (layer+1))-1)/(growthfactor-1)); } + blp.thickness = thickness; GenerateBoundaryLayer (*mesh, blp); return TCL_OK; }