Apply blayer_as_meshing_parameter changes

This commit is contained in:
Matthias Hochsteger 2024-09-27 13:12:14 +02:00
parent 7774b6ccd7
commit 663eaffbb2
8 changed files with 301 additions and 157 deletions

View File

@ -152,7 +152,7 @@ struct GrowthVectorLimiter {
growthvectors(tool_.growthvectors),
map_from(mesh.Points().Size()),
debug("debug.txt") {
changed_domains = params.domains;
changed_domains = tool.domains;
if (!params.outside) changed_domains.Invert();
map_from = tool.mapfrom;
@ -851,12 +851,13 @@ BoundaryLayerTool::BoundaryLayerTool(Mesh& mesh_,
: mesh(mesh_), topo(mesh_.GetTopology()), params(params_) {
static Timer timer("BoundaryLayerTool::ctor");
RegionTimer regt(timer);
ProcessParameters();
// for(auto & seg : mesh.LineSegments())
// seg.edgenr = seg.epgeominfo[1].edgenr;
total_height = 0.0;
for (auto h : params.heights) total_height += h;
for (auto h : par_heights) total_height += h;
max_edge_nr = -1;
for (const auto& seg : mesh.LineSegments())
@ -867,7 +868,7 @@ BoundaryLayerTool::BoundaryLayerTool(Mesh& mesh_,
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);
for (auto i : Range(1, mesh.GetNFD() + 1)) {
@ -876,7 +877,6 @@ BoundaryLayerTool::BoundaryLayerTool(Mesh& mesh_,
}
}
domains = params.domains;
if (!params.outside) domains.Invert();
topo.SetBuildVertex2Element(true);
@ -910,7 +910,7 @@ void BoundaryLayerTool ::CreateNewFaceDescriptors() {
for (auto i : Range(1, nfd_old + 1)) {
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())) {
int new_si = mesh.GetNFD() + 1;
@ -933,7 +933,7 @@ void BoundaryLayerTool ::CreateNewFaceDescriptors() {
}
}
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!");
@ -986,7 +986,7 @@ void BoundaryLayerTool ::CalculateGrowthVectors() {
for (auto sei : p2sel[pi]) {
const auto& sel = mesh[sei];
auto facei = sel.GetIndex();
if (!params.surfid.Contains(facei)) continue;
if (!par_surfid.Contains(facei)) continue;
auto n = surfacefacs[sel.GetIndex()] * getNormal(sel);
@ -1107,7 +1107,7 @@ BitArray BoundaryLayerTool ::ProjectGrowthVectorsOnSurface() {
else
continue;
if (!params.project_boundaries.Contains(sel.GetIndex())) continue;
if (!par_project_boundaries.Contains(sel.GetIndex())) continue;
auto& g = growthvectors[pi];
auto ng = n * g;
auto gg = g * g;
@ -1124,7 +1124,7 @@ BitArray BoundaryLayerTool ::ProjectGrowthVectorsOnSurface() {
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))
par_surfid.Contains(seg2.si))
count++;
if (count == 1) {
growthvectors[seg[0]] = {0., 0., 0.};
@ -1284,8 +1284,8 @@ void BoundaryLayerTool ::InsertNewElements(
Point<3> p = mesh[pi];
PointIndex pi_last = pi;
double height = 0.0;
for (auto i : Range(params.heights)) {
height += params.heights[i];
for (auto i : Range(par_heights)) {
height += par_heights[i];
auto pi_new = mesh.AddPoint(p);
mapfrom.Append(pi);
new_points.Append(pi_new);
@ -1309,7 +1309,7 @@ void BoundaryLayerTool ::InsertNewElements(
// get point from mapto (or the group if point is mapped to multiple new
// points) layer = -1 means last point (top of boundary layer)
auto newPoint = [&](PointIndex pi, int layer = -1, int group = 0) {
if (layer == -1) layer = params.heights.Size() - 1;
if (layer == -1) layer = par_heights.Size() - 1;
if (special_boundary_points.count(pi))
return special_boundary_points[pi].growth_groups[group].new_points[layer];
else
@ -1408,7 +1408,7 @@ void BoundaryLayerTool ::InsertNewElements(
s0.si = segj.si;
new_segments.Append(s0);
for (auto i : Range(params.heights)) {
for (auto i : Range(par_heights)) {
Element2d sel(QUAD);
p3 = newPoint(pp2, i);
p4 = newPoint(pp1, i);
@ -1487,7 +1487,7 @@ void BoundaryLayerTool ::InsertNewElements(
bool add_volume_element = true;
for (auto pi : sel.PNums())
if (numGroups(pi) > 1) add_volume_element = false;
for (auto j : Range(params.heights)) {
for (auto j : Range(par_heights)) {
auto eltype = points.Size() == 3 ? PRISM : HEX;
Element el(eltype);
for (auto i : Range(points)) el[i] = points[i];
@ -1700,7 +1700,7 @@ void BoundaryLayerTool ::InsertNewElements(
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;
@ -1719,7 +1719,7 @@ void BoundaryLayerTool ::InsertNewElements(
if (fixed.Size() == 1) {
PointIndex p1 = moved[0];
PointIndex p2 = moved[1];
for (auto i : Range(params.heights)) {
for (auto i : Range(par_heights)) {
PointIndex p3 = newPoint(moved[1], i);
PointIndex p4 = newPoint(moved[0], i);
Element nel(PYRAMID);
@ -1742,7 +1742,7 @@ void BoundaryLayerTool ::InsertNewElements(
}
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 = newPoint(moved[0], i);
for (auto& p : nel.PNums()) {
@ -1762,7 +1762,7 @@ void BoundaryLayerTool ::InsertNewElements(
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 = newPoint(moved[1], i);
PointIndex p4 = newPoint(moved[0], i);
Element nel(PYRAMID);
@ -1875,6 +1875,94 @@ void BoundaryLayerTool ::FixVolumeElements() {
}
}
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() {
CreateNewFaceDescriptors();
CalculateGrowthVectors();

View File

@ -14,23 +14,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<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;
};
struct SpecialBoundaryPoint {
struct GrowthGroup {
Array<int> faces;
@ -59,6 +42,7 @@ class BoundaryLayerTool
{
public:
BoundaryLayerTool(Mesh & mesh_, const BoundaryLayerParameters & params_);
void ProcessParameters();
void Perform();
Mesh & mesh;
@ -75,6 +59,12 @@ class BoundaryLayerTool
int np, nseg, nse, ne;
double total_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;
Array<Segment> segments, new_segments;
Array<Array<PointIndex>, PointIndex> mapto;

View File

@ -593,8 +593,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;

View File

@ -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 ()
{

View File

@ -8,6 +8,8 @@
/* Date: 01. Okt. 95 */
/**************************************************************************/
#include <variant>
#include <mydefs.hpp>
#include <general/template.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
{
@ -1397,6 +1417,8 @@ namespace netgen
int nthreads = 4;
Flags geometrySpecificParameters;
Array<BoundaryLayerParameters> boundary_layers;
///
MeshingParameters ();
///

View File

@ -1469,8 +1469,8 @@ py::arg("point_tolerance") = -1.)
.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 ("BoundaryLayer", [](Mesh & self, variant<string, int> boundary,
variant<double, py::list> thickness,
.def ("BoundaryLayer", [](Mesh & self, variant<string, int, std::vector<int>> boundary,
variant<double, std::vector<double>> thickness,
variant<string, map<string, string>> material,
variant<string, int> domain, bool outside,
optional<string> 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<int>(&boundary); bc)
{
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.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,82 @@ project_boundaries : Optional[str] = None
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_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)
;
#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"),

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"));
if(kwargs.contains("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(throw_if_not_all_parsed)

View File

@ -1108,7 +1108,7 @@ namespace netgen
// Use an array to support creation of boundary
// layers for multiple surfaces in the future...
Array<int> surfid;
std::vector<int> 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<double> 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;
}