ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-27 14:20:34 +05:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'blayer_as_meshing_parameter' into 'master'

Browse Source 
Add BoundarylayerParameters to MeshingParameters

See merge request ngsolve/netgen!677
This commit is contained in:
Schöberl, Joachim 2024-09-27 16:30:50 +02:00
commit 2ff62bc283
9 changed files with 317 additions and 161 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
libsrc/meshing/basegeom.cpp
View File

@ -1245,8 +1245,7 @@ namespace netgen
void NetgenGeometry :: FinalizeMesh(Mesh& mesh) const void NetgenGeometry :: FinalizeMesh(Mesh& mesh) const
{ {
if(solids.Size()) for (int i = 0; i < std::min(solids.Size(), (size_t)mesh.GetNDomains()); i++)
for (int i = 0; i < mesh.GetNDomains(); i++)
if (auto name = solids[i]->properties.name) if (auto name = solids[i]->properties.name)
mesh.SetMaterial (i+1, *name); mesh.SetMaterial (i+1, *name);

128
libsrc/meshing/boundarylayer.cpp
View File

@ -267,7 +267,7 @@ namespace netgen
for(SurfaceElementIndex sei : mesh.SurfaceElements().Range()) for(SurfaceElementIndex sei : mesh.SurfaceElements().Range())
{ {
auto facei = mesh[sei].GetIndex(); auto facei = mesh[sei].GetIndex();
if(facei < nfd_old && !params.surfid.Contains(facei)) if(facei < nfd_old && !par_surfid.Contains(facei))
continue; continue;
auto sel = mesh[sei]; auto sel = mesh[sei];
@ -322,7 +322,7 @@ namespace netgen
for(auto pi : sel.PNums()) for(auto pi : sel.PNums())
box.Add(mesh[pi]); box.Add(mesh[pi]);
// also add moved points to bounding box // also add moved points to bounding box
if(params.surfid.Contains(sel.GetIndex())) if(par_surfid.Contains(sel.GetIndex()))
for(auto pi : sel.PNums()) for(auto pi : sel.PNums())
box.Add(mesh[pi]+limits[pi]*height*growthvectors[pi]); box.Add(mesh[pi]+limits[pi]*height*growthvectors[pi]);
tree.Insert(box, sei); tree.Insert(box, sei);
@ -346,7 +346,7 @@ namespace netgen
return false; return false;
auto face = GetMappedFace(sei, -2); auto face = GetMappedFace(sei, -2);
double lam_ = 999; double lam_ = 999;
bool is_bl_sel = params.surfid.Contains(sel.GetIndex()); bool is_bl_sel = par_surfid.Contains(sel.GetIndex());
if (step == 0) if (step == 0)
{ {
@ -521,7 +521,7 @@ namespace netgen
for(SurfaceElementIndex sei : myrange) for(SurfaceElementIndex sei : myrange)
{ {
auto facei = mesh[sei].GetIndex(); 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()) for(auto pi : mesh[sei].PNums())
if(mesh[pi].Type() == SURFACEPOINT) if(mesh[pi].Type() == SURFACEPOINT)
@ -593,12 +593,13 @@ namespace netgen
{ {
static Timer timer("BoundaryLayerTool::ctor"); static Timer timer("BoundaryLayerTool::ctor");
RegionTimer regt(timer); RegionTimer regt(timer);
ProcessParameters();
//for(auto & seg : mesh.LineSegments()) //for(auto & seg : mesh.LineSegments())
//seg.edgenr = seg.epgeominfo[1].edgenr; //seg.edgenr = seg.epgeominfo[1].edgenr;
height = 0.0; height = 0.0;
for (auto h : params.heights) for (auto h : par_heights)
height += h; height += h;
max_edge_nr = -1; max_edge_nr = -1;
@ -611,7 +612,7 @@ namespace netgen
new_mat_nrs.SetSize(mesh.FaceDescriptors().Size() + 1); new_mat_nrs.SetSize(mesh.FaceDescriptors().Size() + 1);
new_mat_nrs = -1; new_mat_nrs = -1;
for(auto [bcname, matname] : params.new_mat) for(auto [bcname, matname] : par_new_mat)
{ {
mesh.SetMaterial(++ndom, matname); mesh.SetMaterial(++ndom, matname);
regex pattern(bcname); regex pattern(bcname);
@ -623,7 +624,6 @@ namespace netgen
} }
} }
domains = params.domains;
if(!params.outside) if(!params.outside)
domains.Invert(); domains.Invert();
@ -660,7 +660,7 @@ namespace netgen
{ {
const auto& fd = mesh.GetFaceDescriptor(i); const auto& fd = mesh.GetFaceDescriptor(i);
string name = fd.GetBCName(); 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())) 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) if(surfacefacs[si] == 0.0)
throw Exception("Surface " + to_string(si) + " is not a boundary of the domain to be grown into!"); 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]; const auto & sel = mesh[sei];
auto facei = sel.GetIndex(); auto facei = sel.GetIndex();
if(!params.surfid.Contains(facei)) if(!par_surfid.Contains(facei))
continue; continue;
auto n = surfacefacs[sel.GetIndex()] * getNormal(sel); auto n = surfacefacs[sel.GetIndex()] * getNormal(sel);
@ -943,7 +943,7 @@ namespace netgen
else else
continue; continue;
if(!params.project_boundaries.Contains(sel.GetIndex())) if(!par_project_boundaries.Contains(sel.GetIndex()))
continue; continue;
auto& g = growthvectors[pi]; auto& g = growthvectors[pi];
auto ng = n * g; auto ng = n * g;
@ -962,7 +962,7 @@ namespace netgen
{ {
int count = 0; int count = 0;
for(const auto& seg2 : segments) 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++; count++;
if(count == 1) if(count == 1)
{ {
@ -1099,9 +1099,9 @@ namespace netgen
if (growthvectors[pi].Length2() != 0) if (growthvectors[pi].Length2() != 0)
{ {
Point<3> p = mesh[pi]; 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)); mapto[pi].Append(mesh.AddPoint(p));
} }
} }
@ -1152,7 +1152,7 @@ namespace netgen
s0.si = segj.si; s0.si = segj.si;
new_segments.Append(s0); new_segments.Append(s0);
for(auto i : Range(params.heights)) for(auto i : Range(par_heights))
{ {
Element2d sel(QUAD); Element2d sel(QUAD);
p3 = mapto[pp2][i]; p3 = mapto[pp2][i];
@ -1220,7 +1220,7 @@ namespace netgen
{ {
Array<PointIndex> points(sel.PNums()); Array<PointIndex> points(sel.PNums());
if(surfacefacs[sel.GetIndex()] > 0) Swap(points[0], points[2]); 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; auto eltype = points.Size() == 3 ? PRISM : HEX;
Element el(eltype); Element el(eltype);
@ -1328,7 +1328,7 @@ namespace netgen
PointIndex p4 = p1; PointIndex p4 = p1;
PointIndex p5 = p2; PointIndex p5 = p2;
PointIndex p6 = p3; PointIndex p6 = p3;
for(auto i : Range(params.heights)) for(auto i : Range(par_heights))
{ {
Element nel(PRISM); Element nel(PRISM);
nel[0] = p4; nel[1] = p5; nel[2] = p6; nel[0] = p4; nel[1] = p5; nel[2] = p6;
@ -1344,7 +1344,7 @@ namespace netgen
{ {
PointIndex p1 = moved[0]; PointIndex p1 = moved[0];
PointIndex p2 = moved[1]; PointIndex p2 = moved[1];
for(auto i : Range(params.heights)) for(auto i : Range(par_heights))
{ {
PointIndex p3 = mapto[moved[1]][i]; PointIndex p3 = mapto[moved[1]][i];
PointIndex p4 = mapto[moved[0]][i]; PointIndex p4 = mapto[moved[0]][i];
@ -1366,7 +1366,7 @@ namespace netgen
if(moved.Size() == 1 && fixed.Size() == 1) if(moved.Size() == 1 && fixed.Size() == 1)
{ {
PointIndex p1 = moved[0]; PointIndex p1 = moved[0];
for(auto i : Range(params.heights)) for(auto i : Range(par_heights))
{ {
Element nel = el; Element nel = el;
PointIndex p2 = mapto[moved[0]][i]; 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())); throw Exception("This case is not implemented yet! Fixed size = " + ToString(fixed.Size()));
PointIndex p1 = moved[0]; PointIndex p1 = moved[0];
PointIndex p2 = moved[1]; PointIndex p2 = moved[1];
for(auto i : Range(params.heights)) for(auto i : Range(par_heights))
{ {
PointIndex p3 = mapto[moved[1]][i]; PointIndex p3 = mapto[moved[1]][i];
PointIndex p4 = mapto[moved[0]][i]; PointIndex p4 = mapto[moved[0]][i];
@ -1517,6 +1517,94 @@ namespace netgen
} }
} }
void BoundaryLayerTool :: ProcessParameters()
{
if(int* bc = get_if<int>(&params.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<string>(&params.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<string>(&params.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<std::vector<int>>(&params.boundary);
for(auto id : surfids)
par_surfid.Append(id);
}
if(string* mat = get_if<string>(&params.new_material); mat)
par_new_mat = { { ".*", *mat } };
else
par_new_mat = *get_if<map<string, string>>(&params.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<double>(&params.thickness); height)
{
par_heights.Append(*height);
}
else
{
auto & heights = *get_if<std::vector<double>>(&params.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<string>(&params.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<int>(&params.domain);
domains.SetBit(idomain);
}
}
void BoundaryLayerTool :: Perform() void BoundaryLayerTool :: Perform()
{ {
CreateNewFaceDescriptors(); CreateNewFaceDescriptors();

24
libsrc/meshing/boundarylayer.hpp
View File

@ -10,23 +10,6 @@ DLL_HEADER extern void InsertVirtualBoundaryLayer (Mesh & mesh);
/// Create a typical prismatic boundary layer on the given /// Create a typical prismatic boundary layer on the given
/// surfaces /// surfaces
class BoundaryLayerParameters
{
public:
// parameters by Philippose ..
Array<int> surfid;
Array<double> heights;
map<string, string> 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<size_t> project_boundaries;
};
DLL_HEADER void GenerateBoundaryLayer (Mesh & mesh, DLL_HEADER void GenerateBoundaryLayer (Mesh & mesh,
const BoundaryLayerParameters & blp); const BoundaryLayerParameters & blp);
@ -36,6 +19,7 @@ class BoundaryLayerTool
{ {
public: public:
BoundaryLayerTool(Mesh & mesh_, const BoundaryLayerParameters & params_); BoundaryLayerTool(Mesh & mesh_, const BoundaryLayerParameters & params_);
void ProcessParameters();
void Perform(); void Perform();
protected: protected:
@ -53,6 +37,12 @@ class BoundaryLayerTool
int np, nseg, nse, ne; int np, nseg, nse, ne;
double height; double height;
// These parameters are derived from given BoundaryLayerParameters and the Mesh
Array<double> par_heights;
Array<int> par_surfid;
map<string, string> par_new_mat;
Array<size_t> par_project_boundaries;
bool have_single_segments; bool have_single_segments;
Array<Segment> segments, new_segments; Array<Segment> segments, new_segments;

4
libsrc/meshing/meshfunc.cpp
View File

@ -592,8 +592,8 @@ namespace netgen
static Timer t("MeshVolume"); RegionTimer reg(t); static Timer t("MeshVolume"); RegionTimer reg(t);
mesh3d.Compress(); mesh3d.Compress();
for (auto bl : mp.boundary_layers)
GenerateBoundaryLayer(mesh3d, bl);
if(mesh3d.GetNDomains()==0) if(mesh3d.GetNDomains()==0)
return MESHING3_OK; return MESHING3_OK;

63
libsrc/meshing/meshtype.cpp
View File

@ -2905,6 +2905,69 @@ namespace netgen
ost << "table: " << endl << idpoints_table << endl; 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 () MeshingParameters :: MeshingParameters ()
{ {

22
libsrc/meshing/meshtype.hpp
View File

@ -8,6 +8,8 @@
/* Date: 01. Okt. 95 */ /* Date: 01. Okt. 95 */
/**************************************************************************/ /**************************************************************************/
#include <variant>
#include <mydefs.hpp> #include <mydefs.hpp>
#include <general/template.hpp> #include <general/template.hpp>
#include <core/mpi_wrapper.hpp> #include <core/mpi_wrapper.hpp>
@ -1272,6 +1274,24 @@ namespace netgen
}; };
struct BoundaryLayerParameters
{
std::variant<string, int, std::vector<int>> boundary;
std::variant<double, std::vector<double>> thickness;
std::variant<string, std::map<string, string>> new_material;
std::variant<string, int> domain;
bool outside;
std::optional<string> 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 class DLL_HEADER MeshingParameters
{ {
@ -1397,6 +1417,8 @@ namespace netgen
int nthreads = 4; int nthreads = 4;
Flags geometrySpecificParameters; Flags geometrySpecificParameters;
Array<BoundaryLayerParameters> boundary_layers;
/// ///
MeshingParameters (); MeshingParameters ();
/// ///

210
libsrc/meshing/python_mesh.cpp
View File

@ -1469,8 +1469,8 @@ py::arg("point_tolerance") = -1.)
.def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard<py::gil_scoped_release>()) .def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard<py::gil_scoped_release>())
.def ("BoundaryLayer2", GenerateBoundaryLayer2, py::arg("domain"), py::arg("thicknesses"), py::arg("make_new_domain")=true, py::arg("boundaries")=Array<int>{}) .def ("BoundaryLayer2", GenerateBoundaryLayer2, py::arg("domain"), py::arg("thicknesses"), py::arg("make_new_domain")=true, py::arg("boundaries")=Array<int>{})
.def ("BoundaryLayer", [](Mesh & self, variant<string, int> boundary, .def ("BoundaryLayer", [](Mesh & self, variant<string, int, std::vector<int>> boundary,
variant<double, py::list> thickness, variant<double, std::vector<double>> thickness,
variant<string, map<string, string>> material, variant<string, map<string, string>> material,
variant<string, int> domain, bool outside, variant<string, int> domain, bool outside,
optional<string> project_boundaries, optional<string> project_boundaries,
@ -1479,127 +1479,22 @@ py::arg("point_tolerance") = -1.)
bool keep_surfaceindex) bool keep_surfaceindex)
{ {
BoundaryLayerParameters blp; BoundaryLayerParameters blp;
BitArray boundaries(self.GetNFD()+1); blp.boundary = boundary;
boundaries.Clear(); blp.thickness = thickness;
if(int* bc = get_if<int>(&boundary); bc) blp.new_material = material;
{ blp.domain = domain;
for (int i = 1; i <= self.GetNFD(); i++)
if(self.GetFaceDescriptor(i).BCProperty() == *bc)
boundaries.SetBit(i);
}
else
{
regex pattern(*get_if<string>(&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<string>(&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<string>(&material); mat)
blp.new_mat = { { ".*", *mat } };
else
blp.new_mat = *get_if<map<string, string>>(&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<double>(&thickness); pthickness)
{
blp.heights.Append(*pthickness);
}
else
{
auto thicknesses = *get_if<py::list>(&thickness);
for(auto val : thicknesses)
blp.heights.Append(val.cast<double>());
}
int nr_domains = self.GetNDomains();
blp.domains.SetSize(nr_domains + 1); // one based
blp.domains.Clear();
if(string* pdomain = get_if<string>(&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<int>(&domain);
blp.domains.SetBit(idomain);
}
blp.outside = outside; blp.outside = outside;
blp.project_boundaries = project_boundaries;
blp.grow_edges = grow_edges; blp.grow_edges = grow_edges;
blp.limit_growth_vectors = limit_growth_vectors; blp.limit_growth_vectors = limit_growth_vectors;
blp.sides_keep_surfaceindex = sides_keep_surfaceindex; blp.sides_keep_surfaceindex = sides_keep_surfaceindex;
blp.keep_surfaceindex = keep_surfaceindex; blp.keep_surfaceindex = keep_surfaceindex;
GenerateBoundaryLayer (self, blp); GenerateBoundaryLayer (self, blp);
self.UpdateTopology(); self.UpdateTopology();
}, py::arg("boundary"), py::arg("thickness"), py::arg("material"), }, py::arg("boundary"), py::arg("thickness"), py::arg("material"),
py::arg("domains") = ".*", py::arg("outside") = false, 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("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("keep_surfaceindex")=false, "Add boundary layer to mesh. see help(BoundaryLayerParameters) for details.")
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_static ("EnableTableClass", [] (string name, bool set) .def_static ("EnableTableClass", [] (string name, bool set)
{ {
@ -1825,6 +1720,95 @@ project_boundaries : Optional[str] = None
m.attr("debugparam") = py::cast(&debugparam); m.attr("debugparam") = py::cast(&debugparam);
py::class_<BoundaryLayerParameters>(m, "BoundaryLayerParameters")
.def(py::init([](
std::variant<string, int, std::vector<int>> boundary,
std::variant<double, std::vector<double>> thickness,
std::variant<string, std::map<string, string>> new_material,
std::variant<string, int> domain,
bool outside,
std::optional<string> 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<BoundaryLayerParameters>();
return cls(**d).cast<BoundaryLayerParameters>();
}
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<py::dict, BoundaryLayerParameters>();
#ifdef NG_CGNS #ifdef NG_CGNS
m.def("ReadCGNSFile", &ReadCGNSFile, py::arg("filename"), py::arg("base")=1, "Read mesh and solution vectors from CGNS file"); 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"), m.def("WriteCGNSFile", &WriteCGNSFile, py::arg("mesh"), py::arg("filename"), py::arg("names"), py::arg("values"), py::arg("locations"),

8
libsrc/meshing/python_mesh.hpp
View File

@ -185,6 +185,14 @@ inline void CreateMPfromKwargs(MeshingParameters& mp, py::kwargs kwargs, bool th
mp.nthreads = py::cast<int>(kwargs.attr("pop")("nthreads")); mp.nthreads = py::cast<int>(kwargs.attr("pop")("nthreads"));
if(kwargs.contains("closeedgefac")) if(kwargs.contains("closeedgefac"))
mp.closeedgefac = py::cast<optional<double>>(kwargs.attr("pop")("closeedgefac")); mp.closeedgefac = py::cast<optional<double>>(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<BoundaryLayerParameters>(layer));
}
if(kwargs.size()) if(kwargs.size())
{ {
if(throw_if_not_all_parsed) if(throw_if_not_all_parsed)

16
ng/ngpkg.cpp
View File

@ -1108,7 +1108,7 @@ namespace netgen
// Use an array to support creation of boundary // Use an array to support creation of boundary
// layers for multiple surfaces in the future... // layers for multiple surfaces in the future...
Array<int> surfid; std::vector<int> surfid;
int surfinp = 0; int surfinp = 0;
int prismlayers = 1; int prismlayers = 1;
double hfirst = 0.01; double hfirst = 0.01;
@ -1119,13 +1119,13 @@ namespace netgen
{ {
cout << "Enter Surface ID (-1 to end list): "; cout << "Enter Surface ID (-1 to end list): ";
cin >> surfinp; 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; cout << "Selected surfaces are:" << endl;
for(auto i : Range(surfid)) for(auto i : Range(surfid.size()))
cout << "Surface " << i << ": " << surfid[i] << endl; cout << "Surface " << i << ": " << surfid[i] << endl;
cout << endl << "Enter number of prism layers: "; cout << endl << "Enter number of prism layers: ";
@ -1141,15 +1141,17 @@ namespace netgen
if(growthfactor <= 0.0) growthfactor = 0.5; if(growthfactor <= 0.0) growthfactor = 0.5;
BoundaryLayerParameters blp; BoundaryLayerParameters blp;
blp.surfid = surfid; blp.boundary = surfid;
std::vector<double> thickness;
for(auto i : Range(prismlayers)) for(auto i : Range(prismlayers))
{ {
auto layer = i+1; auto layer = i+1;
if(growthfactor == 1) if(growthfactor == 1)
blp.heights.Append(layer * hfirst); thickness.push_back(layer * hfirst);
else 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); GenerateBoundaryLayer (*mesh, blp);
return TCL_OK; return TCL_OK;
} }