ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-02-07 02:04:19 +05:00
Code Issues Packages Projects Releases Wiki Activity

format boundarylayer.cpp

Browse Source
This commit is contained in:
Matthias Hochsteger 2024-09-27 17:30:36 +02:00
parent 2124ab5980
commit bb2646945a
1 changed files with 406 additions and 340 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

334
libsrc/meshing/boundarylayer.cpp
View File

@ -33,27 +33,32 @@ std::tuple<int, int> FindCloseVectors(FlatArray<Vec<3>> ns,
} }
Vec<3> CalcGrowthVector(FlatArray<Vec<3>> ns) { Vec<3> CalcGrowthVector(FlatArray<Vec<3>> ns) {
if (ns.Size() == 0) return {0, 0, 0}; if (ns.Size() == 0)
if (ns.Size() == 1) return ns[0]; return {0, 0, 0};
if (ns.Size() == 1)
return ns[0];
if (ns.Size() == 2) { if (ns.Size() == 2) {
auto gw = ns[0]; auto gw = ns[0];
auto n = ns[1]; auto n = ns[1];
auto npn = gw * n; auto npn = gw * n;
auto npnp = gw * gw; auto npnp = gw * gw;
auto nn = n * n; auto nn = n * n;
if (fabs(nn - npn * npn / npnp) < 1e-6) return n; if (fabs(nn - npn * npn / npnp) < 1e-6)
return n;
gw += (nn - npn) / (nn - npn * npn / npnp) * (n - npn / npnp * gw); gw += (nn - npn) / (nn - npn * npn / npnp) * (n - npn / npnp * gw);
return gw; return gw;
} }
if (ns.Size() == 3) { if (ns.Size() == 3) {
DenseMatrix mat(3, 3); DenseMatrix mat(3, 3);
for (auto i : Range(3)) for (auto i : Range(3))
for (auto j : Range(3)) mat(i, j) = ns[i][j]; for (auto j : Range(3))
mat(i, j) = ns[i][j];
if (fabs(mat.Det()) > 1e-6) { if (fabs(mat.Det()) > 1e-6) {
DenseMatrix mat(3, 3); DenseMatrix mat(3, 3);
for (auto i : Range(3)) for (auto i : Range(3))
for (auto j : Range(3)) mat(i, j) = ns[i] * ns[j]; for (auto j : Range(3))
mat(i, j) = ns[i] * ns[j];
Vector rhs(3); Vector rhs(3);
rhs = 1.; rhs = 1.;
Vector res(3); Vector res(3);
@ -61,7 +66,8 @@ Vec<3> CalcGrowthVector(FlatArray<Vec<3>> ns) {
CalcInverse(mat, inv); CalcInverse(mat, inv);
inv.Mult(rhs, res); inv.Mult(rhs, res);
Vec<3> growth = 0.; Vec<3> growth = 0.;
for (auto i : Range(ns)) growth += res[i] * ns[i]; for (auto i : Range(ns))
growth += res[i] * ns[i];
return growth; return growth;
} }
} }
@ -114,8 +120,10 @@ SpecialBoundaryPoint ::SpecialBoundaryPoint(
g1_faces.Append(facei); g1_faces.Append(facei);
g2_faces.Append(facei); g2_faces.Append(facei);
} }
for (auto fi : g1_faces) normals1.Append(normals.at(fi)); for (auto fi : g1_faces)
for (auto fi : g2_faces) normals2.Append(normals.at(fi)); normals1.Append(normals.at(fi));
for (auto fi : g2_faces)
normals2.Append(normals.at(fi));
growth_groups.Append(GrowthGroup(g1_faces, normals1)); growth_groups.Append(GrowthGroup(g1_faces, normals1));
growth_groups.Append(GrowthGroup(g2_faces, normals2)); growth_groups.Append(GrowthGroup(g2_faces, normals2));
} }
@ -125,13 +133,16 @@ Vec<3> BoundaryLayerTool ::getEdgeTangent(PointIndex pi, int edgenr) {
ArrayMem<PointIndex, 2> pts; ArrayMem<PointIndex, 2> pts;
for (auto segi : topo.GetVertexSegments(pi)) { for (auto segi : topo.GetVertexSegments(pi)) {
auto &seg = mesh[segi]; auto &seg = mesh[segi];
if (seg.edgenr != edgenr + 1) continue; if (seg.edgenr != edgenr + 1)
continue;
PointIndex other = seg[0] + seg[1] - pi; PointIndex other = seg[0] + seg[1] - pi;
if (!pts.Contains(other)) pts.Append(other); if (!pts.Contains(other))
pts.Append(other);
} }
if (pts.Size() != 2) { if (pts.Size() != 2) {
cout << "getEdgeTangent pi = " << pi << ", edgenr = " << edgenr << endl; cout << "getEdgeTangent pi = " << pi << ", edgenr = " << edgenr << endl;
for (auto segi : topo.GetVertexSegments(pi)) cout << mesh[segi] << endl; for (auto segi : topo.GetVertexSegments(pi))
cout << mesh[segi] << endl;
throw Exception("Something went wrong in getEdgeTangent!"); throw Exception("Something went wrong in getEdgeTangent!");
} }
tangent = mesh[pts[1]] - mesh[pts[0]]; tangent = mesh[pts[1]] - mesh[pts[0]];
@ -144,7 +155,6 @@ void BoundaryLayerTool ::LimitGrowthVectorLengths() {
GrowthVectorLimiter limiter(*this); GrowthVectorLimiter limiter(*this);
limiter.Perform(); limiter.Perform();
} }
// depending on the geometry type, the mesh contains segments multiple times // depending on the geometry type, the mesh contains segments multiple times
@ -155,7 +165,8 @@ bool HaveSingleSegments(const Mesh& mesh) {
for (auto segi : Range(mesh.LineSegments())) { for (auto segi : Range(mesh.LineSegments())) {
mesh.GetTopology().GetSegmentSurfaceElements(segi + 1, surf_els); mesh.GetTopology().GetSegmentSurfaceElements(segi + 1, surf_els);
if (surf_els.Size() < 2) continue; if (surf_els.Size() < 2)
continue;
auto seg = mesh[segi]; auto seg = mesh[segi];
auto pi0 = min(seg[0], seg[1]); auto pi0 = min(seg[0], seg[1]);
@ -163,12 +174,14 @@ bool HaveSingleSegments(const Mesh& mesh) {
auto p0_segs = topo.GetVertexSegments(seg[0]); auto p0_segs = topo.GetVertexSegments(seg[0]);
for (auto segi_other : p0_segs) { for (auto segi_other : p0_segs) {
if (segi_other == segi) continue; if (segi_other == segi)
continue;
auto seg_other = mesh[segi_other]; auto seg_other = mesh[segi_other];
auto pi0_other = min(seg_other[0], seg_other[1]); auto pi0_other = min(seg_other[0], seg_other[1]);
auto pi1_other = max(seg_other[0], seg_other[1]); auto pi1_other = max(seg_other[0], seg_other[1]);
if (pi0_other == pi0 && pi1_other == pi1) return false; if (pi0_other == pi0 && pi1_other == pi1)
return false;
} }
// found segment with multiple adjacent surface elements but no other // found segment with multiple adjacent surface elements but no other
@ -197,7 +210,8 @@ Array<Segment> BuildSegments(Mesh& mesh) {
auto np = sel.GetNP(); auto np = sel.GetNP();
for (auto i : Range(np)) { for (auto i : Range(np)) {
if (sel[i] == seg[0]) { if (sel[i] == seg[0]) {
if (sel[(i + 1) % np] != seg[1]) swap(seg[0], seg[1]); if (sel[(i + 1) % np] != seg[1])
swap(seg[0], seg[1]);
break; break;
} }
} }
@ -224,12 +238,15 @@ void MergeAndAddSegments(Mesh& mesh, FlatArray<Segment> segments,
} }
}; };
for (const auto& seg : segments) addSegment(seg); for (const auto &seg : segments)
addSegment(seg);
for (const auto& seg : new_segments) addSegment(seg); for (const auto &seg : new_segments)
addSegment(seg);
} }
// TODO: Hack, move this to the header or restructure the whole growth_vectors storage // TODO: Hack, move this to the header or restructure the whole growth_vectors
// storage
static std::map<PointIndex, Vec<3>> non_bl_growth_vectors; static std::map<PointIndex, Vec<3>> non_bl_growth_vectors;
void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() { void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
@ -273,8 +290,9 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
for (SurfaceElementIndex sei : myrange) { for (SurfaceElementIndex sei : myrange) {
for (auto pi : mesh[sei].PNums()) { for (auto pi : mesh[sei].PNums()) {
auto pi_from = mapfrom[pi]; auto pi_from = mapfrom[pi];
if((pi_from.IsValid() && mesh[pi_from].Type() == SURFACEPOINT) if ((pi_from.IsValid() && mesh[pi_from].Type() == SURFACEPOINT) ||
|| (!pi_from.IsValid() && mapto[pi].Size()==0 && mesh[pi].Type() == SURFACEPOINT)) (!pi_from.IsValid() && mapto[pi].Size() == 0 &&
mesh[pi].Type() == SURFACEPOINT))
points_set.insert(pi); points_set.insert(pi);
} }
} }
@ -292,7 +310,8 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
for (auto seg : segments) for (auto seg : segments)
if (has_moved_points[seg.edgenr]) if (has_moved_points[seg.edgenr])
for (auto pi : seg.PNums()) for (auto pi : seg.PNums())
if (mesh[pi].Type() == EDGEPOINT) points_set.insert(pi); if (mesh[pi].Type() == EDGEPOINT)
points_set.insert(pi);
Array<PointIndex> points; Array<PointIndex> points;
for (auto pi : points_set) for (auto pi : points_set)
@ -311,21 +330,23 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
std::set<PointIndex> suround; std::set<PointIndex> suround;
suround.insert(pi); suround.insert(pi);
// average only tangent component on new bl points, average whole growth vector otherwise // average only tangent component on new bl points, average whole growth
// vector otherwise
bool do_average_tangent = mapfrom[pi].IsValid(); bool do_average_tangent = mapfrom[pi].IsValid();
correction = 0.0; correction = 0.0;
for (auto sei : sels) { for (auto sei : sels) {
const auto &sel = mesh[sei]; const auto &sel = mesh[sei];
for (auto pi1 : sel.PNums()) { for (auto pi1 : sel.PNums()) {
if (suround.count(pi1)) continue; if (suround.count(pi1))
continue;
suround.insert(pi1); suround.insert(pi1);
auto gw_other = getGW(pi1) + corrections[pi1]; auto gw_other = getGW(pi1) + corrections[pi1];
if (do_average_tangent) { if (do_average_tangent) {
auto normal_other = getNormal(mesh[sei]); auto normal_other = getNormal(mesh[sei]);
auto tangent_part = gw_other - (gw_other * normal_other) * normal_other; auto tangent_part =
gw_other - (gw_other * normal_other) * normal_other;
correction += tangent_part; correction += tangent_part;
} } else {
else {
correction += gw_other; correction += gw_other;
} }
} }
@ -350,11 +371,13 @@ BoundaryLayerTool::BoundaryLayerTool(Mesh& mesh_,
// seg.edgenr = seg.epgeominfo[1].edgenr; // seg.edgenr = seg.epgeominfo[1].edgenr;
total_height = 0.0; total_height = 0.0;
for (auto h : par_heights) total_height += h; for (auto h : par_heights)
total_height += h;
max_edge_nr = -1; max_edge_nr = -1;
for (const auto &seg : mesh.LineSegments()) for (const auto &seg : mesh.LineSegments())
if (seg.edgenr > max_edge_nr) max_edge_nr = seg.edgenr; if (seg.edgenr > max_edge_nr)
max_edge_nr = seg.edgenr;
int ndom = mesh.GetNDomains(); int ndom = mesh.GetNDomains();
ndom_old = ndom; ndom_old = ndom;
@ -366,11 +389,13 @@ BoundaryLayerTool::BoundaryLayerTool(Mesh& mesh_,
regex pattern(bcname); regex pattern(bcname);
for (auto i : Range(1, mesh.GetNFD() + 1)) { for (auto i : Range(1, mesh.GetNFD() + 1)) {
auto &fd = mesh.GetFaceDescriptor(i); auto &fd = mesh.GetFaceDescriptor(i);
if (regex_match(fd.GetBCName(), pattern)) new_mat_nrs[i] = ndom; if (regex_match(fd.GetBCName(), pattern))
new_mat_nrs[i] = ndom;
} }
} }
if (!params.outside) domains.Invert(); if (!params.outside)
domains.Invert();
topo.SetBuildVertex2Element(true); topo.SetBuildVertex2Element(true);
mesh.UpdateTopology(); mesh.UpdateTopology();
@ -393,7 +418,8 @@ BoundaryLayerTool::BoundaryLayerTool(Mesh& mesh_,
moved_surfaces.SetSize(nfd_old + 1); moved_surfaces.SetSize(nfd_old + 1);
moved_surfaces.Clear(); moved_surfaces.Clear();
si_map.SetSize(nfd_old + 1); si_map.SetSize(nfd_old + 1);
for (auto i : Range(nfd_old + 1)) si_map[i] = i; for (auto i : Range(nfd_old + 1))
si_map[i] = i;
} }
void BoundaryLayerTool ::CreateNewFaceDescriptors() { void BoundaryLayerTool ::CreateNewFaceDescriptors() {
@ -422,7 +448,8 @@ void BoundaryLayerTool ::CreateNewFaceDescriptors() {
// result in pushed through elements, since we do not (yet) // result in pushed through elements, since we do not (yet)
// curvature through layers. // curvature through layers.
// Therefore we disable curving for these surfaces. // Therefore we disable curving for these surfaces.
if (!params.keep_surfaceindex) mesh.GetFaceDescriptor(i).SetSurfNr(-1); if (!params.keep_surfaceindex)
mesh.GetFaceDescriptor(i).SetSurfNr(-1);
} }
} }
@ -437,11 +464,13 @@ void BoundaryLayerTool ::CreateFaceDescriptorsSides() {
face_done.Clear(); face_done.Clear();
for (const auto &sel : mesh.SurfaceElements()) { for (const auto &sel : mesh.SurfaceElements()) {
auto facei = sel.GetIndex(); auto facei = sel.GetIndex();
if (face_done.Test(facei)) continue; if (face_done.Test(facei))
continue;
bool point_moved = false; bool point_moved = false;
// bool point_fixed = false; // bool point_fixed = false;
for (auto pi : sel.PNums()) { for (auto pi : sel.PNums()) {
if (growthvectors[pi].Length() > 0) point_moved = true; if (growthvectors[pi].Length() > 0)
point_moved = true;
/* /*
else else
point_fixed = true; point_fixed = true;
@ -470,7 +499,8 @@ void BoundaryLayerTool ::CalculateGrowthVectors() {
for (auto pi : mesh.Points().Range()) { for (auto pi : mesh.Points().Range()) {
const auto &p = mesh[pi]; const auto &p = mesh[pi];
if (p.Type() == INNERPOINT) continue; if (p.Type() == INNERPOINT)
continue;
std::map<int, Vec<3>> normals; std::map<int, Vec<3>> normals;
@ -479,7 +509,8 @@ void BoundaryLayerTool ::CalculateGrowthVectors() {
for (auto sei : p2sel[pi]) { for (auto sei : p2sel[pi]) {
const auto &sel = mesh[sei]; const auto &sel = mesh[sei];
auto facei = sel.GetIndex(); auto facei = sel.GetIndex();
if (!par_surfid.Contains(facei)) continue; if (!par_surfid.Contains(facei))
continue;
auto n = surfacefacs[sel.GetIndex()] * getNormal(sel); auto n = surfacefacs[sel.GetIndex()] * getNormal(sel);
@ -487,11 +518,13 @@ void BoundaryLayerTool ::CalculateGrowthVectors() {
itrig += sel.GetNP(); itrig += sel.GetNP();
auto v0 = (mesh[sel.PNumMod(itrig + 1)] - mesh[pi]).Normalize(); auto v0 = (mesh[sel.PNumMod(itrig + 1)] - mesh[pi]).Normalize();
auto v1 = (mesh[sel.PNumMod(itrig - 1)] - mesh[pi]).Normalize(); auto v1 = (mesh[sel.PNumMod(itrig - 1)] - mesh[pi]).Normalize();
if (normals.count(facei) == 0) normals[facei] = {0., 0., 0.}; if (normals.count(facei) == 0)
normals[facei] = {0., 0., 0.};
normals[facei] += acos(v0 * v1) * n; normals[facei] += acos(v0 * v1) * n;
} }
for (auto& [facei, n] : normals) n *= 1.0 / n.Length(); for (auto &[facei, n] : normals)
n *= 1.0 / n.Length();
// combine normal vectors for each face to keep uniform distances // combine normal vectors for each face to keep uniform distances
ArrayMem<Vec<3>, 5> ns; ArrayMem<Vec<3>, 5> ns;
@ -536,15 +569,18 @@ BoundaryLayerTool ::BuildSegMap() {
is_boundary_moved.Clear(); is_boundary_moved.Clear();
for (auto si : Range(segments)) { for (auto si : Range(segments)) {
if (segs_done[si]) continue; if (segs_done[si])
continue;
const auto &segi = segments[si]; const auto &segi = segments[si];
if (!moved_surfaces.Test(segi.si)) continue; if (!moved_surfaces.Test(segi.si))
continue;
segs_done.SetBit(si); segs_done.SetBit(si);
segmap[si].Append(make_pair(si, 0)); segmap[si].Append(make_pair(si, 0));
moved_segs.Append(si); moved_segs.Append(si);
is_edge_moved.SetBit(segi.edgenr); is_edge_moved.SetBit(segi.edgenr);
for (auto sj : Range(segments)) { for (auto sj : Range(segments)) {
if (segs_done.Test(sj)) continue; if (segs_done.Test(sj))
continue;
const auto &segj = segments[sj]; const auto &segj = segments[sj];
if ((segi[0] == segj[0] && segi[1] == segj[1]) || if ((segi[0] == segj[0] && segi[1] == segj[1]) ||
(segi[0] == segj[1] && segi[1] == segj[0])) { (segi[0] == segj[1] && segi[1] == segj[0])) {
@ -587,7 +623,8 @@ BitArray BoundaryLayerTool ::ProjectGrowthVectorsOnSurface() {
auto n = getNormal(sel); auto n = getNormal(sel);
for (auto i : Range(sel.PNums())) { for (auto i : Range(sel.PNums())) {
auto pi = sel.PNums()[i]; auto pi = sel.PNums()[i];
if (growthvectors[pi].Length2() == 0.) continue; if (growthvectors[pi].Length2() == 0.)
continue;
auto next = sel.PNums()[(i + 1) % sel.GetNV()]; auto next = sel.PNums()[(i + 1) % sel.GetNV()];
auto prev = sel.PNums()[i == 0 ? sel.GetNV() - 1 : i - 1]; auto prev = sel.PNums()[i == 0 ? sel.GetNV() - 1 : i - 1];
auto v1 = (mesh[next] - mesh[pi]).Normalize(); auto v1 = (mesh[next] - mesh[pi]).Normalize();
@ -600,7 +637,8 @@ BitArray BoundaryLayerTool ::ProjectGrowthVectorsOnSurface() {
else else
continue; continue;
if (!par_project_boundaries.Contains(sel.GetIndex())) continue; if (!par_project_boundaries.Contains(sel.GetIndex()))
continue;
auto &g = growthvectors[pi]; auto &g = growthvectors[pi];
auto ng = n * g; auto ng = n * g;
auto gg = g * g; auto gg = g * g;
@ -632,9 +670,11 @@ BitArray BoundaryLayerTool ::ProjectGrowthVectorsOnSurface() {
void BoundaryLayerTool ::InterpolateGrowthVectors() { void BoundaryLayerTool ::InterpolateGrowthVectors() {
int new_max_edge_nr = max_edge_nr; int new_max_edge_nr = max_edge_nr;
for (const auto &seg : segments) for (const auto &seg : segments)
if (seg.edgenr > new_max_edge_nr) new_max_edge_nr = seg.edgenr; if (seg.edgenr > new_max_edge_nr)
new_max_edge_nr = seg.edgenr;
for (const auto &seg : new_segments) for (const auto &seg : new_segments)
if (seg.edgenr > new_max_edge_nr) new_max_edge_nr = seg.edgenr; if (seg.edgenr > new_max_edge_nr)
new_max_edge_nr = seg.edgenr;
auto getGW = [&](PointIndex pi) -> Vec<3> { auto getGW = [&](PointIndex pi) -> Vec<3> {
if (growth_vector_map.count(pi) == 0) if (growth_vector_map.count(pi) == 0)
@ -664,10 +704,12 @@ void BoundaryLayerTool ::InterpolateGrowthVectors() {
// return true; // return true;
// return false; // return false;
auto segs = topo.GetVertexSegments(pi); auto segs = topo.GetVertexSegments(pi);
if (segs.Size() == 1) return true; if (segs.Size() == 1)
return true;
auto first_edgenr = mesh[segs[0]].edgenr; auto first_edgenr = mesh[segs[0]].edgenr;
for (auto segi : segs) for (auto segi : segs)
if (mesh[segi].edgenr != first_edgenr) return true; if (mesh[segi].edgenr != first_edgenr)
return true;
return false; return false;
}; };
@ -680,7 +722,8 @@ void BoundaryLayerTool ::InterpolateGrowthVectors() {
if (points.Size() == 0 && if (points.Size() == 0 &&
(is_end_point(seg[0]) || is_end_point(seg[1]))) { (is_end_point(seg[0]) || is_end_point(seg[1]))) {
PointIndex seg0 = seg[0], seg1 = seg[1]; PointIndex seg0 = seg[0], seg1 = seg[1];
if (is_end_point(seg[1])) Swap(seg0, seg1); if (is_end_point(seg[1]))
Swap(seg0, seg1);
points.Append(seg0); points.Append(seg0);
points.Append(seg1); points.Append(seg1);
edge_len += (mesh[seg[1]] - mesh[seg[0]]).Length(); edge_len += (mesh[seg[1]] - mesh[seg[0]]).Length();
@ -694,13 +737,15 @@ void BoundaryLayerTool ::InterpolateGrowthVectors() {
} }
if (!points.Size()) if (!points.Size())
throw Exception("Could not find startpoint for edge " + ToString(edgenr)); throw Exception("Could not find startpoint for edge " +
ToString(edgenr));
while (true) { while (true) {
bool point_found = false; bool point_found = false;
for (auto si : topo.GetVertexSegments(points.Last())) { for (auto si : topo.GetVertexSegments(points.Last())) {
const auto &seg = mesh[si]; const auto &seg = mesh[si];
if (seg.edgenr - 1 != edgenr) continue; if (seg.edgenr - 1 != edgenr)
continue;
if (seg[0] == points.Last() && points[points.Size() - 2] != seg[1]) { if (seg[0] == points.Last() && points[points.Size() - 2] != seg[1]) {
edge_len += (mesh[points.Last()] - mesh[seg[1]]).Length(); edge_len += (mesh[points.Last()] - mesh[seg[1]]).Length();
points.Append(seg[1]); points.Append(seg[1]);
@ -714,15 +759,18 @@ void BoundaryLayerTool ::InterpolateGrowthVectors() {
break; break;
} }
} }
if (is_end_point(points.Last())) break; if (is_end_point(points.Last()))
break;
if (!point_found) { if (!point_found) {
throw Exception( throw Exception(
string("Could not find connected list of line segments for edge ") + string(
"Could not find connected list of line segments for edge ") +
edgenr); edgenr);
} }
} }
if (getGW(points[0]).Length2() == 0 && getGW(points.Last()).Length2() == 0) if (getGW(points[0]).Length2() == 0 &&
getGW(points.Last()).Length2() == 0)
continue; continue;
// cout << "Points to average " << endl << points << endl; // cout << "Points to average " << endl << points << endl;
@ -767,7 +815,8 @@ void BoundaryLayerTool ::InsertNewElements(
mapfrom = PointIndex::INVALID; mapfrom = PointIndex::INVALID;
auto changed_domains = domains; auto changed_domains = domains;
if (!params.outside) changed_domains.Invert(); if (!params.outside)
changed_domains.Invert();
auto &identifications = mesh.GetIdentifications(); auto &identifications = mesh.GetIdentifications();
const int identnr = identifications.GetNr("boundarylayer"); const int identnr = identifications.GetNr("boundarylayer");
@ -783,7 +832,8 @@ void BoundaryLayerTool ::InsertNewElements(
mapfrom.Append(pi); mapfrom.Append(pi);
new_points.Append(pi_new); new_points.Append(pi_new);
growth_vector_map[pi_new] = {&growth_vector, height}; growth_vector_map[pi_new] = {&growth_vector, height};
if (special_boundary_points.count(pi) > 0) mesh.AddLockedPoint(pi_new); if (special_boundary_points.count(pi) > 0)
mesh.AddLockedPoint(pi_new);
pi_last = pi_new; pi_last = pi_new;
} }
}; };
@ -802,7 +852,8 @@ void BoundaryLayerTool ::InsertNewElements(
// get point from mapto (or the group if point is mapped to multiple new // get point from mapto (or the group if point is mapped to multiple new
// points) layer = -1 means last point (top of boundary layer) // points) layer = -1 means last point (top of boundary layer)
auto newPoint = [&](PointIndex pi, int layer = -1, int group = 0) { auto newPoint = [&](PointIndex pi, int layer = -1, int group = 0) {
if (layer == -1) layer = par_heights.Size() - 1; if (layer == -1)
layer = par_heights.Size() - 1;
if (special_boundary_points.count(pi)) if (special_boundary_points.count(pi))
return special_boundary_points[pi].growth_groups[group].new_points[layer]; return special_boundary_points[pi].growth_groups[group].new_points[layer];
else else
@ -829,7 +880,8 @@ void BoundaryLayerTool ::InsertNewElements(
} }
const auto &all_groups = special_boundary_points[pi].growth_groups; const auto &all_groups = special_boundary_points[pi].growth_groups;
for (auto i : Range(n)) for (auto i : Range(n))
if (all_groups[i].faces.Contains(face_index)) groups.Append(i); if (all_groups[i].faces.Contains(face_index))
groups.Append(i);
// cout << "groups " << pi << ", " << face_index << endl << groups; // cout << "groups " << pi << ", " << face_index << endl << groups;
return groups; return groups;
}; };
@ -839,7 +891,8 @@ void BoundaryLayerTool ::InsertNewElements(
map<int, int> edge_map; map<int, int> edge_map;
int edge_nr = max_edge_nr; int edge_nr = max_edge_nr;
auto getEdgeNr = [&](int ei) { auto getEdgeNr = [&](int ei) {
if (edge_map.count(ei) == 0) edge_map[ei] = ++edge_nr; if (edge_map.count(ei) == 0)
edge_map[ei] = ++edge_nr;
return edge_map[ei]; return edge_map[ei];
}; };
if (params.grow_edges) { if (params.grow_edges) {
@ -951,10 +1004,12 @@ void BoundaryLayerTool ::InsertNewElements(
auto getClosestGroup = [&](PointIndex pi, SurfaceElementIndex sei) { auto getClosestGroup = [&](PointIndex pi, SurfaceElementIndex sei) {
auto n = numGroups(pi); auto n = numGroups(pi);
if (n == 1) return 0; if (n == 1)
return 0;
const auto &sel = mesh[sei]; const auto &sel = mesh[sei];
auto groups = getGroups(pi, sel.GetIndex()); auto groups = getGroups(pi, sel.GetIndex());
if (groups.Size() == 1) return groups[0]; if (groups.Size() == 1)
return groups[0];
auto &growth_groups = special_boundary_points[pi].growth_groups; auto &growth_groups = special_boundary_points[pi].growth_groups;
@ -974,20 +1029,26 @@ void BoundaryLayerTool ::InsertNewElements(
const auto sel = mesh[si]; const auto sel = mesh[si];
if (moved_surfaces.Test(sel.GetIndex())) { if (moved_surfaces.Test(sel.GetIndex())) {
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]);
ArrayMem<int, 4> groups(points.Size()); ArrayMem<int, 4> groups(points.Size());
for (auto i : Range(points)) groups[i] = getClosestGroup(sel[i], si); for (auto i : Range(points))
groups[i] = getClosestGroup(sel[i], si);
bool add_volume_element = true; bool add_volume_element = true;
for (auto pi : sel.PNums()) for (auto pi : sel.PNums())
if (numGroups(pi) > 1) add_volume_element = false; if (numGroups(pi) > 1)
add_volume_element = false;
for (auto j : Range(par_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);
for (auto i : Range(points)) el[i] = points[i]; for (auto i : Range(points))
el[i] = points[i];
for (auto i : Range(points)) for (auto i : Range(points))
points[i] = newPoint(sel.PNums()[i], j, groups[i]); points[i] = newPoint(sel.PNums()[i], j, groups[i]);
if (surfacefacs[sel.GetIndex()] > 0) Swap(points[0], points[2]); if (surfacefacs[sel.GetIndex()] > 0)
for (auto i : Range(points)) el[sel.PNums().Size() + i] = points[i]; Swap(points[0], points[2]);
for (auto i : Range(points))
el[sel.PNums().Size() + i] = points[i];
auto new_index = new_mat_nrs[sel.GetIndex()]; auto new_index = new_mat_nrs[sel.GetIndex()];
if (new_index == -1) if (new_index == -1)
throw Exception("Boundary " + ToString(sel.GetIndex()) + throw Exception("Boundary " + ToString(sel.GetIndex()) +
@ -1006,7 +1067,8 @@ void BoundaryLayerTool ::InsertNewElements(
} }
} }
Element2d newel = sel; Element2d newel = sel;
for (auto i : Range(points)) newel[i] = newPoint(sel[i], -1, groups[i]); for (auto i : Range(points))
newel[i] = newPoint(sel[i], -1, groups[i]);
newel.SetIndex(si_map[sel.GetIndex()]); newel.SetIndex(si_map[sel.GetIndex()]);
mesh.AddSurfaceElement(newel); mesh.AddSurfaceElement(newel);
@ -1036,7 +1098,8 @@ void BoundaryLayerTool ::InsertNewElements(
} }
} else { } else {
bool has_moved = false; bool has_moved = false;
for (auto p : sel.PNums()) has_moved |= hasMoved(p); for (auto p : sel.PNums())
has_moved |= hasMoved(p);
if (has_moved) if (has_moved)
for (auto p : sel.PNums()) { for (auto p : sel.PNums()) {
if (hasMoved(p)) { if (hasMoved(p)) {
@ -1048,7 +1111,8 @@ void BoundaryLayerTool ::InsertNewElements(
} }
if (is_boundary_moved.Test(sel.GetIndex())) { if (is_boundary_moved.Test(sel.GetIndex())) {
for (auto &p : mesh[si].PNums()) for (auto &p : mesh[si].PNums())
if (hasMoved(p)) p = newPoint(p); if (hasMoved(p))
p = newPoint(p);
} }
} }
@ -1056,7 +1120,8 @@ void BoundaryLayerTool ::InsertNewElements(
auto &seg = segments[sei]; auto &seg = segments[sei];
if (is_boundary_moved.Test(seg.si)) if (is_boundary_moved.Test(seg.si))
for (auto &p : seg.PNums()) for (auto &p : seg.PNums())
if (hasMoved(p)) p = newPoint(p); if (hasMoved(p))
p = newPoint(p);
// else if(hasMoved(seg[0]) || hasMoved(seg[1])) // else if(hasMoved(seg[0]) || hasMoved(seg[1]))
// { // {
// auto tangent = mesh[seg[1]] - mesh[seg[0]]; // auto tangent = mesh[seg[1]] - mesh[seg[0]];
@ -1067,25 +1132,28 @@ void BoundaryLayerTool ::InsertNewElements(
// } // }
} }
// fill holes in surface mesh at special boundary points (i.e. points with >=4 adjacent // fill holes in surface mesh at special boundary points (i.e. points with >=4
// boundary faces) // adjacent boundary faces)
auto p2sel = mesh.CreatePoint2SurfaceElementTable(); auto p2sel = mesh.CreatePoint2SurfaceElementTable();
for (auto &[special_pi, special_point] : special_boundary_points) { for (auto &[special_pi, special_point] : special_boundary_points) {
if (special_point.growth_groups.Size() != 2) if (special_point.growth_groups.Size() != 2)
throw Exception("special_point.growth_groups.Size() != 2"); throw Exception("special_point.growth_groups.Size() != 2");
// Special points are split into two new points, when mapping a surface element, we choose the closer one to the center. // Special points are split into two new points, when mapping a surface
// Now, find points which are mapped to both new points (for different surface elements they belong to). // element, we choose the closer one to the center. Now, find points which
// At exactly these points we need to insert new surface elements to fill the hole. // are mapped to both new points (for different surface elements they belong
// to). At exactly these points we need to insert new surface elements to
// fill the hole.
std::map<int, std::array<std::set<PointIndex>, 2>> close_group; std::map<int, std::array<std::set<PointIndex>, 2>> close_group;
for (auto sei : p2sel[special_pi]) { for (auto sei : p2sel[special_pi]) {
const auto &sel = mesh[sei]; const auto &sel = mesh[sei];
for (auto p : sel.PNums()) for (auto p : sel.PNums())
if (p != special_pi) close_group[sel.GetIndex()][getClosestGroup(special_pi, sei)].insert(p); if (p != special_pi)
close_group[sel.GetIndex()][getClosestGroup(special_pi, sei)].insert(
p);
} }
for( auto [fi, groups] : close_group ) for (auto [fi, groups] : close_group) {
{
const auto mapped_fi = si_map[fi]; const auto mapped_fi = si_map[fi];
std::set<PointIndex> common_points; std::set<PointIndex> common_points;
for (auto pi : groups[0]) for (auto pi : groups[0])
@ -1096,7 +1164,8 @@ void BoundaryLayerTool ::InsertNewElements(
auto new_special_pi0 = special_point.growth_groups[0].new_points.Last(); auto new_special_pi0 = special_point.growth_groups[0].new_points.Last();
auto new_special_pi1 = special_point.growth_groups[1].new_points.Last(); auto new_special_pi1 = special_point.growth_groups[1].new_points.Last();
for (auto sei : p2sel[pi_common]) { for (auto sei : p2sel[pi_common]) {
if(mesh[sei].GetIndex() == mapped_fi && mesh[sei].PNums().Contains(new_special_pi0)) { if (mesh[sei].GetIndex() == mapped_fi &&
mesh[sei].PNums().Contains(new_special_pi0)) {
auto sel = mesh[sei]; auto sel = mesh[sei];
sel.Invert(); sel.Invert();
for (auto &pi : sel.PNums()) for (auto &pi : sel.PNums())
@ -1109,19 +1178,19 @@ void BoundaryLayerTool ::InsertNewElements(
} }
} }
for (auto &[pi, special_point] : special_boundary_points) { for (auto &[pi, special_point] : special_boundary_points) {
if (special_point.growth_groups.Size() != 2) if (special_point.growth_groups.Size() != 2)
throw Exception("special_point.growth_groups.Size() != 2"); throw Exception("special_point.growth_groups.Size() != 2");
for (auto igroup : Range(2)) { for (auto igroup : Range(2)) {
auto &group = special_point.growth_groups[igroup]; auto &group = special_point.growth_groups[igroup];
std::set<int> faces; std::set<int> faces;
for (auto face : group.faces) faces.insert(si_map[face]); for (auto face : group.faces)
faces.insert(si_map[face]);
auto pi_new = group.new_points.Last(); auto pi_new = group.new_points.Last();
auto pi_new_other = auto pi_new_other =
special_point.growth_groups[1 - igroup].new_points.Last(); special_point.growth_groups[1 - igroup].new_points.Last();
for (auto sei : p2sel[pi_new]) faces.erase(mesh[sei].GetIndex()); for (auto sei : p2sel[pi_new])
faces.erase(mesh[sei].GetIndex());
for (auto face : faces) for (auto face : faces)
for (auto seg : new_segments) { for (auto seg : new_segments) {
if ( // seg.si == face if ( // seg.si == face
@ -1154,9 +1223,12 @@ void BoundaryLayerTool ::InsertNewElements(
ArrayMem<PointIndex, 4> moved; ArrayMem<PointIndex, 4> moved;
bool moved_bnd = false; bool moved_bnd = false;
for (const auto &p : el.PNums()) { for (const auto &p : el.PNums()) {
if (fixed_points.Test(p)) fixed.Append(p); if (fixed_points.Test(p))
if (hasMoved(p)) moved.Append(p); fixed.Append(p);
if (moveboundarypoint.Test(p)) moved_bnd = true; if (hasMoved(p))
moved.Append(p);
if (moveboundarypoint.Test(p))
moved_bnd = true;
} }
bool do_move, do_insert; bool do_move, do_insert;
@ -1189,7 +1261,8 @@ void BoundaryLayerTool ::InsertNewElements(
auto v1 = mesh[p1]; auto v1 = mesh[p1];
auto n = Cross(mesh[p2] - v1, mesh[p3] - v1); auto n = Cross(mesh[p2] - v1, mesh[p3] - v1);
auto d = mesh[newPoint(p1, 0)] - v1; auto d = mesh[newPoint(p1, 0)] - v1;
if (n * d > 0) Swap(p2, p3); if (n * d > 0)
Swap(p2, p3);
PointIndex p4 = p1; PointIndex p4 = p1;
PointIndex p5 = p2; PointIndex p5 = p2;
PointIndex p6 = p3; PointIndex p6 = p3;
@ -1303,10 +1376,12 @@ void BoundaryLayerTool ::SetDomInOutSides() {
for (auto sei : Range(mesh.SurfaceElements())) { for (auto sei : Range(mesh.SurfaceElements())) {
auto &sel = mesh[sei]; auto &sel = mesh[sei];
auto index = sel.GetIndex(); auto index = sel.GetIndex();
if (done.Test(index)) continue; if (done.Test(index))
continue;
done.SetBit(index); done.SetBit(index);
auto &fd = mesh.GetFaceDescriptor(index); auto &fd = mesh.GetFaceDescriptor(index);
if (fd.DomainIn() != -1) continue; if (fd.DomainIn() != -1)
continue;
int e1, e2; int e1, e2;
mesh.GetTopology().GetSurface2VolumeElement(sei + 1, e1, e2); mesh.GetTopology().GetSurface2VolumeElement(sei + 1, e1, e2);
if (e1 == 0) if (e1 == 0)
@ -1325,7 +1400,8 @@ void BoundaryLayerTool ::AddSegments() {
MergeAndAddSegments(mesh, segments, new_segments); MergeAndAddSegments(mesh, segments, new_segments);
else { else {
mesh.LineSegments() = segments; mesh.LineSegments() = segments;
for (auto& seg : new_segments) mesh.AddSegment(seg); for (auto &seg : new_segments)
mesh.AddSegment(seg);
} }
} }
@ -1336,16 +1412,19 @@ void BoundaryLayerTool ::FixVolumeElements() {
is_inner_point.Clear(); is_inner_point.Clear();
auto changed_domains = domains; auto changed_domains = domains;
if (!params.outside) changed_domains.Invert(); if (!params.outside)
changed_domains.Invert();
for (ElementIndex ei : Range(ne)) for (ElementIndex ei : Range(ne))
if (changed_domains.Test(mesh[ei].GetIndex())) if (changed_domains.Test(mesh[ei].GetIndex()))
for (auto pi : mesh[ei].PNums()) for (auto pi : mesh[ei].PNums())
if (mesh[pi].Type() == INNERPOINT) is_inner_point.SetBit(pi); if (mesh[pi].Type() == INNERPOINT)
is_inner_point.SetBit(pi);
Array<PointIndex> points; Array<PointIndex> points;
for (auto pi : mesh.Points().Range()) for (auto pi : mesh.Points().Range())
if (is_inner_point.Test(pi)) points.Append(pi); if (is_inner_point.Test(pi))
points.Append(pi);
auto p2el = mesh.CreatePoint2ElementTable(is_inner_point); auto p2el = mesh.CreatePoint2ElementTable(is_inner_point);
@ -1368,36 +1447,36 @@ void BoundaryLayerTool ::FixVolumeElements() {
} }
} }
void BoundaryLayerTool :: ProcessParameters() void BoundaryLayerTool ::ProcessParameters() {
{ if (int *bc = get_if<int>(&params.boundary); bc) {
if(int* bc = get_if<int>(&params.boundary); bc)
{
for (int i = 1; i <= mesh.GetNFD(); i++) for (int i = 1; i <= mesh.GetNFD(); i++)
if (mesh.GetFaceDescriptor(i).BCProperty() == *bc) if (mesh.GetFaceDescriptor(i).BCProperty() == *bc)
par_surfid.Append(i); par_surfid.Append(i);
} } else if (string *s = get_if<string>(&params.boundary); s) {
else if(string* s = get_if<string>(&params.boundary); s)
{
regex pattern(*s); regex pattern(*s);
BitArray boundaries(mesh.GetNFD() + 1); BitArray boundaries(mesh.GetNFD() + 1);
boundaries.Clear(); boundaries.Clear();
for(int i = 1; i<=mesh.GetNFD(); i++) for (int i = 1; i <= mesh.GetNFD(); i++) {
{
auto &fd = mesh.GetFaceDescriptor(i); auto &fd = mesh.GetFaceDescriptor(i);
if(regex_match(fd.GetBCName(), pattern)) if (regex_match(fd.GetBCName(), pattern)) {
{
boundaries.SetBit(i); boundaries.SetBit(i);
auto dom_pattern = get_if<string>(&params.domain); auto dom_pattern = get_if<string>(&params.domain);
// only add if adjacent to domain // only add if adjacent to domain
if(dom_pattern) if (dom_pattern) {
{
regex pattern(*dom_pattern); regex pattern(*dom_pattern);
bool mat1_match = fd.DomainIn() > 0 && regex_match(mesh.GetMaterial(fd.DomainIn()), pattern); bool mat1_match =
bool mat2_match = fd.DomainOut() > 0 && regex_match(mesh.GetMaterial(fd.DomainOut()), pattern); 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 boundary is inner or outer remove from list
if (mat1_match == mat2_match) if (mat1_match == mat2_match)
boundaries.Clear(i); boundaries.Clear(i);
// if((fd.DomainIn() > 0 && regex_match(mesh.GetMaterial(fd.DomainIn()), pattern)) || (fd.DomainOut() > 0 && regex_match(self.GetMaterial(fd.DomainOut()), pattern))) // if((fd.DomainIn() > 0 &&
// regex_match(mesh.GetMaterial(fd.DomainIn()), pattern)) ||
// (fd.DomainOut() > 0 &&
// regex_match(self.GetMaterial(fd.DomainOut()), pattern)))
// boundaries.Clear(i); // boundaries.Clear(i);
// par_surfid.Append(i); // par_surfid.Append(i);
} }
@ -1408,9 +1487,7 @@ void BoundaryLayerTool ::FixVolumeElements() {
for (int i = 1; i <= mesh.GetNFD(); i++) for (int i = 1; i <= mesh.GetNFD(); i++)
if (boundaries.Test(i)) if (boundaries.Test(i))
par_surfid.Append(i); par_surfid.Append(i);
} } else {
else
{
auto &surfids = *get_if<std::vector<int>>(&params.boundary); auto &surfids = *get_if<std::vector<int>>(&params.boundary);
for (auto id : surfids) for (auto id : surfids)
par_surfid.Append(id); par_surfid.Append(id);
@ -1420,29 +1497,22 @@ void BoundaryLayerTool ::FixVolumeElements() {
else else
par_new_mat = *get_if<map<string, string>>(&params.new_material); par_new_mat = *get_if<map<string, string>>(&params.new_material);
if(params.project_boundaries.has_value()) if (params.project_boundaries.has_value()) {
{
auto proj_bnd = *params.project_boundaries; auto proj_bnd = *params.project_boundaries;
if(string* s = get_if<string>(&proj_bnd); s) if (string *s = get_if<string>(&proj_bnd); s) {
{
regex pattern(*s); regex pattern(*s);
for (int i = 1; i <= mesh.GetNFD(); i++) for (int i = 1; i <= mesh.GetNFD(); i++)
if (regex_match(mesh.GetFaceDescriptor(i).GetBCName(), pattern)) if (regex_match(mesh.GetFaceDescriptor(i).GetBCName(), pattern))
par_project_boundaries.Append(i); par_project_boundaries.Append(i);
} } else {
else
{
for (auto id : *get_if<std::vector<int>>(&proj_bnd)) for (auto id : *get_if<std::vector<int>>(&proj_bnd))
par_project_boundaries.Append(id); par_project_boundaries.Append(id);
} }
} }
if(double* height = get_if<double>(&params.thickness); height) if (double *height = get_if<double>(&params.thickness); height) {
{
par_heights.Append(*height); par_heights.Append(*height);
} } else {
else
{
auto &heights = *get_if<std::vector<double>>(&params.thickness); auto &heights = *get_if<std::vector<double>>(&params.thickness);
for (auto val : heights) for (auto val : heights)
par_heights.Append(val); par_heights.Append(val);
@ -1451,19 +1521,14 @@ void BoundaryLayerTool ::FixVolumeElements() {
int nr_domains = mesh.GetNDomains(); int nr_domains = mesh.GetNDomains();
domains.SetSize(nr_domains + 1); // one based domains.SetSize(nr_domains + 1); // one based
domains.Clear(); domains.Clear();
if(string* pdomain = get_if<string>(&params.domain); pdomain) if (string *pdomain = get_if<string>(&params.domain); pdomain) {
{
regex pattern(*pdomain); regex pattern(*pdomain);
for (auto i : Range(1, nr_domains + 1)) for (auto i : Range(1, nr_domains + 1))
if (regex_match(mesh.GetMaterial(i), pattern)) if (regex_match(mesh.GetMaterial(i), pattern))
domains.SetBit(i); domains.SetBit(i);
} } else if (int *idomain = get_if<int>(&params.domain); idomain) {
else if(int *idomain = get_if<int>(&params.domain); idomain)
{
domains.SetBit(*idomain); domains.SetBit(*idomain);
} } else {
else
{
for (auto i : *get_if<std::vector<int>>(&params.domain)) for (auto i : *get_if<std::vector<int>>(&params.domain))
domains.SetBit(i); domains.SetBit(i);
} }
@ -1487,7 +1552,8 @@ void BoundaryLayerTool ::Perform() {
mesh.UpdateTopology(); mesh.UpdateTopology();
InterpolateGrowthVectors(); InterpolateGrowthVectors();
if (params.limit_growth_vectors) LimitGrowthVectorLengths(); if (params.limit_growth_vectors)
LimitGrowthVectorLengths();
for (auto [pi, data] : growth_vector_map) { for (auto [pi, data] : growth_vector_map) {
auto [gw, height] = data; auto [gw, height] = data;