2024-04-05 20:09:34 +05:00
|
|
|
#include "boundarylayer.hpp"
|
2024-09-27 20:03:26 +05:00
|
|
|
#include "boundarylayer_limiter.hpp"
|
2009-01-13 04:40:13 +05:00
|
|
|
|
2022-11-10 18:35:58 +05:00
|
|
|
#include <regex>
|
2024-04-05 20:09:34 +05:00
|
|
|
#include <set>
|
2023-09-04 16:43:47 +05:00
|
|
|
|
|
|
|
#include "debugging.hpp"
|
2024-04-05 20:09:34 +05:00
|
|
|
#include "global.hpp"
|
2023-09-04 16:43:47 +05:00
|
|
|
#include "meshfunc.hpp"
|
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
namespace netgen {
|
2022-03-04 19:42:16 +05:00
|
|
|
|
2024-09-04 18:49:02 +05:00
|
|
|
struct SpecialPointException : public Exception {
|
|
|
|
SpecialPointException() : Exception("") {}
|
|
|
|
};
|
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
std::tuple<int, int> FindCloseVectors(FlatArray<Vec<3>> ns,
|
|
|
|
bool find_max = true) {
|
2023-11-29 13:19:27 +05:00
|
|
|
int maxpos1;
|
|
|
|
int maxpos2;
|
2024-04-05 20:09:34 +05:00
|
|
|
|
2023-11-29 13:19:27 +05:00
|
|
|
double val = find_max ? -1e99 : 1e99;
|
|
|
|
for (auto i : Range(ns))
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto j : Range(i + 1, ns.Size())) {
|
|
|
|
double ip = ns[i] * ns[j];
|
|
|
|
if ((find_max && (ip > val)) || (!find_max && (ip < val))) {
|
|
|
|
val = ip;
|
|
|
|
maxpos1 = i;
|
|
|
|
maxpos2 = j;
|
|
|
|
}
|
|
|
|
}
|
2023-11-29 13:19:27 +05:00
|
|
|
return {maxpos1, maxpos2};
|
|
|
|
}
|
|
|
|
|
|
|
|
Vec<3> CalcGrowthVector(FlatArray<Vec<3>> ns) {
|
2024-09-27 20:30:36 +05:00
|
|
|
if (ns.Size() == 0)
|
|
|
|
return {0, 0, 0};
|
|
|
|
if (ns.Size() == 1)
|
|
|
|
return ns[0];
|
2024-04-05 20:09:34 +05:00
|
|
|
if (ns.Size() == 2) {
|
|
|
|
auto gw = ns[0];
|
|
|
|
auto n = ns[1];
|
|
|
|
auto npn = gw * n;
|
|
|
|
auto npnp = gw * gw;
|
|
|
|
auto nn = n * n;
|
2024-09-27 20:30:36 +05:00
|
|
|
if (fabs(nn - npn * npn / npnp) < 1e-6)
|
|
|
|
return n;
|
2024-04-05 20:09:34 +05:00
|
|
|
gw += (nn - npn) / (nn - npn * npn / npnp) * (n - npn / npnp * gw);
|
|
|
|
return gw;
|
|
|
|
}
|
|
|
|
if (ns.Size() == 3) {
|
|
|
|
DenseMatrix mat(3, 3);
|
|
|
|
for (auto i : Range(3))
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto j : Range(3))
|
|
|
|
mat(i, j) = ns[i][j];
|
2024-04-05 20:09:34 +05:00
|
|
|
|
|
|
|
if (fabs(mat.Det()) > 1e-6) {
|
|
|
|
DenseMatrix mat(3, 3);
|
|
|
|
for (auto i : Range(3))
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto j : Range(3))
|
|
|
|
mat(i, j) = ns[i] * ns[j];
|
2024-04-05 20:09:34 +05:00
|
|
|
Vector rhs(3);
|
|
|
|
rhs = 1.;
|
|
|
|
Vector res(3);
|
|
|
|
DenseMatrix inv(3, ns.Size());
|
|
|
|
CalcInverse(mat, inv);
|
|
|
|
inv.Mult(rhs, res);
|
|
|
|
Vec<3> growth = 0.;
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto i : Range(ns))
|
|
|
|
growth += res[i] * ns[i];
|
2024-04-05 20:09:34 +05:00
|
|
|
return growth;
|
2023-11-29 13:19:27 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2023-11-29 13:19:27 +05:00
|
|
|
auto [maxpos1, maxpos2] = FindCloseVectors(ns);
|
|
|
|
Array<Vec<3>> new_normals;
|
|
|
|
new_normals = ns;
|
|
|
|
const auto dot = ns[maxpos1] * ns[maxpos2];
|
2024-04-05 20:09:34 +05:00
|
|
|
auto average = 0.5 * (ns[maxpos1] + ns[maxpos2]);
|
2023-11-29 13:19:27 +05:00
|
|
|
average.Normalize();
|
|
|
|
new_normals[maxpos1] = average;
|
|
|
|
new_normals.DeleteElement(maxpos2);
|
|
|
|
auto gw = CalcGrowthVector(new_normals);
|
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto n : ns)
|
|
|
|
if (n * gw < 0)
|
2024-09-04 18:49:02 +05:00
|
|
|
throw SpecialPointException();
|
2023-11-29 13:19:27 +05:00
|
|
|
return gw;
|
|
|
|
}
|
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
SpecialBoundaryPoint ::GrowthGroup ::GrowthGroup(FlatArray<int> faces_,
|
|
|
|
FlatArray<Vec<3>> normals) {
|
2023-11-29 13:19:27 +05:00
|
|
|
faces = faces_;
|
|
|
|
growth_vector = CalcGrowthVector(normals);
|
|
|
|
}
|
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
SpecialBoundaryPoint ::SpecialBoundaryPoint(
|
2024-09-27 20:30:36 +05:00
|
|
|
const std::map<int, Vec<3>> &normals) {
|
2023-11-29 13:19:27 +05:00
|
|
|
// find opposing face normals
|
|
|
|
Array<Vec<3>> ns;
|
|
|
|
Array<int> faces;
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto [face, normal] : normals) {
|
2023-11-29 13:19:27 +05:00
|
|
|
ns.Append(normal);
|
|
|
|
faces.Append(face);
|
|
|
|
}
|
|
|
|
|
|
|
|
auto [minface1, minface2] = FindCloseVectors(ns, false);
|
|
|
|
minface1 = faces[minface1];
|
|
|
|
minface2 = faces[minface2];
|
|
|
|
Array<int> g1_faces;
|
|
|
|
g1_faces.Append(minface1);
|
|
|
|
Array<int> g2_faces;
|
|
|
|
g2_faces.Append(minface2);
|
2024-09-04 18:49:02 +05:00
|
|
|
auto n1 = normals.at(minface1);
|
|
|
|
auto n2 = normals.at(minface2);
|
2024-09-27 20:30:36 +05:00
|
|
|
separating_direction = 0.5 * (n2 - n1);
|
2023-11-29 13:19:27 +05:00
|
|
|
|
|
|
|
Array<Vec<3>> normals1, normals2;
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto [facei, normali] : normals)
|
|
|
|
if (facei != minface1 && facei != minface2) {
|
2023-11-29 13:19:27 +05:00
|
|
|
g1_faces.Append(facei);
|
|
|
|
g2_faces.Append(facei);
|
|
|
|
}
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto fi : g1_faces)
|
|
|
|
normals1.Append(normals.at(fi));
|
|
|
|
for (auto fi : g2_faces)
|
|
|
|
normals2.Append(normals.at(fi));
|
2023-11-29 13:19:27 +05:00
|
|
|
growth_groups.Append(GrowthGroup(g1_faces, normals1));
|
|
|
|
growth_groups.Append(GrowthGroup(g2_faces, normals2));
|
|
|
|
}
|
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
Vec<3> BoundaryLayerTool ::getEdgeTangent(PointIndex pi, int edgenr) {
|
|
|
|
Vec<3> tangent = 0.0;
|
|
|
|
ArrayMem<PointIndex, 2> pts;
|
|
|
|
for (auto segi : topo.GetVertexSegments(pi)) {
|
2024-09-27 20:30:36 +05:00
|
|
|
auto &seg = mesh[segi];
|
|
|
|
if (seg.edgenr != edgenr + 1)
|
|
|
|
continue;
|
2024-04-05 20:09:34 +05:00
|
|
|
PointIndex other = seg[0] + seg[1] - pi;
|
2024-09-27 20:30:36 +05:00
|
|
|
if (!pts.Contains(other))
|
|
|
|
pts.Append(other);
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
|
|
|
if (pts.Size() != 2) {
|
|
|
|
cout << "getEdgeTangent pi = " << pi << ", edgenr = " << edgenr << endl;
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto segi : topo.GetVertexSegments(pi))
|
|
|
|
cout << mesh[segi] << endl;
|
2024-04-05 20:09:34 +05:00
|
|
|
throw Exception("Something went wrong in getEdgeTangent!");
|
|
|
|
}
|
|
|
|
tangent = mesh[pts[1]] - mesh[pts[0]];
|
|
|
|
return tangent.Normalize();
|
|
|
|
}
|
2023-11-22 23:21:26 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
void BoundaryLayerTool ::LimitGrowthVectorLengths() {
|
|
|
|
static Timer tall("BoundaryLayerTool::LimitGrowthVectorLengths");
|
|
|
|
RegionTimer rtall(tall);
|
2022-03-07 19:58:09 +05:00
|
|
|
|
2024-09-27 20:03:26 +05:00
|
|
|
GrowthVectorLimiter limiter(*this);
|
|
|
|
limiter.Perform();
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2022-03-07 19:58:09 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
// depending on the geometry type, the mesh contains segments multiple times
|
|
|
|
// (once for each face)
|
2024-09-27 20:30:36 +05:00
|
|
|
bool HaveSingleSegments(const Mesh &mesh) {
|
|
|
|
auto &topo = mesh.GetTopology();
|
2024-04-05 20:09:34 +05:00
|
|
|
NgArray<SurfaceElementIndex> surf_els;
|
2022-02-16 23:51:54 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto segi : Range(mesh.LineSegments())) {
|
|
|
|
mesh.GetTopology().GetSegmentSurfaceElements(segi + 1, surf_els);
|
2024-09-27 20:30:36 +05:00
|
|
|
if (surf_els.Size() < 2)
|
|
|
|
continue;
|
2022-02-16 23:51:54 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
auto seg = mesh[segi];
|
|
|
|
auto pi0 = min(seg[0], seg[1]);
|
|
|
|
auto pi1 = max(seg[0], seg[1]);
|
|
|
|
auto p0_segs = topo.GetVertexSegments(seg[0]);
|
2022-02-16 23:51:54 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto segi_other : p0_segs) {
|
2024-09-27 20:30:36 +05:00
|
|
|
if (segi_other == segi)
|
|
|
|
continue;
|
2022-02-16 23:51:54 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
auto seg_other = mesh[segi_other];
|
|
|
|
auto pi0_other = min(seg_other[0], seg_other[1]);
|
|
|
|
auto pi1_other = max(seg_other[0], seg_other[1]);
|
2024-09-27 20:30:36 +05:00
|
|
|
if (pi0_other == pi0 && pi1_other == pi1)
|
|
|
|
return false;
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
|
|
|
|
|
|
|
// found segment with multiple adjacent surface elements but no other
|
|
|
|
// segments with same points -> have single segments
|
|
|
|
return true;
|
2022-02-16 23:51:54 +05:00
|
|
|
}
|
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
return true;
|
|
|
|
}
|
2022-02-16 23:51:54 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
// duplicates segments (and sets seg.si accordingly) to have a unified data
|
|
|
|
// structure for all geometry types
|
2024-09-27 20:30:36 +05:00
|
|
|
Array<Segment> BuildSegments(Mesh &mesh) {
|
2024-04-05 20:09:34 +05:00
|
|
|
Array<Segment> segments;
|
|
|
|
// auto& topo = mesh.GetTopology();
|
2023-04-27 18:35:10 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
NgArray<SurfaceElementIndex> surf_els;
|
2024-02-28 22:08:08 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto segi : Range(mesh.LineSegments())) {
|
|
|
|
auto seg = mesh[segi];
|
|
|
|
mesh.GetTopology().GetSegmentSurfaceElements(segi + 1, surf_els);
|
|
|
|
for (auto seli : surf_els) {
|
2024-09-27 20:30:36 +05:00
|
|
|
const auto &sel = mesh[seli];
|
2024-04-05 20:09:34 +05:00
|
|
|
seg.si = sel.GetIndex();
|
2024-02-28 22:08:08 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
auto np = sel.GetNP();
|
|
|
|
for (auto i : Range(np)) {
|
|
|
|
if (sel[i] == seg[0]) {
|
2024-09-27 20:30:36 +05:00
|
|
|
if (sel[(i + 1) % np] != seg[1])
|
|
|
|
swap(seg[0], seg[1]);
|
2024-04-05 20:09:34 +05:00
|
|
|
break;
|
|
|
|
}
|
|
|
|
}
|
2022-02-16 23:51:54 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
segments.Append(seg);
|
2023-11-29 13:19:27 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
|
|
|
return segments;
|
|
|
|
}
|
2022-03-04 19:42:16 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
void MergeAndAddSegments(Mesh &mesh, FlatArray<Segment> segments,
|
2024-04-05 20:09:34 +05:00
|
|
|
FlatArray<Segment> new_segments) {
|
|
|
|
INDEX_2_HASHTABLE<bool> already_added(segments.Size() +
|
|
|
|
2 * new_segments.Size());
|
2024-02-28 22:08:08 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
mesh.LineSegments().SetSize0();
|
2024-02-28 22:08:08 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
auto addSegment = [&](const auto &seg) {
|
2024-04-05 20:09:34 +05:00
|
|
|
INDEX_2 i2(seg[0], seg[1]);
|
|
|
|
i2.Sort();
|
|
|
|
if (!already_added.Used(i2)) {
|
|
|
|
mesh.AddSegment(seg);
|
|
|
|
already_added.Set(i2, true);
|
|
|
|
}
|
|
|
|
};
|
2022-02-28 21:24:44 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
for (const auto &seg : segments)
|
|
|
|
addSegment(seg);
|
2024-04-05 20:09:34 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
for (const auto &seg : new_segments)
|
|
|
|
addSegment(seg);
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2024-02-28 22:08:08 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
// TODO: Hack, move this to the header or restructure the whole growth_vectors
|
|
|
|
// storage
|
2024-04-10 20:34:39 +05:00
|
|
|
static std::map<PointIndex, Vec<3>> non_bl_growth_vectors;
|
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
|
|
|
|
static Timer tall("InterpolateSurfaceGrowthVectors");
|
|
|
|
RegionTimer rtall(tall);
|
|
|
|
static Timer tsmooth("InterpolateSurfaceGrowthVectors-Smoothing");
|
|
|
|
auto np_old = this->np;
|
|
|
|
auto np = mesh.GetNP();
|
2024-04-10 20:34:39 +05:00
|
|
|
|
|
|
|
non_bl_growth_vectors.clear();
|
2024-04-05 20:09:34 +05:00
|
|
|
|
|
|
|
auto getGW = [&](PointIndex pi) -> Vec<3> {
|
2024-04-10 20:34:39 +05:00
|
|
|
if (growth_vector_map.count(pi) == 0) {
|
|
|
|
non_bl_growth_vectors[pi] = .0;
|
|
|
|
growth_vector_map[pi] = {&non_bl_growth_vectors[pi], 1.0};
|
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
auto [gw, height] = growth_vector_map[pi];
|
|
|
|
return height * (*gw);
|
|
|
|
};
|
2024-04-10 19:04:16 +05:00
|
|
|
auto addGW = [&](PointIndex pi, Vec<3> vec) {
|
2024-04-10 20:34:39 +05:00
|
|
|
if (growth_vector_map.count(pi) == 0) {
|
|
|
|
non_bl_growth_vectors[pi] = .0;
|
|
|
|
growth_vector_map[pi] = {&non_bl_growth_vectors[pi], 1.0};
|
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
auto [gw, height] = growth_vector_map[pi];
|
2024-04-10 19:04:16 +05:00
|
|
|
*gw += 1.0 / height * vec;
|
2024-04-05 20:09:34 +05:00
|
|
|
};
|
|
|
|
|
|
|
|
Array<Vec<3>, PointIndex> normals(np);
|
|
|
|
for (auto pi = np_old; pi < np; pi++) {
|
|
|
|
normals[pi + PointIndex::BASE] = getGW(pi + PointIndex::BASE);
|
2024-02-28 22:08:08 +05:00
|
|
|
}
|
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
auto hasMoved = [&](PointIndex pi) {
|
|
|
|
return (pi - PointIndex::BASE >= np_old) || mapto[pi].Size() > 0 ||
|
|
|
|
special_boundary_points.count(pi);
|
|
|
|
};
|
|
|
|
|
2024-04-10 19:04:16 +05:00
|
|
|
std::set<PointIndex> points_set;
|
2024-04-10 20:34:39 +05:00
|
|
|
ParallelForRange(mesh.SurfaceElements().Range(), [&](auto myrange) {
|
2024-04-05 20:09:34 +05:00
|
|
|
for (SurfaceElementIndex sei : myrange) {
|
2024-04-10 20:34:39 +05:00
|
|
|
for (auto pi : mesh[sei].PNums()) {
|
|
|
|
auto pi_from = mapfrom[pi];
|
2024-09-27 20:30:36 +05:00
|
|
|
if ((pi_from.IsValid() && mesh[pi_from].Type() == SURFACEPOINT) ||
|
|
|
|
(!pi_from.IsValid() && mapto[pi].Size() == 0 &&
|
|
|
|
mesh[pi].Type() == SURFACEPOINT))
|
2024-04-10 19:04:16 +05:00
|
|
|
points_set.insert(pi);
|
2024-04-10 20:34:39 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
|
|
|
});
|
2024-02-28 22:08:08 +05:00
|
|
|
|
2024-04-10 20:34:39 +05:00
|
|
|
Array<bool> has_moved_points(max_edge_nr + 1);
|
|
|
|
has_moved_points = false;
|
|
|
|
std::set<PointIndex> moved_edge_points;
|
|
|
|
|
|
|
|
for (auto seg : segments) {
|
|
|
|
if (hasMoved(seg[0]) != hasMoved(seg[1]))
|
|
|
|
has_moved_points[seg.edgenr] = true;
|
|
|
|
}
|
|
|
|
|
|
|
|
for (auto seg : segments)
|
|
|
|
if (has_moved_points[seg.edgenr])
|
|
|
|
for (auto pi : seg.PNums())
|
2024-09-27 20:30:36 +05:00
|
|
|
if (mesh[pi].Type() == EDGEPOINT)
|
|
|
|
points_set.insert(pi);
|
2024-04-10 20:34:39 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
Array<PointIndex> points;
|
2024-04-10 19:04:16 +05:00
|
|
|
for (auto pi : points_set)
|
|
|
|
points.Append(pi);
|
|
|
|
QuickSort(points);
|
2024-04-05 20:09:34 +05:00
|
|
|
|
|
|
|
auto p2sel = mesh.CreatePoint2SurfaceElementTable();
|
|
|
|
// smooth tangential part of growth vectors from edges to surface elements
|
2024-04-10 19:04:16 +05:00
|
|
|
Array<Vec<3>, PointIndex> corrections(mesh.GetNP());
|
|
|
|
corrections = 0.0;
|
2024-04-05 20:09:34 +05:00
|
|
|
RegionTimer rtsmooth(tsmooth);
|
|
|
|
for ([[maybe_unused]] auto i : Range(10)) {
|
|
|
|
for (auto pi : points) {
|
|
|
|
auto sels = p2sel[pi];
|
2024-09-27 20:30:36 +05:00
|
|
|
auto &correction = corrections[pi];
|
2024-04-05 20:09:34 +05:00
|
|
|
std::set<PointIndex> suround;
|
|
|
|
suround.insert(pi);
|
2024-04-10 20:34:39 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
// average only tangent component on new bl points, average whole growth
|
|
|
|
// vector otherwise
|
2024-04-10 20:34:39 +05:00
|
|
|
bool do_average_tangent = mapfrom[pi].IsValid();
|
|
|
|
correction = 0.0;
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto sei : sels) {
|
2024-09-27 20:30:36 +05:00
|
|
|
const auto &sel = mesh[sei];
|
2024-04-10 20:34:39 +05:00
|
|
|
for (auto pi1 : sel.PNums()) {
|
2024-09-27 20:30:36 +05:00
|
|
|
if (suround.count(pi1))
|
|
|
|
continue;
|
2024-04-10 20:34:39 +05:00
|
|
|
suround.insert(pi1);
|
2024-09-27 20:30:36 +05:00
|
|
|
auto gw_other = getGW(pi1) + corrections[pi1];
|
|
|
|
if (do_average_tangent) {
|
2024-04-05 20:09:34 +05:00
|
|
|
auto normal_other = getNormal(mesh[sei]);
|
2024-09-27 20:30:36 +05:00
|
|
|
auto tangent_part =
|
|
|
|
gw_other - (gw_other * normal_other) * normal_other;
|
2024-04-10 20:34:39 +05:00
|
|
|
correction += tangent_part;
|
2024-09-27 20:30:36 +05:00
|
|
|
} else {
|
2024-04-10 20:34:39 +05:00
|
|
|
correction += gw_other;
|
|
|
|
}
|
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2024-04-10 20:34:39 +05:00
|
|
|
correction *= 1.0 / suround.size();
|
2024-09-27 20:30:36 +05:00
|
|
|
if (!do_average_tangent)
|
2024-04-10 20:34:39 +05:00
|
|
|
correction -= getGW(pi);
|
2022-03-04 19:42:16 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto pi : points)
|
2024-04-10 19:04:16 +05:00
|
|
|
addGW(pi, corrections[pi]);
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2022-03-04 19:42:16 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
BoundaryLayerTool::BoundaryLayerTool(Mesh &mesh_,
|
|
|
|
const BoundaryLayerParameters ¶ms_)
|
2024-04-05 20:09:34 +05:00
|
|
|
: mesh(mesh_), topo(mesh_.GetTopology()), params(params_) {
|
|
|
|
static Timer timer("BoundaryLayerTool::ctor");
|
|
|
|
RegionTimer regt(timer);
|
2024-09-27 16:12:14 +05:00
|
|
|
ProcessParameters();
|
2024-04-05 20:09:34 +05:00
|
|
|
|
|
|
|
// for(auto & seg : mesh.LineSegments())
|
|
|
|
// seg.edgenr = seg.epgeominfo[1].edgenr;
|
|
|
|
|
|
|
|
total_height = 0.0;
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto h : par_heights)
|
|
|
|
total_height += h;
|
2024-04-05 20:09:34 +05:00
|
|
|
|
|
|
|
max_edge_nr = -1;
|
2024-09-27 20:30:36 +05:00
|
|
|
for (const auto &seg : mesh.LineSegments())
|
|
|
|
if (seg.edgenr > max_edge_nr)
|
|
|
|
max_edge_nr = seg.edgenr;
|
2024-04-05 20:09:34 +05:00
|
|
|
|
|
|
|
int ndom = mesh.GetNDomains();
|
|
|
|
ndom_old = ndom;
|
|
|
|
|
|
|
|
new_mat_nrs.SetSize(mesh.FaceDescriptors().Size() + 1);
|
|
|
|
new_mat_nrs = -1;
|
2024-09-27 16:12:14 +05:00
|
|
|
for (auto [bcname, matname] : par_new_mat) {
|
2024-04-05 20:09:34 +05:00
|
|
|
mesh.SetMaterial(++ndom, matname);
|
|
|
|
regex pattern(bcname);
|
|
|
|
for (auto i : Range(1, mesh.GetNFD() + 1)) {
|
2024-09-27 20:30:36 +05:00
|
|
|
auto &fd = mesh.GetFaceDescriptor(i);
|
|
|
|
if (regex_match(fd.GetBCName(), pattern))
|
|
|
|
new_mat_nrs[i] = ndom;
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2022-02-28 21:24:44 +05:00
|
|
|
}
|
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
if (!params.outside)
|
|
|
|
domains.Invert();
|
2022-03-07 19:58:09 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
topo.SetBuildVertex2Element(true);
|
|
|
|
mesh.UpdateTopology();
|
2009-06-19 11:43:23 +06:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
have_single_segments = HaveSingleSegments(mesh);
|
2009-06-19 11:43:23 +06:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
if (have_single_segments)
|
|
|
|
segments = BuildSegments(mesh);
|
|
|
|
else
|
|
|
|
segments = mesh.LineSegments();
|
2022-02-16 23:51:54 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
np = mesh.GetNP();
|
|
|
|
ne = mesh.GetNE();
|
|
|
|
nse = mesh.GetNSE();
|
|
|
|
nseg = segments.Size();
|
2022-02-16 23:51:54 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
p2sel = mesh.CreatePoint2SurfaceElementTable();
|
2009-06-19 11:43:23 +06:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
nfd_old = mesh.GetNFD();
|
|
|
|
moved_surfaces.SetSize(nfd_old + 1);
|
|
|
|
moved_surfaces.Clear();
|
|
|
|
si_map.SetSize(nfd_old + 1);
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto i : Range(nfd_old + 1))
|
|
|
|
si_map[i] = i;
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2020-04-23 18:44:32 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
void BoundaryLayerTool ::CreateNewFaceDescriptors() {
|
|
|
|
surfacefacs.SetSize(nfd_old + 1);
|
|
|
|
surfacefacs = 0.0;
|
|
|
|
// create new FaceDescriptors
|
|
|
|
for (auto i : Range(1, nfd_old + 1)) {
|
2024-09-27 20:30:36 +05:00
|
|
|
const auto &fd = mesh.GetFaceDescriptor(i);
|
2024-04-05 20:09:34 +05:00
|
|
|
string name = fd.GetBCName();
|
2024-09-27 16:12:14 +05:00
|
|
|
if (par_surfid.Contains(i)) {
|
2024-04-05 20:09:34 +05:00
|
|
|
if (auto isIn = domains.Test(fd.DomainIn());
|
|
|
|
isIn != domains.Test(fd.DomainOut())) {
|
|
|
|
int new_si = mesh.GetNFD() + 1;
|
|
|
|
surfacefacs[i] = isIn ? 1. : -1.;
|
|
|
|
// -1 surf nr is so that curving does not do anything
|
|
|
|
FaceDescriptor new_fd(-1, isIn ? new_mat_nrs[i] : fd.DomainIn(),
|
|
|
|
isIn ? fd.DomainOut() : new_mat_nrs[i], -1);
|
|
|
|
new_fd.SetBCProperty(new_si);
|
2024-09-04 18:53:47 +05:00
|
|
|
new_fd.SetSurfColour(fd.SurfColour());
|
2024-04-05 20:09:34 +05:00
|
|
|
mesh.AddFaceDescriptor(new_fd);
|
|
|
|
si_map[i] = new_si;
|
|
|
|
moved_surfaces.SetBit(i);
|
|
|
|
mesh.SetBCName(new_si - 1, "mapped_" + name);
|
|
|
|
}
|
|
|
|
// curving of surfaces with boundary layers will often
|
|
|
|
// result in pushed through elements, since we do not (yet)
|
|
|
|
// curvature through layers.
|
|
|
|
// Therefore we disable curving for these surfaces.
|
2024-09-27 20:30:36 +05:00
|
|
|
if (!params.keep_surfaceindex)
|
|
|
|
mesh.GetFaceDescriptor(i).SetSurfNr(-1);
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2022-03-07 19:58:09 +05:00
|
|
|
}
|
2020-04-19 23:00:06 +05:00
|
|
|
|
2024-09-27 16:12:14 +05:00
|
|
|
for (auto si : par_surfid)
|
2024-04-05 20:09:34 +05:00
|
|
|
if (surfacefacs[si] == 0.0)
|
|
|
|
throw Exception("Surface " + to_string(si) +
|
|
|
|
" is not a boundary of the domain to be grown into!");
|
|
|
|
}
|
2023-09-04 13:42:02 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
void BoundaryLayerTool ::CreateFaceDescriptorsSides() {
|
|
|
|
BitArray face_done(mesh.GetNFD() + 1);
|
|
|
|
face_done.Clear();
|
2024-09-27 20:30:36 +05:00
|
|
|
for (const auto &sel : mesh.SurfaceElements()) {
|
2024-04-05 20:09:34 +05:00
|
|
|
auto facei = sel.GetIndex();
|
2024-09-27 20:30:36 +05:00
|
|
|
if (face_done.Test(facei))
|
|
|
|
continue;
|
2024-04-05 20:09:34 +05:00
|
|
|
bool point_moved = false;
|
|
|
|
// bool point_fixed = false;
|
|
|
|
for (auto pi : sel.PNums()) {
|
2024-09-27 20:30:36 +05:00
|
|
|
if (growthvectors[pi].Length() > 0)
|
|
|
|
point_moved = true;
|
2024-04-05 20:09:34 +05:00
|
|
|
/*
|
|
|
|
else
|
|
|
|
point_fixed = true;
|
|
|
|
*/
|
|
|
|
}
|
|
|
|
if (point_moved && !moved_surfaces.Test(facei)) {
|
|
|
|
int new_si = mesh.GetNFD() + 1;
|
2024-09-27 20:30:36 +05:00
|
|
|
const auto &fd = mesh.GetFaceDescriptor(facei);
|
2024-04-05 20:09:34 +05:00
|
|
|
// auto isIn = domains.Test(fd.DomainIn());
|
|
|
|
// auto isOut = domains.Test(fd.DomainOut());
|
|
|
|
int si = params.sides_keep_surfaceindex ? facei : -1;
|
|
|
|
// domin and domout can only be set later
|
|
|
|
FaceDescriptor new_fd(si, -1, -1, si);
|
|
|
|
new_fd.SetBCProperty(new_si);
|
|
|
|
mesh.AddFaceDescriptor(new_fd);
|
|
|
|
si_map[facei] = new_si;
|
|
|
|
mesh.SetBCName(new_si - 1, fd.GetBCName());
|
|
|
|
face_done.SetBit(facei);
|
|
|
|
}
|
2022-03-07 19:58:09 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2020-04-19 23:00:06 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
void BoundaryLayerTool ::CalculateGrowthVectors() {
|
|
|
|
growthvectors.SetSize(np);
|
|
|
|
growthvectors = 0.;
|
2023-03-30 20:19:34 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto pi : mesh.Points().Range()) {
|
2024-09-27 20:30:36 +05:00
|
|
|
const auto &p = mesh[pi];
|
|
|
|
if (p.Type() == INNERPOINT)
|
|
|
|
continue;
|
2022-03-02 00:18:05 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
std::map<int, Vec<3>> normals;
|
2022-03-02 00:18:05 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
// calculate one normal vector per face (average with angles as weights for
|
|
|
|
// multiple surface elements within a face)
|
|
|
|
for (auto sei : p2sel[pi]) {
|
2024-09-27 20:30:36 +05:00
|
|
|
const auto &sel = mesh[sei];
|
2024-04-05 20:09:34 +05:00
|
|
|
auto facei = sel.GetIndex();
|
2024-09-27 20:30:36 +05:00
|
|
|
if (!par_surfid.Contains(facei))
|
|
|
|
continue;
|
2022-03-02 00:18:05 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
auto n = surfacefacs[sel.GetIndex()] * getNormal(sel);
|
2022-03-02 00:18:05 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
int itrig = sel.PNums().Pos(pi);
|
|
|
|
itrig += sel.GetNP();
|
|
|
|
auto v0 = (mesh[sel.PNumMod(itrig + 1)] - mesh[pi]).Normalize();
|
|
|
|
auto v1 = (mesh[sel.PNumMod(itrig - 1)] - mesh[pi]).Normalize();
|
2024-09-27 20:30:36 +05:00
|
|
|
if (normals.count(facei) == 0)
|
|
|
|
normals[facei] = {0., 0., 0.};
|
2024-04-05 20:09:34 +05:00
|
|
|
normals[facei] += acos(v0 * v1) * n;
|
|
|
|
}
|
2022-03-02 00:18:05 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto &[facei, n] : normals)
|
|
|
|
n *= 1.0 / n.Length();
|
2023-02-13 16:04:30 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
// combine normal vectors for each face to keep uniform distances
|
|
|
|
ArrayMem<Vec<3>, 5> ns;
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto &[facei, n] : normals) {
|
2024-04-05 20:09:34 +05:00
|
|
|
ns.Append(n);
|
|
|
|
}
|
2023-02-13 16:04:30 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
try {
|
|
|
|
growthvectors[pi] = CalcGrowthVector(ns);
|
2024-09-27 20:30:36 +05:00
|
|
|
} catch (const SpecialPointException &e) {
|
2024-04-05 20:09:34 +05:00
|
|
|
special_boundary_points.emplace(pi, normals);
|
|
|
|
growthvectors[pi] =
|
|
|
|
special_boundary_points[pi].growth_groups[0].growth_vector;
|
2022-03-02 00:18:05 +05:00
|
|
|
}
|
2022-03-07 19:58:09 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2020-04-19 23:00:06 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
Array<Array<pair<SegmentIndex, int>>, SegmentIndex>
|
|
|
|
BoundaryLayerTool ::BuildSegMap() {
|
|
|
|
// Bit array to keep track of segments already processed
|
|
|
|
BitArray segs_done(nseg + 1);
|
|
|
|
segs_done.Clear();
|
|
|
|
|
|
|
|
// map for all segments with same points
|
|
|
|
// points to pair of SegmentIndex, int
|
|
|
|
// int is type of other segment, either:
|
|
|
|
// 0 == adjacent surface grows layer
|
|
|
|
// 1 == adjacent surface doesn't grow layer, but layer ends on it
|
|
|
|
// 2 == adjacent surface is interior surface that ends on layer
|
|
|
|
// 3 == adjacent surface is exterior surface that ends on layer (not allowed
|
|
|
|
// yet)
|
|
|
|
Array<Array<pair<SegmentIndex, int>>, SegmentIndex> segmap(segments.Size());
|
|
|
|
|
|
|
|
// moved segments
|
|
|
|
is_edge_moved.SetSize(max_edge_nr + 1);
|
|
|
|
is_edge_moved = false;
|
|
|
|
|
|
|
|
// boundaries to project endings to
|
|
|
|
is_boundary_projected.SetSize(nfd_old + 1);
|
|
|
|
is_boundary_projected.Clear();
|
|
|
|
is_boundary_moved.SetSize(nfd_old + 1);
|
|
|
|
is_boundary_moved.Clear();
|
|
|
|
|
|
|
|
for (auto si : Range(segments)) {
|
2024-09-27 20:30:36 +05:00
|
|
|
if (segs_done[si])
|
|
|
|
continue;
|
|
|
|
const auto &segi = segments[si];
|
|
|
|
if (!moved_surfaces.Test(segi.si))
|
|
|
|
continue;
|
2024-04-05 20:09:34 +05:00
|
|
|
segs_done.SetBit(si);
|
|
|
|
segmap[si].Append(make_pair(si, 0));
|
|
|
|
moved_segs.Append(si);
|
|
|
|
is_edge_moved.SetBit(segi.edgenr);
|
|
|
|
for (auto sj : Range(segments)) {
|
2024-09-27 20:30:36 +05:00
|
|
|
if (segs_done.Test(sj))
|
|
|
|
continue;
|
|
|
|
const auto &segj = segments[sj];
|
2024-04-05 20:09:34 +05:00
|
|
|
if ((segi[0] == segj[0] && segi[1] == segj[1]) ||
|
|
|
|
(segi[0] == segj[1] && segi[1] == segj[0])) {
|
|
|
|
segs_done.SetBit(sj);
|
|
|
|
int type;
|
|
|
|
if (moved_surfaces.Test(segj.si)) {
|
|
|
|
type = 0;
|
|
|
|
moved_segs.Append(sj);
|
2024-09-27 20:30:36 +05:00
|
|
|
} else if (const auto &fd = mesh.GetFaceDescriptor(segj.si);
|
2024-04-05 20:09:34 +05:00
|
|
|
domains.Test(fd.DomainIn()) &&
|
|
|
|
domains.Test(fd.DomainOut())) {
|
|
|
|
type = 2;
|
|
|
|
if (fd.DomainIn() == 0 || fd.DomainOut() == 0)
|
|
|
|
is_boundary_projected.SetBit(segj.si);
|
2024-09-27 20:30:36 +05:00
|
|
|
} else if (const auto &fd = mesh.GetFaceDescriptor(segj.si);
|
2024-04-05 20:09:34 +05:00
|
|
|
!domains.Test(fd.DomainIn()) &&
|
|
|
|
!domains.Test(fd.DomainOut())) {
|
|
|
|
type = 3;
|
|
|
|
is_boundary_moved.SetBit(segj.si);
|
|
|
|
} else {
|
|
|
|
type = 1;
|
|
|
|
// in case 1 we project the growthvector onto the surface
|
|
|
|
is_boundary_projected.SetBit(segj.si);
|
|
|
|
}
|
|
|
|
segmap[si].Append(make_pair(sj, type));
|
2020-11-17 22:43:39 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2022-03-07 19:58:09 +05:00
|
|
|
}
|
2020-04-19 23:00:06 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
return segmap;
|
|
|
|
}
|
|
|
|
|
|
|
|
BitArray BoundaryLayerTool ::ProjectGrowthVectorsOnSurface() {
|
|
|
|
BitArray in_surface_direction(nfd_old + 1);
|
|
|
|
in_surface_direction.Clear();
|
|
|
|
// project growthvector on surface for inner angles
|
|
|
|
if (params.grow_edges) {
|
2024-09-27 20:30:36 +05:00
|
|
|
for (const auto &sel : mesh.SurfaceElements())
|
2024-04-05 20:09:34 +05:00
|
|
|
if (is_boundary_projected.Test(sel.GetIndex())) {
|
|
|
|
auto n = getNormal(sel);
|
|
|
|
for (auto i : Range(sel.PNums())) {
|
|
|
|
auto pi = sel.PNums()[i];
|
2024-09-27 20:30:36 +05:00
|
|
|
if (growthvectors[pi].Length2() == 0.)
|
|
|
|
continue;
|
2024-04-05 20:09:34 +05:00
|
|
|
auto next = sel.PNums()[(i + 1) % sel.GetNV()];
|
|
|
|
auto prev = sel.PNums()[i == 0 ? sel.GetNV() - 1 : i - 1];
|
|
|
|
auto v1 = (mesh[next] - mesh[pi]).Normalize();
|
|
|
|
auto v2 = (mesh[prev] - mesh[pi]).Normalize();
|
|
|
|
auto v3 = growthvectors[pi];
|
|
|
|
v3.Normalize();
|
|
|
|
auto tol = v1.Length() * 1e-12;
|
|
|
|
if ((v1 * v3 > -tol) && (v2 * v3 > -tol))
|
|
|
|
in_surface_direction.SetBit(sel.GetIndex());
|
|
|
|
else
|
|
|
|
continue;
|
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
if (!par_project_boundaries.Contains(sel.GetIndex()))
|
|
|
|
continue;
|
|
|
|
auto &g = growthvectors[pi];
|
2024-04-05 20:09:34 +05:00
|
|
|
auto ng = n * g;
|
|
|
|
auto gg = g * g;
|
|
|
|
auto nn = n * n;
|
|
|
|
// if(fabs(ng*ng-nn*gg) < 1e-12 || fabs(ng) < 1e-12) continue;
|
|
|
|
auto a = -ng * ng / (ng * ng - nn * gg);
|
|
|
|
auto b = ng * gg / (ng * ng - nn * gg);
|
|
|
|
g += a * g + b * n;
|
|
|
|
}
|
2020-11-17 22:43:39 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
} else {
|
2024-09-27 20:30:36 +05:00
|
|
|
for (const auto &seg : segments) {
|
2024-04-05 20:09:34 +05:00
|
|
|
int count = 0;
|
2024-09-27 20:30:36 +05:00
|
|
|
for (const auto &seg2 : segments)
|
2024-04-05 20:09:34 +05:00
|
|
|
if (((seg[0] == seg2[0] && seg[1] == seg2[1]) ||
|
|
|
|
(seg[0] == seg2[1] && seg[1] == seg2[0])) &&
|
2024-09-27 16:12:14 +05:00
|
|
|
par_surfid.Contains(seg2.si))
|
2024-04-05 20:09:34 +05:00
|
|
|
count++;
|
|
|
|
if (count == 1) {
|
|
|
|
growthvectors[seg[0]] = {0., 0., 0.};
|
|
|
|
growthvectors[seg[1]] = {0., 0., 0.};
|
2020-11-17 22:43:39 +05:00
|
|
|
}
|
2022-03-07 19:58:09 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2022-03-07 19:58:09 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
return in_surface_direction;
|
|
|
|
}
|
|
|
|
|
|
|
|
void BoundaryLayerTool ::InterpolateGrowthVectors() {
|
|
|
|
int new_max_edge_nr = max_edge_nr;
|
2024-09-27 20:30:36 +05:00
|
|
|
for (const auto &seg : segments)
|
|
|
|
if (seg.edgenr > new_max_edge_nr)
|
|
|
|
new_max_edge_nr = seg.edgenr;
|
|
|
|
for (const auto &seg : new_segments)
|
|
|
|
if (seg.edgenr > new_max_edge_nr)
|
|
|
|
new_max_edge_nr = seg.edgenr;
|
2024-04-05 20:09:34 +05:00
|
|
|
|
|
|
|
auto getGW = [&](PointIndex pi) -> Vec<3> {
|
|
|
|
if (growth_vector_map.count(pi) == 0)
|
|
|
|
growth_vector_map[pi] = {&growthvectors[pi], total_height};
|
|
|
|
auto [gw, height] = growth_vector_map[pi];
|
|
|
|
return height * (*gw);
|
|
|
|
};
|
|
|
|
auto addGW = [&](PointIndex pi, Vec<3> vec) {
|
|
|
|
if (growth_vector_map.count(pi) == 0)
|
|
|
|
growth_vector_map[pi] = {&growthvectors[pi], total_height};
|
|
|
|
auto [gw, height] = growth_vector_map[pi];
|
|
|
|
*gw += 1.0 / height * vec;
|
|
|
|
};
|
2023-11-29 13:19:27 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
// interpolate tangential component of growth vector along edge
|
2024-09-27 20:30:36 +05:00
|
|
|
if (max_edge_nr < new_max_edge_nr)
|
|
|
|
for (auto edgenr : Range(max_edge_nr + 1, new_max_edge_nr)) {
|
|
|
|
// cout << "SEARCH EDGE " << edgenr +1 << endl;
|
|
|
|
// if(!is_edge_moved[edgenr+1]) continue;
|
|
|
|
|
|
|
|
// build sorted list of edge
|
|
|
|
Array<PointIndex> points;
|
|
|
|
// find first vertex on edge
|
|
|
|
double edge_len = 0.;
|
|
|
|
auto is_end_point = [&](PointIndex pi) {
|
|
|
|
// if(mesh[pi].Type() == FIXEDPOINT)
|
|
|
|
// return true;
|
|
|
|
// return false;
|
|
|
|
auto segs = topo.GetVertexSegments(pi);
|
|
|
|
if (segs.Size() == 1)
|
|
|
|
return true;
|
|
|
|
auto first_edgenr = mesh[segs[0]].edgenr;
|
|
|
|
for (auto segi : segs)
|
|
|
|
if (mesh[segi].edgenr != first_edgenr)
|
|
|
|
return true;
|
|
|
|
return false;
|
|
|
|
};
|
|
|
|
|
|
|
|
bool any_grows = false;
|
|
|
|
|
|
|
|
for (const auto &seg : segments) {
|
|
|
|
if (seg.edgenr - 1 == edgenr) {
|
|
|
|
if (getGW(seg[0]).Length2() != 0 || getGW(seg[1]).Length2() != 0)
|
|
|
|
any_grows = true;
|
|
|
|
if (points.Size() == 0 &&
|
|
|
|
(is_end_point(seg[0]) || is_end_point(seg[1]))) {
|
|
|
|
PointIndex seg0 = seg[0], seg1 = seg[1];
|
|
|
|
if (is_end_point(seg[1]))
|
|
|
|
Swap(seg0, seg1);
|
|
|
|
points.Append(seg0);
|
|
|
|
points.Append(seg1);
|
|
|
|
edge_len += (mesh[seg[1]] - mesh[seg[0]]).Length();
|
|
|
|
}
|
2023-11-29 13:19:27 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2022-06-07 17:51:41 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
if (!any_grows) {
|
|
|
|
// cout << "skip edge " << edgenr+1 << endl;
|
|
|
|
continue;
|
|
|
|
}
|
2023-11-29 13:19:27 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
if (!points.Size())
|
|
|
|
throw Exception("Could not find startpoint for edge " +
|
|
|
|
ToString(edgenr));
|
|
|
|
|
|
|
|
while (true) {
|
|
|
|
bool point_found = false;
|
|
|
|
for (auto si : topo.GetVertexSegments(points.Last())) {
|
|
|
|
const auto &seg = mesh[si];
|
|
|
|
if (seg.edgenr - 1 != edgenr)
|
|
|
|
continue;
|
|
|
|
if (seg[0] == points.Last() && points[points.Size() - 2] != seg[1]) {
|
|
|
|
edge_len += (mesh[points.Last()] - mesh[seg[1]]).Length();
|
|
|
|
points.Append(seg[1]);
|
|
|
|
point_found = true;
|
|
|
|
break;
|
|
|
|
} else if (seg[1] == points.Last() &&
|
|
|
|
points[points.Size() - 2] != seg[0]) {
|
|
|
|
edge_len += (mesh[points.Last()] - mesh[seg[0]]).Length();
|
|
|
|
points.Append(seg[0]);
|
|
|
|
point_found = true;
|
|
|
|
break;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
if (is_end_point(points.Last()))
|
2024-04-05 20:09:34 +05:00
|
|
|
break;
|
2024-09-27 20:30:36 +05:00
|
|
|
if (!point_found) {
|
|
|
|
throw Exception(
|
|
|
|
string(
|
|
|
|
"Could not find connected list of line segments for edge ") +
|
|
|
|
edgenr);
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2022-02-28 12:29:22 +05:00
|
|
|
}
|
2023-11-29 13:19:27 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
if (getGW(points[0]).Length2() == 0 &&
|
|
|
|
getGW(points.Last()).Length2() == 0)
|
|
|
|
continue;
|
|
|
|
// cout << "Points to average " << endl << points << endl;
|
|
|
|
|
|
|
|
// tangential part of growth vectors
|
|
|
|
auto t1 = (mesh[points[1]] - mesh[points[0]]).Normalize();
|
|
|
|
auto gt1 = getGW(points[0]) * t1 * t1;
|
|
|
|
auto t2 =
|
|
|
|
(mesh[points.Last()] - mesh[points[points.Size() - 2]]).Normalize();
|
|
|
|
auto gt2 = getGW(points.Last()) * t2 * t2;
|
|
|
|
|
|
|
|
// if(!is_edge_moved[edgenr+1])
|
|
|
|
// {
|
|
|
|
// if(getGW(points[0]) * (mesh[points[1]] - mesh[points[0]]) < 0)
|
|
|
|
// gt1 = 0.;
|
|
|
|
// if(getGW(points.Last()) * (mesh[points[points.Size()-2]] -
|
|
|
|
// mesh[points.Last()]) < 0)
|
|
|
|
// gt2 = 0.;
|
|
|
|
// }
|
|
|
|
|
|
|
|
double len = 0.;
|
|
|
|
for (auto i : IntRange(1, points.Size() - 1)) {
|
|
|
|
auto pi = points[i];
|
|
|
|
len += (mesh[pi] - mesh[points[i - 1]]).Length();
|
|
|
|
auto t = getEdgeTangent(pi, edgenr);
|
|
|
|
auto lam = len / edge_len;
|
|
|
|
auto interpol = (1 - lam) * (gt1 * t) * t + lam * (gt2 * t) * t;
|
|
|
|
addGW(pi, interpol);
|
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2023-12-11 23:05:04 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
InterpolateSurfaceGrowthVectors();
|
|
|
|
}
|
2023-11-29 13:19:27 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
void BoundaryLayerTool ::InsertNewElements(
|
|
|
|
FlatArray<Array<pair<SegmentIndex, int>>, SegmentIndex> segmap,
|
2024-09-27 20:30:36 +05:00
|
|
|
const BitArray &in_surface_direction) {
|
2024-04-05 20:09:34 +05:00
|
|
|
static Timer timer("BoundaryLayerTool::InsertNewElements");
|
|
|
|
RegionTimer rt(timer);
|
|
|
|
mapto.SetSize(0);
|
|
|
|
mapto.SetSize(np);
|
2024-09-04 18:49:02 +05:00
|
|
|
mapfrom.SetSize(mesh.GetNP());
|
|
|
|
mapfrom = PointIndex::INVALID;
|
2024-04-05 20:09:34 +05:00
|
|
|
|
|
|
|
auto changed_domains = domains;
|
2024-09-27 20:30:36 +05:00
|
|
|
if (!params.outside)
|
|
|
|
changed_domains.Invert();
|
2024-04-05 20:09:34 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
auto &identifications = mesh.GetIdentifications();
|
2024-04-05 20:09:34 +05:00
|
|
|
const int identnr = identifications.GetNr("boundarylayer");
|
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
auto add_points = [&](PointIndex pi, Vec<3> &growth_vector,
|
|
|
|
Array<PointIndex> &new_points) {
|
2024-04-05 20:09:34 +05:00
|
|
|
Point<3> p = mesh[pi];
|
|
|
|
PointIndex pi_last = pi;
|
2024-07-26 21:09:37 +05:00
|
|
|
double height = 0.0;
|
2024-09-27 16:12:14 +05:00
|
|
|
for (auto i : Range(par_heights)) {
|
|
|
|
height += par_heights[i];
|
2024-04-05 20:09:34 +05:00
|
|
|
auto pi_new = mesh.AddPoint(p);
|
2024-09-04 18:49:02 +05:00
|
|
|
mapfrom.Append(pi);
|
2024-04-05 20:09:34 +05:00
|
|
|
new_points.Append(pi_new);
|
2024-07-26 21:09:37 +05:00
|
|
|
growth_vector_map[pi_new] = {&growth_vector, height};
|
2024-09-27 20:30:36 +05:00
|
|
|
if (special_boundary_points.count(pi) > 0)
|
|
|
|
mesh.AddLockedPoint(pi_new);
|
2024-04-05 20:09:34 +05:00
|
|
|
pi_last = pi_new;
|
2023-11-29 13:19:27 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
};
|
2023-11-29 13:19:27 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
// insert new points
|
|
|
|
for (PointIndex pi = 1; pi <= np; pi++) {
|
|
|
|
if (growthvectors[pi].Length2() != 0) {
|
|
|
|
if (special_boundary_points.count(pi)) {
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto &group : special_boundary_points[pi].growth_groups)
|
2024-04-05 20:09:34 +05:00
|
|
|
add_points(pi, group.growth_vector, group.new_points);
|
|
|
|
} else
|
|
|
|
add_points(pi, growthvectors[pi], mapto[pi]);
|
|
|
|
}
|
|
|
|
}
|
2023-11-29 13:19:27 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
// 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) {
|
2024-09-27 20:30:36 +05:00
|
|
|
if (layer == -1)
|
|
|
|
layer = par_heights.Size() - 1;
|
2024-04-05 20:09:34 +05:00
|
|
|
if (special_boundary_points.count(pi))
|
|
|
|
return special_boundary_points[pi].growth_groups[group].new_points[layer];
|
|
|
|
else
|
|
|
|
return mapto[pi][layer];
|
|
|
|
};
|
2023-11-29 13:19:27 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
auto hasMoved = [&](PointIndex pi) {
|
|
|
|
return mapto[pi].Size() > 0 || special_boundary_points.count(pi);
|
|
|
|
};
|
2023-11-29 13:19:27 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
auto numGroups = [&](PointIndex pi) -> size_t {
|
|
|
|
if (special_boundary_points.count(pi))
|
|
|
|
return special_boundary_points[pi].growth_groups.Size();
|
|
|
|
else
|
|
|
|
return 1;
|
|
|
|
};
|
2023-11-29 13:19:27 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
auto getGroups = [&](PointIndex pi, int face_index) -> Array<int> {
|
|
|
|
auto n = numGroups(pi);
|
|
|
|
Array<int> groups;
|
|
|
|
if (n == 1) {
|
|
|
|
groups.Append(0);
|
2023-11-29 13:19:27 +05:00
|
|
|
return groups;
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2024-09-27 20:30:36 +05:00
|
|
|
const auto &all_groups = special_boundary_points[pi].growth_groups;
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto i : Range(n))
|
2024-09-27 20:30:36 +05:00
|
|
|
if (all_groups[i].faces.Contains(face_index))
|
|
|
|
groups.Append(i);
|
2024-04-05 20:09:34 +05:00
|
|
|
// cout << "groups " << pi << ", " << face_index << endl << groups;
|
|
|
|
return groups;
|
|
|
|
};
|
2020-04-19 23:00:06 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
// add 2d quads on required surfaces
|
|
|
|
map<pair<PointIndex, PointIndex>, int> seg2edge;
|
|
|
|
map<int, int> edge_map;
|
|
|
|
int edge_nr = max_edge_nr;
|
|
|
|
auto getEdgeNr = [&](int ei) {
|
2024-09-27 20:30:36 +05:00
|
|
|
if (edge_map.count(ei) == 0)
|
|
|
|
edge_map[ei] = ++edge_nr;
|
2024-04-05 20:09:34 +05:00
|
|
|
return edge_map[ei];
|
|
|
|
};
|
|
|
|
if (params.grow_edges) {
|
|
|
|
for (auto sei : moved_segs) {
|
|
|
|
// copy here since we will add segments and this would
|
|
|
|
// invalidate a reference!
|
|
|
|
// auto segi = segments[sei];
|
|
|
|
for (auto [sej, type] : segmap[sei]) {
|
|
|
|
auto segj = segments[sej];
|
|
|
|
if (type == 0) {
|
|
|
|
auto addSegment = [&](PointIndex p0, PointIndex p1,
|
|
|
|
bool extra_edge_nr = false) {
|
|
|
|
Segment s;
|
|
|
|
s[0] = p0;
|
|
|
|
s[1] = p1;
|
|
|
|
s[2] = PointIndex::INVALID;
|
|
|
|
auto pair =
|
|
|
|
s[0] < s[1] ? make_pair(s[0], s[1]) : make_pair(s[1], s[0]);
|
|
|
|
if (extra_edge_nr)
|
|
|
|
s.edgenr = ++edge_nr;
|
|
|
|
else
|
|
|
|
s.edgenr = getEdgeNr(segj.edgenr);
|
|
|
|
s.si = si_map[segj.si];
|
|
|
|
new_segments.Append(s);
|
|
|
|
// cout << __LINE__ <<"\t" << s << endl;
|
|
|
|
return s;
|
|
|
|
};
|
|
|
|
|
|
|
|
auto p0 = segj[0], p1 = segj[1];
|
|
|
|
auto g0 = getGroups(p0, segj.si);
|
|
|
|
auto g1 = getGroups(p1, segj.si);
|
|
|
|
|
|
|
|
if (g0.Size() == 1 && g1.Size() == 1)
|
|
|
|
auto s =
|
|
|
|
addSegment(newPoint(p0, -1, g0[0]), newPoint(p1, -1, g1[0]));
|
|
|
|
else {
|
|
|
|
if (g0.Size() == 2)
|
|
|
|
addSegment(newPoint(p0, -1, g0[0]), newPoint(p0, -1, g0[1]));
|
|
|
|
if (g1.Size() == 2)
|
|
|
|
addSegment(newPoint(p1, -1, g1[0]), newPoint(p1, -1, g1[1]));
|
2020-11-17 22:43:39 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
|
|
|
// here we need to grow the quad elements
|
|
|
|
else if (type == 1) {
|
|
|
|
PointIndex pp1 = segj[1];
|
|
|
|
PointIndex pp2 = segj[0];
|
|
|
|
if (in_surface_direction.Test(segj.si)) {
|
|
|
|
Swap(pp1, pp2);
|
|
|
|
is_boundary_moved.SetBit(segj.si);
|
|
|
|
}
|
|
|
|
PointIndex p1 = pp1;
|
|
|
|
PointIndex p2 = pp2;
|
|
|
|
PointIndex p3, p4;
|
|
|
|
Segment s0;
|
|
|
|
s0[0] = p1;
|
|
|
|
s0[1] = p2;
|
|
|
|
s0[2] = PointIndex::INVALID;
|
|
|
|
s0.edgenr = segj.edgenr;
|
|
|
|
s0.si = segj.si;
|
|
|
|
new_segments.Append(s0);
|
|
|
|
|
2024-09-27 16:12:14 +05:00
|
|
|
for (auto i : Range(par_heights)) {
|
2024-04-05 20:09:34 +05:00
|
|
|
Element2d sel(QUAD);
|
|
|
|
p3 = newPoint(pp2, i);
|
|
|
|
p4 = newPoint(pp1, i);
|
|
|
|
sel[0] = p1;
|
|
|
|
sel[1] = p2;
|
|
|
|
sel[2] = p3;
|
|
|
|
sel[3] = p4;
|
|
|
|
for (auto i : Range(4)) {
|
|
|
|
sel.GeomInfo()[i].u = 0.0;
|
|
|
|
sel.GeomInfo()[i].v = 0.0;
|
|
|
|
}
|
|
|
|
sel.SetIndex(si_map[segj.si]);
|
|
|
|
mesh.AddSurfaceElement(sel);
|
|
|
|
|
|
|
|
// TODO: Too many, would be enough to only add outermost ones
|
|
|
|
Segment s1;
|
|
|
|
s1[0] = p2;
|
|
|
|
s1[1] = p3;
|
|
|
|
s1[2] = PointIndex::INVALID;
|
|
|
|
auto pair = make_pair(p2, p3);
|
|
|
|
s1.edgenr = getEdgeNr(segj.edgenr);
|
|
|
|
s1.si = segj.si;
|
|
|
|
// new_segments.Append(s1);
|
|
|
|
Segment s2;
|
|
|
|
s2[0] = p4;
|
|
|
|
s2[1] = p1;
|
|
|
|
s2[2] = PointIndex::INVALID;
|
|
|
|
pair = make_pair(p1, p4);
|
|
|
|
s2.edgenr = getEdgeNr(segj.edgenr);
|
|
|
|
s2.si = segj.si;
|
|
|
|
// new_segments.Append(s2);
|
|
|
|
p1 = p4;
|
|
|
|
p2 = p3;
|
|
|
|
}
|
|
|
|
Segment s3;
|
|
|
|
s3[0] = p3;
|
|
|
|
s3[1] = p4;
|
|
|
|
s3[2] = PointIndex::INVALID;
|
|
|
|
auto pair = p3 < p4 ? make_pair(p3, p4) : make_pair(p4, p3);
|
|
|
|
s3.edgenr = getEdgeNr(segj.edgenr);
|
|
|
|
s3.si = segj.si;
|
|
|
|
new_segments.Append(s3);
|
|
|
|
}
|
2020-11-17 22:43:39 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
|
|
|
}
|
2020-04-19 23:00:06 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
auto getClosestGroup = [&](PointIndex pi, SurfaceElementIndex sei) {
|
|
|
|
auto n = numGroups(pi);
|
2024-09-27 20:30:36 +05:00
|
|
|
if (n == 1)
|
|
|
|
return 0;
|
|
|
|
const auto &sel = mesh[sei];
|
2024-09-04 18:49:02 +05:00
|
|
|
auto groups = getGroups(pi, sel.GetIndex());
|
2024-09-27 20:30:36 +05:00
|
|
|
if (groups.Size() == 1)
|
|
|
|
return groups[0];
|
2024-09-04 18:49:02 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
auto &growth_groups = special_boundary_points[pi].growth_groups;
|
2024-09-04 18:49:02 +05:00
|
|
|
|
|
|
|
auto vdir = Center(mesh[sel[0]], mesh[sel[1]], mesh[sel[2]]) - mesh[pi];
|
|
|
|
auto dot = vdir * special_boundary_points[pi].separating_direction;
|
|
|
|
|
|
|
|
return dot > 0 ? 1 : 0;
|
2024-04-05 20:09:34 +05:00
|
|
|
};
|
2023-12-11 23:05:04 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
BitArray fixed_points(np + 1);
|
|
|
|
fixed_points.Clear();
|
|
|
|
BitArray moveboundarypoint(np + 1);
|
|
|
|
moveboundarypoint.Clear();
|
|
|
|
auto p2el = mesh.CreatePoint2ElementTable();
|
|
|
|
for (SurfaceElementIndex si = 0; si < nse; si++) {
|
|
|
|
// copy because surfaceels array will be resized!
|
|
|
|
const auto sel = mesh[si];
|
|
|
|
if (moved_surfaces.Test(sel.GetIndex())) {
|
|
|
|
Array<PointIndex> points(sel.PNums());
|
2024-09-27 20:30:36 +05:00
|
|
|
if (surfacefacs[sel.GetIndex()] > 0)
|
|
|
|
Swap(points[0], points[2]);
|
2024-04-05 20:09:34 +05:00
|
|
|
ArrayMem<int, 4> groups(points.Size());
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto i : Range(points))
|
|
|
|
groups[i] = getClosestGroup(sel[i], si);
|
2024-04-05 20:09:34 +05:00
|
|
|
bool add_volume_element = true;
|
|
|
|
for (auto pi : sel.PNums())
|
2024-09-27 20:30:36 +05:00
|
|
|
if (numGroups(pi) > 1)
|
|
|
|
add_volume_element = false;
|
2024-09-27 16:12:14 +05:00
|
|
|
for (auto j : Range(par_heights)) {
|
2024-04-05 20:09:34 +05:00
|
|
|
auto eltype = points.Size() == 3 ? PRISM : HEX;
|
|
|
|
Element el(eltype);
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto i : Range(points))
|
|
|
|
el[i] = points[i];
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto i : Range(points))
|
|
|
|
points[i] = newPoint(sel.PNums()[i], j, groups[i]);
|
2024-09-27 20:30:36 +05:00
|
|
|
if (surfacefacs[sel.GetIndex()] > 0)
|
|
|
|
Swap(points[0], points[2]);
|
|
|
|
for (auto i : Range(points))
|
|
|
|
el[sel.PNums().Size() + i] = points[i];
|
2024-04-05 20:09:34 +05:00
|
|
|
auto new_index = new_mat_nrs[sel.GetIndex()];
|
|
|
|
if (new_index == -1)
|
|
|
|
throw Exception("Boundary " + ToString(sel.GetIndex()) +
|
|
|
|
" with name " + mesh.GetBCName(sel.GetIndex() - 1) +
|
|
|
|
" extruded, but no new material specified for it!");
|
|
|
|
el.SetIndex(new_mat_nrs[sel.GetIndex()]);
|
|
|
|
if (add_volume_element)
|
|
|
|
mesh.AddVolumeElement(el);
|
|
|
|
else {
|
|
|
|
// Let the volume mesher fill the hole with pyramids/tets
|
|
|
|
// To insert pyramids, we need close surface identifications on open
|
|
|
|
// quads
|
|
|
|
for (auto i : Range(points))
|
|
|
|
if (numGroups(sel[i]) == 1)
|
|
|
|
identifications.Add(el[i], el[i + points.Size()], identnr);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
Element2d newel = sel;
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto i : Range(points))
|
|
|
|
newel[i] = newPoint(sel[i], -1, groups[i]);
|
2024-04-05 20:09:34 +05:00
|
|
|
newel.SetIndex(si_map[sel.GetIndex()]);
|
|
|
|
mesh.AddSurfaceElement(newel);
|
|
|
|
|
|
|
|
// also move volume element adjacent to this surface element accordingly
|
|
|
|
ElementIndex ei = -1;
|
|
|
|
// if(groups[0] || groups[1] || groups[2])
|
|
|
|
// for(auto ei_ : p2el[sel.PNums()[0]])
|
|
|
|
// {
|
|
|
|
// const auto & el = mesh[ei_];
|
|
|
|
// // if(!domains.Test(el.GetIndex())) continue;
|
|
|
|
// cout << "check " << ei_ << "\t" << el << "\t" << sel << endl;
|
|
|
|
// auto pnums = el.PNums();
|
|
|
|
// if(pnums.Contains(sel[1]) && pnums.Contains(sel[2])) {
|
|
|
|
// ei = ei_;
|
|
|
|
// break;
|
|
|
|
// }
|
|
|
|
// }
|
|
|
|
if (ei != -1) {
|
2024-09-27 20:30:36 +05:00
|
|
|
auto &el = mesh[ei];
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto i : Range(el.GetNP()))
|
|
|
|
for (auto j : Range(3)) {
|
|
|
|
if (groups[j] && el[i] == sel[j]) {
|
|
|
|
el[i] = newel[j];
|
|
|
|
break;
|
2023-12-11 23:05:04 +05:00
|
|
|
}
|
2020-11-17 22:43:39 +05:00
|
|
|
}
|
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
} else {
|
|
|
|
bool has_moved = false;
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto p : sel.PNums())
|
|
|
|
has_moved |= hasMoved(p);
|
2024-04-05 20:09:34 +05:00
|
|
|
if (has_moved)
|
|
|
|
for (auto p : sel.PNums()) {
|
|
|
|
if (hasMoved(p)) {
|
|
|
|
fixed_points.SetBit(p);
|
|
|
|
if (is_boundary_moved.Test(sel.GetIndex()))
|
|
|
|
moveboundarypoint.SetBit(p);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
if (is_boundary_moved.Test(sel.GetIndex())) {
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto &p : mesh[si].PNums())
|
|
|
|
if (hasMoved(p))
|
|
|
|
p = newPoint(p);
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
|
|
|
}
|
2015-01-09 02:18:22 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
for (SegmentIndex sei = 0; sei < nseg; sei++) {
|
2024-09-27 20:30:36 +05:00
|
|
|
auto &seg = segments[sei];
|
2024-04-05 20:09:34 +05:00
|
|
|
if (is_boundary_moved.Test(seg.si))
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto &p : seg.PNums())
|
|
|
|
if (hasMoved(p))
|
|
|
|
p = newPoint(p);
|
2024-04-05 20:09:34 +05:00
|
|
|
// else if(hasMoved(seg[0]) || hasMoved(seg[1]))
|
|
|
|
// {
|
|
|
|
// auto tangent = mesh[seg[1]] - mesh[seg[0]];
|
|
|
|
// if(hasMoved(seg[0]) && growthvectors[seg[0]] * tangent > 0)
|
|
|
|
// seg[0] = newPoint(seg[0]);
|
|
|
|
// if(hasMoved(seg[1]) && growthvectors[seg[1]] * tangent < 0)
|
|
|
|
// seg[1] = newPoint(seg[1]);
|
|
|
|
// }
|
|
|
|
}
|
2023-11-29 13:19:27 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
// fill holes in surface mesh at special boundary points (i.e. points with >=4
|
|
|
|
// adjacent boundary faces)
|
2024-04-05 20:09:34 +05:00
|
|
|
auto p2sel = mesh.CreatePoint2SurfaceElementTable();
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto &[special_pi, special_point] : special_boundary_points) {
|
2024-09-04 18:49:02 +05:00
|
|
|
if (special_point.growth_groups.Size() != 2)
|
|
|
|
throw Exception("special_point.growth_groups.Size() != 2");
|
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
// Special points are split into two new points, when mapping a surface
|
|
|
|
// element, we choose the closer one to the center. Now, find points which
|
|
|
|
// 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.
|
2024-09-04 18:49:02 +05:00
|
|
|
std::map<int, std::array<std::set<PointIndex>, 2>> close_group;
|
|
|
|
for (auto sei : p2sel[special_pi]) {
|
2024-09-27 20:30:36 +05:00
|
|
|
const auto &sel = mesh[sei];
|
2024-09-04 18:49:02 +05:00
|
|
|
for (auto p : sel.PNums())
|
2024-09-27 20:30:36 +05:00
|
|
|
if (p != special_pi)
|
|
|
|
close_group[sel.GetIndex()][getClosestGroup(special_pi, sei)].insert(
|
|
|
|
p);
|
2024-09-04 18:49:02 +05:00
|
|
|
}
|
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto [fi, groups] : close_group) {
|
2024-09-04 18:49:02 +05:00
|
|
|
const auto mapped_fi = si_map[fi];
|
|
|
|
std::set<PointIndex> common_points;
|
|
|
|
for (auto pi : groups[0])
|
2024-09-27 20:30:36 +05:00
|
|
|
if (groups[1].count(pi) == 1)
|
2024-09-04 18:49:02 +05:00
|
|
|
common_points.insert(pi);
|
2024-09-27 20:30:36 +05:00
|
|
|
if (common_points.size() > 0) {
|
2024-09-04 18:49:02 +05:00
|
|
|
auto pi_common = mapto[*common_points.begin()].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();
|
|
|
|
for (auto sei : p2sel[pi_common]) {
|
2024-09-27 20:30:36 +05:00
|
|
|
if (mesh[sei].GetIndex() == mapped_fi &&
|
|
|
|
mesh[sei].PNums().Contains(new_special_pi0)) {
|
2024-09-04 18:49:02 +05:00
|
|
|
auto sel = mesh[sei];
|
|
|
|
sel.Invert();
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto &pi : sel.PNums())
|
|
|
|
if (pi != pi_common && pi != new_special_pi0)
|
|
|
|
pi = new_special_pi1;
|
2024-09-04 18:49:02 +05:00
|
|
|
mesh.AddSurfaceElement(sel);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto &[pi, special_point] : special_boundary_points) {
|
2024-04-05 20:09:34 +05:00
|
|
|
if (special_point.growth_groups.Size() != 2)
|
2023-11-29 13:19:27 +05:00
|
|
|
throw Exception("special_point.growth_groups.Size() != 2");
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto igroup : Range(2)) {
|
2024-09-27 20:30:36 +05:00
|
|
|
auto &group = special_point.growth_groups[igroup];
|
2023-11-29 13:19:27 +05:00
|
|
|
std::set<int> faces;
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto face : group.faces)
|
|
|
|
faces.insert(si_map[face]);
|
2023-11-29 13:19:27 +05:00
|
|
|
auto pi_new = group.new_points.Last();
|
2024-04-05 20:09:34 +05:00
|
|
|
auto pi_new_other =
|
|
|
|
special_point.growth_groups[1 - igroup].new_points.Last();
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto sei : p2sel[pi_new])
|
|
|
|
faces.erase(mesh[sei].GetIndex());
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto face : faces)
|
|
|
|
for (auto seg : new_segments) {
|
2024-09-27 20:30:36 +05:00
|
|
|
if ( // seg.si == face
|
2024-04-05 20:09:34 +05:00
|
|
|
(seg[0] == pi_new || seg[1] == pi_new) &&
|
|
|
|
(seg[0] != pi_new_other && seg[1] != pi_new_other)) {
|
2023-11-29 13:19:27 +05:00
|
|
|
bool is_correct_face = false;
|
|
|
|
auto pi_other = seg[0] == pi_new ? seg[1] : seg[0];
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto sei : p2sel[pi_other]) {
|
|
|
|
if (mesh[sei].GetIndex() == face) {
|
2023-11-29 13:19:27 +05:00
|
|
|
is_correct_face = true;
|
|
|
|
break;
|
|
|
|
}
|
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
if (is_correct_face) {
|
2023-11-29 13:19:27 +05:00
|
|
|
Element2d sel;
|
|
|
|
sel[0] = seg[1];
|
|
|
|
sel[1] = seg[0];
|
|
|
|
sel[2] = pi_new_other;
|
|
|
|
sel.SetIndex(face);
|
|
|
|
mesh.AddSurfaceElement(sel);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2015-01-09 02:18:22 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
for (ElementIndex ei = 0; ei < ne; ei++) {
|
|
|
|
auto el = mesh[ei];
|
|
|
|
ArrayMem<PointIndex, 4> fixed;
|
|
|
|
ArrayMem<PointIndex, 4> moved;
|
|
|
|
bool moved_bnd = false;
|
2024-09-27 20:30:36 +05:00
|
|
|
for (const auto &p : el.PNums()) {
|
|
|
|
if (fixed_points.Test(p))
|
|
|
|
fixed.Append(p);
|
|
|
|
if (hasMoved(p))
|
|
|
|
moved.Append(p);
|
|
|
|
if (moveboundarypoint.Test(p))
|
|
|
|
moved_bnd = true;
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2020-11-24 03:48:49 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
bool do_move, do_insert;
|
|
|
|
if (changed_domains.Test(el.GetIndex())) {
|
|
|
|
do_move = fixed.Size() && moved_bnd;
|
|
|
|
do_insert = do_move;
|
|
|
|
} else {
|
|
|
|
do_move = !fixed.Size() || moved_bnd;
|
|
|
|
do_insert = !do_move;
|
|
|
|
}
|
|
|
|
|
2024-04-10 14:13:57 +05:00
|
|
|
// if (do_move) {
|
|
|
|
// for (auto& p : mesh[ei].PNums())
|
|
|
|
// if (hasMoved(p)) {
|
|
|
|
// if (special_boundary_points.count(p)) {
|
|
|
|
// auto& special_point = special_boundary_points[p];
|
|
|
|
// auto& group = special_point.growth_groups[0];
|
|
|
|
// p = group.new_points.Last();
|
|
|
|
// } else
|
|
|
|
// p = newPoint(p);
|
|
|
|
// }
|
|
|
|
// }
|
2024-04-05 20:09:34 +05:00
|
|
|
if (do_insert) {
|
|
|
|
if (el.GetType() == TET) {
|
2024-09-27 20:30:36 +05:00
|
|
|
if (moved.Size() == 3) // inner corner
|
2024-04-05 20:09:34 +05:00
|
|
|
{
|
|
|
|
PointIndex p1 = moved[0];
|
|
|
|
PointIndex p2 = moved[1];
|
|
|
|
PointIndex p3 = moved[2];
|
|
|
|
auto v1 = mesh[p1];
|
|
|
|
auto n = Cross(mesh[p2] - v1, mesh[p3] - v1);
|
|
|
|
auto d = mesh[newPoint(p1, 0)] - v1;
|
2024-09-27 20:30:36 +05:00
|
|
|
if (n * d > 0)
|
|
|
|
Swap(p2, p3);
|
2024-04-05 20:09:34 +05:00
|
|
|
PointIndex p4 = p1;
|
|
|
|
PointIndex p5 = p2;
|
|
|
|
PointIndex p6 = p3;
|
2024-09-27 16:12:14 +05:00
|
|
|
for (auto i : Range(par_heights)) {
|
2024-04-05 20:09:34 +05:00
|
|
|
Element nel(PRISM);
|
|
|
|
nel[0] = p4;
|
|
|
|
nel[1] = p5;
|
|
|
|
nel[2] = p6;
|
|
|
|
p4 = newPoint(p1, i);
|
|
|
|
p5 = newPoint(p2, i);
|
|
|
|
p6 = newPoint(p3, i);
|
|
|
|
nel[3] = p4;
|
|
|
|
nel[4] = p5;
|
|
|
|
nel[5] = p6;
|
|
|
|
nel.SetIndex(el.GetIndex());
|
|
|
|
mesh.AddVolumeElement(nel);
|
2020-11-24 03:48:49 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
|
|
|
if (moved.Size() == 2) {
|
|
|
|
if (fixed.Size() == 1) {
|
|
|
|
PointIndex p1 = moved[0];
|
|
|
|
PointIndex p2 = moved[1];
|
2024-09-27 16:12:14 +05:00
|
|
|
for (auto i : Range(par_heights)) {
|
2024-04-05 20:09:34 +05:00
|
|
|
PointIndex p3 = newPoint(moved[1], i);
|
|
|
|
PointIndex p4 = newPoint(moved[0], i);
|
|
|
|
Element nel(PYRAMID);
|
|
|
|
nel[0] = p1;
|
|
|
|
nel[1] = p2;
|
|
|
|
nel[2] = p3;
|
|
|
|
nel[3] = p4;
|
|
|
|
nel[4] = el[0] + el[1] + el[2] + el[3] - fixed[0] - moved[0] -
|
|
|
|
moved[1];
|
|
|
|
if (Cross(mesh[p2] - mesh[p1], mesh[p4] - mesh[p1]) *
|
|
|
|
(mesh[nel[4]] - mesh[nel[1]]) >
|
|
|
|
0)
|
|
|
|
Swap(nel[1], nel[3]);
|
|
|
|
nel.SetIndex(el.GetIndex());
|
|
|
|
mesh.AddVolumeElement(nel);
|
|
|
|
p1 = p4;
|
|
|
|
p2 = p3;
|
|
|
|
}
|
2020-11-24 03:48:49 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
|
|
|
if (moved.Size() == 1 && fixed.Size() == 1) {
|
|
|
|
PointIndex p1 = moved[0];
|
2024-09-27 16:12:14 +05:00
|
|
|
for (auto i : Range(par_heights)) {
|
2024-04-05 20:09:34 +05:00
|
|
|
Element nel = el;
|
|
|
|
PointIndex p2 = newPoint(moved[0], i);
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto &p : nel.PNums()) {
|
2024-04-05 20:09:34 +05:00
|
|
|
if (p == moved[0])
|
|
|
|
p = p1;
|
|
|
|
else if (p == fixed[0])
|
|
|
|
p = p2;
|
|
|
|
}
|
|
|
|
p1 = p2;
|
|
|
|
mesh.AddVolumeElement(nel);
|
2020-11-17 22:43:39 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
|
|
|
} else if (el.GetType() == PYRAMID) {
|
|
|
|
if (moved.Size() == 2) {
|
|
|
|
if (fixed.Size() != 2)
|
|
|
|
throw Exception("This case is not implemented yet! Fixed size = " +
|
|
|
|
ToString(fixed.Size()));
|
|
|
|
PointIndex p1 = moved[0];
|
|
|
|
PointIndex p2 = moved[1];
|
2024-09-27 16:12:14 +05:00
|
|
|
for (auto i : Range(par_heights)) {
|
2024-04-05 20:09:34 +05:00
|
|
|
PointIndex p3 = newPoint(moved[1], i);
|
|
|
|
PointIndex p4 = newPoint(moved[0], i);
|
|
|
|
Element nel(PYRAMID);
|
|
|
|
nel[0] = p1;
|
|
|
|
nel[1] = p2;
|
|
|
|
nel[2] = p3;
|
|
|
|
nel[3] = p4;
|
|
|
|
nel[4] = el[0] + el[1] + el[2] + el[3] + el[4] - fixed[0] -
|
|
|
|
fixed[1] - moved[0] - moved[1];
|
|
|
|
if (Cross(mesh[p2] - mesh[p1], mesh[p4] - mesh[p1]) *
|
|
|
|
(mesh[nel[4]] - mesh[nel[1]]) >
|
|
|
|
0)
|
|
|
|
Swap(nel[1], nel[3]);
|
|
|
|
nel.SetIndex(el.GetIndex());
|
|
|
|
mesh.AddVolumeElement(nel);
|
|
|
|
p1 = p4;
|
|
|
|
p2 = p3;
|
2020-11-24 03:48:49 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
} else if (moved.Size() == 1)
|
|
|
|
throw Exception("This case is not implemented yet!");
|
|
|
|
} else if (do_move) {
|
|
|
|
throw Exception(
|
|
|
|
"Boundarylayer only implemented for tets and pyramids outside "
|
|
|
|
"yet!");
|
2020-11-17 22:43:39 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2022-03-07 19:58:09 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2009-06-19 11:43:23 +06:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
void BoundaryLayerTool ::SetDomInOut() {
|
|
|
|
for (auto i : Range(1, nfd_old + 1))
|
|
|
|
if (moved_surfaces.Test(i)) {
|
|
|
|
if (auto dom = mesh.GetFaceDescriptor(si_map[i]).DomainIn();
|
|
|
|
dom > ndom_old)
|
|
|
|
mesh.GetFaceDescriptor(i).SetDomainOut(dom);
|
|
|
|
else
|
|
|
|
mesh.GetFaceDescriptor(i).SetDomainIn(
|
|
|
|
mesh.GetFaceDescriptor(si_map[i]).DomainOut());
|
|
|
|
}
|
|
|
|
}
|
2023-03-31 02:21:58 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
void BoundaryLayerTool ::SetDomInOutSides() {
|
|
|
|
BitArray done(mesh.GetNFD() + 1);
|
|
|
|
done.Clear();
|
|
|
|
for (auto sei : Range(mesh.SurfaceElements())) {
|
2024-09-27 20:30:36 +05:00
|
|
|
auto &sel = mesh[sei];
|
2024-04-05 20:09:34 +05:00
|
|
|
auto index = sel.GetIndex();
|
2024-09-27 20:30:36 +05:00
|
|
|
if (done.Test(index))
|
|
|
|
continue;
|
2024-04-05 20:09:34 +05:00
|
|
|
done.SetBit(index);
|
2024-09-27 20:30:36 +05:00
|
|
|
auto &fd = mesh.GetFaceDescriptor(index);
|
|
|
|
if (fd.DomainIn() != -1)
|
|
|
|
continue;
|
2024-04-05 20:09:34 +05:00
|
|
|
int e1, e2;
|
|
|
|
mesh.GetTopology().GetSurface2VolumeElement(sei + 1, e1, e2);
|
|
|
|
if (e1 == 0)
|
|
|
|
fd.SetDomainIn(0);
|
2022-02-16 23:51:54 +05:00
|
|
|
else
|
2024-04-05 20:09:34 +05:00
|
|
|
fd.SetDomainIn(mesh.VolumeElement(e1).GetIndex());
|
|
|
|
if (e2 == 0)
|
|
|
|
fd.SetDomainOut(0);
|
|
|
|
else
|
|
|
|
fd.SetDomainOut(mesh.VolumeElement(e2).GetIndex());
|
2022-03-07 19:58:09 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2022-03-07 20:04:21 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
void BoundaryLayerTool ::AddSegments() {
|
|
|
|
if (have_single_segments)
|
|
|
|
MergeAndAddSegments(mesh, segments, new_segments);
|
|
|
|
else {
|
|
|
|
mesh.LineSegments() = segments;
|
2024-09-27 20:30:36 +05:00
|
|
|
for (auto &seg : new_segments)
|
|
|
|
mesh.AddSegment(seg);
|
2022-03-07 19:58:09 +05:00
|
|
|
}
|
2024-04-05 20:09:34 +05:00
|
|
|
}
|
2022-03-07 19:58:09 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
void BoundaryLayerTool ::FixVolumeElements() {
|
|
|
|
static Timer timer("BoundaryLayerTool::FixVolumeElements");
|
|
|
|
RegionTimer rt(timer);
|
|
|
|
BitArray is_inner_point(mesh.GetNP() + 1);
|
|
|
|
is_inner_point.Clear();
|
|
|
|
|
|
|
|
auto changed_domains = domains;
|
2024-09-27 20:30:36 +05:00
|
|
|
if (!params.outside)
|
|
|
|
changed_domains.Invert();
|
2024-04-05 20:09:34 +05:00
|
|
|
|
|
|
|
for (ElementIndex ei : Range(ne))
|
|
|
|
if (changed_domains.Test(mesh[ei].GetIndex()))
|
|
|
|
for (auto pi : mesh[ei].PNums())
|
2024-09-27 20:30:36 +05:00
|
|
|
if (mesh[pi].Type() == INNERPOINT)
|
|
|
|
is_inner_point.SetBit(pi);
|
2024-04-05 20:09:34 +05:00
|
|
|
|
|
|
|
Array<PointIndex> points;
|
|
|
|
for (auto pi : mesh.Points().Range())
|
2024-09-27 20:30:36 +05:00
|
|
|
if (is_inner_point.Test(pi))
|
|
|
|
points.Append(pi);
|
2024-04-05 20:09:34 +05:00
|
|
|
|
|
|
|
auto p2el = mesh.CreatePoint2ElementTable(is_inner_point);
|
|
|
|
|
|
|
|
// smooth growth vectors to shift additional element layers to the inside and
|
|
|
|
// fix flipped tets
|
|
|
|
for ([[maybe_unused]] auto step : Range(0)) {
|
|
|
|
for (auto pi : points) {
|
|
|
|
Vec<3> average_gw = 0.0;
|
2024-09-27 20:30:36 +05:00
|
|
|
auto &els = p2el[pi];
|
2024-04-05 20:09:34 +05:00
|
|
|
size_t cnt = 0;
|
|
|
|
for (auto ei : els)
|
|
|
|
if (ei < ne)
|
|
|
|
for (auto pi1 : mesh[ei].PNums())
|
|
|
|
if (pi1 <= np) {
|
|
|
|
average_gw += growthvectors[pi1];
|
|
|
|
cnt++;
|
|
|
|
}
|
|
|
|
growthvectors[pi] = 1.0 / cnt * average_gw;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2023-09-04 14:48:15 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
void BoundaryLayerTool ::ProcessParameters() {
|
|
|
|
if (int *bc = get_if<int>(¶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<string>(¶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<string>(¶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);
|
2024-09-27 16:12:14 +05:00
|
|
|
}
|
2024-09-27 20:30:36 +05:00
|
|
|
// 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>>(¶ms.boundary);
|
|
|
|
for (auto id : surfids)
|
|
|
|
par_surfid.Append(id);
|
|
|
|
}
|
|
|
|
if (string *mat = get_if<string>(¶ms.new_material); mat)
|
|
|
|
par_new_mat = {{".*", *mat}};
|
|
|
|
else
|
|
|
|
par_new_mat = *get_if<map<string, string>>(¶ms.new_material);
|
|
|
|
|
|
|
|
if (params.project_boundaries.has_value()) {
|
|
|
|
auto proj_bnd = *params.project_boundaries;
|
|
|
|
if (string *s = get_if<string>(&proj_bnd); s) {
|
|
|
|
regex pattern(*s);
|
|
|
|
for (int i = 1; i <= mesh.GetNFD(); i++)
|
|
|
|
if (regex_match(mesh.GetFaceDescriptor(i).GetBCName(), pattern))
|
|
|
|
par_project_boundaries.Append(i);
|
|
|
|
} else {
|
|
|
|
for (auto id : *get_if<std::vector<int>>(&proj_bnd))
|
|
|
|
par_project_boundaries.Append(id);
|
|
|
|
}
|
|
|
|
}
|
2024-09-27 16:12:14 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
if (double *height = get_if<double>(¶ms.thickness); height) {
|
|
|
|
par_heights.Append(*height);
|
|
|
|
} else {
|
|
|
|
auto &heights = *get_if<std::vector<double>>(¶ms.thickness);
|
|
|
|
for (auto val : heights)
|
|
|
|
par_heights.Append(val);
|
|
|
|
}
|
2024-09-27 16:12:14 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
int nr_domains = mesh.GetNDomains();
|
|
|
|
domains.SetSize(nr_domains + 1); // one based
|
|
|
|
domains.Clear();
|
|
|
|
if (string *pdomain = get_if<string>(¶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 if (int *idomain = get_if<int>(¶ms.domain); idomain) {
|
|
|
|
domains.SetBit(*idomain);
|
|
|
|
} else {
|
|
|
|
for (auto i : *get_if<std::vector<int>>(¶ms.domain))
|
|
|
|
domains.SetBit(i);
|
2024-09-27 16:12:14 +05:00
|
|
|
}
|
2024-09-27 20:30:36 +05:00
|
|
|
}
|
2024-09-27 16:12:14 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
void BoundaryLayerTool ::Perform() {
|
|
|
|
CreateNewFaceDescriptors();
|
|
|
|
CalculateGrowthVectors();
|
|
|
|
CreateFaceDescriptorsSides();
|
|
|
|
auto segmap = BuildSegMap();
|
2023-12-11 23:05:04 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
auto in_surface_direction = ProjectGrowthVectorsOnSurface();
|
2024-03-01 16:01:08 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
InsertNewElements(segmap, in_surface_direction);
|
2023-11-29 13:19:27 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
SetDomInOut();
|
|
|
|
AddSegments();
|
2023-11-29 13:19:27 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
mesh.CalcSurfacesOfNode();
|
|
|
|
topo.SetBuildVertex2Element(true);
|
|
|
|
mesh.UpdateTopology();
|
|
|
|
InterpolateGrowthVectors();
|
2023-11-29 13:19:27 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
if (params.limit_growth_vectors)
|
|
|
|
LimitGrowthVectorLengths();
|
2024-02-20 14:25:31 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
for (auto [pi, data] : growth_vector_map) {
|
|
|
|
auto [gw, height] = data;
|
|
|
|
mesh[pi] += height * (*gw);
|
2022-03-07 19:58:09 +05:00
|
|
|
}
|
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
mesh.GetTopology().ClearEdges();
|
|
|
|
mesh.SetNextMajorTimeStamp();
|
|
|
|
mesh.UpdateTopology();
|
|
|
|
SetDomInOutSides();
|
|
|
|
}
|
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
void GenerateBoundaryLayer(Mesh &mesh, const BoundaryLayerParameters &blp) {
|
2024-04-05 20:09:34 +05:00
|
|
|
static Timer timer("Create Boundarylayers");
|
|
|
|
RegionTimer regt(timer);
|
2022-03-07 19:58:09 +05:00
|
|
|
|
2024-04-05 20:09:34 +05:00
|
|
|
BoundaryLayerTool tool(mesh, blp);
|
|
|
|
tool.Perform();
|
|
|
|
}
|
2021-03-16 22:20:40 +05:00
|
|
|
|
2024-09-27 20:30:36 +05:00
|
|
|
} // namespace netgen
|