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

formatting

Browse Source
This commit is contained in:
Matthias Hochsteger 2024-10-08 16:20:53 +02:00
parent e5af7bca42
commit 82965f63b0
5 changed files with 1709 additions and 1437 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

64
.clang-format Normal file
View File

@ -0,0 +1,64 @@
Language: Cpp
BasedOnStyle: LLVM
AlignAfterOpenBracket: Align
AlignConsecutiveAssignments: false
AlignConsecutiveDeclarations: false
AlignEscapedNewlines: Left
AlignOperands: true
AlignTrailingComments: true
AllowAllParametersOfDeclarationOnNextLine: false
AllowShortBlocksOnASingleLine: false
AllowShortCaseLabelsOnASingleLine: true
AllowShortFunctionsOnASingleLine: InlineOnly
AllowShortIfStatementsOnASingleLine: false
AllowShortLoopsOnASingleLine: false
AlwaysBreakAfterDefinitionReturnType: None
AlwaysBreakTemplateDeclarations: Yes
BinPackArguments: false
BinPackParameters: false
BreakBeforeBinaryOperators: NonAssignment
BreakBeforeBraces: Custom
BraceWrapping:
AfterClass: true
AfterControlStatement: true
AfterEnum: true
AfterFunction: true
AfterNamespace: true
AfterStruct: true
AfterUnion: true
BeforeCatch: true
BeforeElse: true
IndentBraces: true
BreakBeforeTernaryOperators: true
BreakConstructorInitializersBeforeComma: false
BreakInheritanceList: AfterColon
ColumnLimit: 0
ConstructorInitializerAllOnOneLineOrOnePerLine: true
ConstructorInitializerIndentWidth: 2
ContinuationIndentWidth: 2
Cpp11BracedListStyle: true
EmptyLineBeforeAccessModifier: Never
IndentCaseLabels: false
IndentPPDirectives: None
IndentWidth: 2
KeepEmptyLinesAtTheStartOfBlocks: false
MaxEmptyLinesToKeep: 1
NamespaceIndentation: None
PointerAlignment: Left
ReflowComments: true
SortIncludes: false
SpaceAfterCStyleCast: false
SpaceAfterLogicalNot: false
SpaceBeforeParens: Custom
SpaceBeforeParensOptions:
AfterControlStatements: true
AfterFunctionDefinitionName: true
AfterFunctionDeclarationName: true
SpacesInAngles: false
SpacesInContainerLiterals: false
SpacesInCStyleCastParentheses: false
SpacesInSquareBrackets: false
Standard: Latest
TabWidth: 2
UseTab: Never

1663
libsrc/meshing/boundarylayer.cpp
View File

File diff suppressed because it is too large Load Diff

127
libsrc/meshing/boundarylayer.hpp
View File

@ -9,101 +9,102 @@ namespace netgen
{ {
/// ///
DLL_HEADER extern void InsertVirtualBoundaryLayer (Mesh & mesh); DLL_HEADER extern void InsertVirtualBoundaryLayer (Mesh& mesh);
/// Create a typical prismatic boundary layer on the given /// Create a typical prismatic boundary layer on the given
/// surfaces /// surfaces
struct SpecialBoundaryPoint { struct SpecialBoundaryPoint
struct GrowthGroup { {
struct GrowthGroup
{
Array<int> faces; Array<int> faces;
Vec<3> growth_vector; Vec<3> growth_vector;
Array<PointIndex> new_points; Array<PointIndex> new_points;
GrowthGroup(FlatArray<int> faces_, FlatArray<Vec<3>> normals); GrowthGroup(FlatArray<int> faces_, FlatArray<Vec<3>> normals);
GrowthGroup(const GrowthGroup &) = default; GrowthGroup(const GrowthGroup&) = default;
GrowthGroup() = default; GrowthGroup() = default;
}; };
// std::map<int, Vec<3>> normals; // std::map<int, Vec<3>> normals;
Array<GrowthGroup> growth_groups; Array<GrowthGroup> growth_groups;
Vec<3> separating_direction; Vec<3> separating_direction;
SpecialBoundaryPoint( const std::map<int, Vec<3>> & normals ); SpecialBoundaryPoint(const std::map<int, Vec<3>>& normals);
SpecialBoundaryPoint() = default; SpecialBoundaryPoint() = default;
}; };
DLL_HEADER void GenerateBoundaryLayer (Mesh & mesh, DLL_HEADER void GenerateBoundaryLayer (Mesh& mesh,
const BoundaryLayerParameters & blp); const BoundaryLayerParameters& blp);
DLL_HEADER int /* new_domain_number */ GenerateBoundaryLayer2 (Mesh & mesh, int domain, const Array<double> & thicknesses, bool should_make_new_domain=true, const Array<int> & boundaries=Array<int>{}); DLL_HEADER int /* new_domain_number */ GenerateBoundaryLayer2 (Mesh& mesh, int domain, const Array<double>& thicknesses, bool should_make_new_domain = true, const Array<int>& boundaries = Array<int>{});
class BoundaryLayerTool class BoundaryLayerTool
{ {
public: public:
BoundaryLayerTool(Mesh & mesh_, const BoundaryLayerParameters & params_); BoundaryLayerTool(Mesh& mesh_, const BoundaryLayerParameters& params_);
void ProcessParameters(); void ProcessParameters ();
void Perform(); void Perform ();
Mesh & mesh; Mesh& mesh;
MeshTopology & topo; MeshTopology& topo;
BoundaryLayerParameters params; BoundaryLayerParameters params;
Array<Vec<3>, PointIndex> growthvectors; Array<Vec<3>, PointIndex> growthvectors;
std::map<PointIndex, Vec<3>> non_bl_growth_vectors; std::map<PointIndex, Vec<3>> non_bl_growth_vectors;
Table<SurfaceElementIndex, PointIndex> p2sel; Table<SurfaceElementIndex, PointIndex> p2sel;
BitArray domains, is_edge_moved, is_boundary_projected, is_boundary_moved; BitArray domains, is_edge_moved, is_boundary_projected, is_boundary_moved;
Array<SegmentIndex> moved_segs; Array<SegmentIndex> moved_segs;
int max_edge_nr, nfd_old, ndom_old; int max_edge_nr, nfd_old, ndom_old;
Array<int> new_mat_nrs; Array<int> new_mat_nrs;
BitArray moved_surfaces; BitArray moved_surfaces;
int np, nseg, nse, ne; int np, nseg, nse, ne;
double total_height; double total_height;
// These parameters are derived from given BoundaryLayerParameters and the Mesh // These parameters are derived from given BoundaryLayerParameters and the Mesh
Array<double> par_heights; Array<double> par_heights;
Array<int> par_surfid; Array<int> par_surfid;
bool insert_only_volume_elements; bool insert_only_volume_elements;
map<string, string> par_new_mat; map<string, string> par_new_mat;
Array<size_t> par_project_boundaries; Array<size_t> par_project_boundaries;
bool have_single_segments; bool have_single_segments;
Array<Segment> segments, new_segments, new_segments_on_moved_bnd; Array<Segment> segments, new_segments, new_segments_on_moved_bnd;
Array<Element2d, SurfaceElementIndex> new_sels, new_sels_on_moved_bnd; Array<Element2d, SurfaceElementIndex> new_sels, new_sels_on_moved_bnd;
Array<Array<PointIndex>, PointIndex> mapto; Array<Array<PointIndex>, PointIndex> mapto;
Array<PointIndex, PointIndex> mapfrom; Array<PointIndex, PointIndex> mapfrom;
Array<double> surfacefacs; Array<double> surfacefacs;
Array<int> si_map; Array<int> si_map;
std::map<PointIndex, SpecialBoundaryPoint> special_boundary_points; std::map<PointIndex, SpecialBoundaryPoint> special_boundary_points;
std::map<PointIndex, std::tuple<Vec<3>*, double>> growth_vector_map; std::map<PointIndex, std::tuple<Vec<3>*, double>> growth_vector_map;
// major steps called in Perform() // major steps called in Perform()
void CreateNewFaceDescriptors(); void CreateNewFaceDescriptors ();
void CreateFaceDescriptorsSides(); void CreateFaceDescriptorsSides ();
void CalculateGrowthVectors(); void CalculateGrowthVectors ();
Array<Array<pair<SegmentIndex, int>>, SegmentIndex> BuildSegMap(); Array<Array<pair<SegmentIndex, int>>, SegmentIndex> BuildSegMap ();
BitArray ProjectGrowthVectorsOnSurface(); BitArray ProjectGrowthVectorsOnSurface ();
void InterpolateSurfaceGrowthVectors(); void InterpolateSurfaceGrowthVectors ();
void InterpolateGrowthVectors(); void InterpolateGrowthVectors ();
void LimitGrowthVectorLengths(); void LimitGrowthVectorLengths ();
void FixSurfaceElements(); void FixSurfaceElements ();
void InsertNewElements(FlatArray<Array<pair<SegmentIndex, int>>, SegmentIndex> segmap, const BitArray & in_surface_direction); void InsertNewElements (FlatArray<Array<pair<SegmentIndex, int>>, SegmentIndex> segmap, const BitArray& in_surface_direction);
void SetDomInOut(); void SetDomInOut ();
void SetDomInOutSides(); void SetDomInOutSides ();
void AddSegments(); void AddSegments ();
void AddSurfaceElements(); void AddSurfaceElements ();
Vec<3> getNormal(const Element2d & el) Vec<3> getNormal (const Element2d& el)
{ {
auto v0 = mesh[el[0]]; auto v0 = mesh[el[0]];
return Cross(mesh[el[1]]-v0, mesh[el[2]]-v0).Normalize(); return Cross(mesh[el[1]] - v0, mesh[el[2]] - v0).Normalize();
} }
Vec<3> getEdgeTangent(PointIndex pi, int edgenr, FlatArray<Segment *> segs); Vec<3> getEdgeTangent (PointIndex pi, int edgenr, FlatArray<Segment*> segs);
}; };
} // namespace netgen } // namespace netgen

556
libsrc/meshing/boundarylayer_interpolate.cpp
View File

@ -1,9 +1,12 @@
#include "boundarylayer.hpp" #include "boundarylayer.hpp"
namespace netgen { namespace netgen
{
namespace detail { namespace detail
struct Neighbor { {
struct Neighbor
{
PointIndex pi; PointIndex pi;
SurfaceElementIndex sei; SurfaceElementIndex sei;
double weight; double weight;
@ -11,72 +14,79 @@ struct Neighbor {
} // namespace detail } // namespace detail
Array<ArrayMem<detail::Neighbor, 20>> Array<ArrayMem<detail::Neighbor, 20>>
BuildNeighbors(FlatArray<PointIndex> points, const Mesh &mesh) { BuildNeighbors (FlatArray<PointIndex> points, const Mesh& mesh)
{
auto p2sel = mesh.CreatePoint2SurfaceElementTable(); auto p2sel = mesh.CreatePoint2SurfaceElementTable();
Array<ArrayMem<detail::Neighbor, 20>> neighbors(points.Size()); Array<ArrayMem<detail::Neighbor, 20>> neighbors(points.Size());
ArrayMem<double, 20> angles; ArrayMem<double, 20> angles;
ArrayMem<double, 20> inv_dists; ArrayMem<double, 20> inv_dists;
for (auto i : points.Range()) { for (auto i : points.Range())
auto &p_neighbors = neighbors[i]; {
auto pi = points[i]; auto& p_neighbors = neighbors[i];
angles.SetSize(0); auto pi = points[i];
inv_dists.SetSize(0); angles.SetSize(0);
for (auto sei : p2sel[pi]) { inv_dists.SetSize(0);
const auto &sel = mesh[sei]; for (auto sei : p2sel[pi])
for (auto pi1 : sel.PNums()) { {
if (pi1 == pi) const auto& sel = mesh[sei];
continue; for (auto pi1 : sel.PNums())
auto pi2 = pi1; {
for (auto pi_ : sel.PNums()) { if (pi1 == pi)
if (pi_ != pi && pi_ != pi1) { continue;
pi2 = pi_; auto pi2 = pi1;
break; for (auto pi_ : sel.PNums())
} {
if (pi_ != pi && pi_ != pi1)
{
pi2 = pi_;
break;
}
}
p_neighbors.Append({pi1, sei, 0.0});
inv_dists.Append(1.0 / (mesh[pi1] - mesh[pi]).Length());
auto dot = (mesh[pi1] - mesh[pi]).Normalize() * (mesh[pi2] - mesh[pi]).Normalize();
angles.Append(acos(dot));
}
} }
p_neighbors.Append({pi1, sei, 0.0}); double sum_inv_dist = 0.0;
inv_dists.Append(1.0 / (mesh[pi1] - mesh[pi]).Length()); for (auto inv_dist : inv_dists)
auto dot = (mesh[pi1] - mesh[pi]).Normalize() * sum_inv_dist += inv_dist;
(mesh[pi2] - mesh[pi]).Normalize(); double sum_angle = 0.0;
angles.Append(acos(dot)); for (auto angle : angles)
} sum_angle += angle;
}
double sum_inv_dist = 0.0;
for (auto inv_dist : inv_dists)
sum_inv_dist += inv_dist;
double sum_angle = 0.0;
for (auto angle : angles)
sum_angle += angle;
double sum_weight = 0.0; double sum_weight = 0.0;
for (auto i : Range(inv_dists)) { for (auto i : Range(inv_dists))
p_neighbors[i].weight = {
inv_dists[i] * angles[i] / sum_inv_dist / sum_angle; p_neighbors[i].weight =
sum_weight += p_neighbors[i].weight; inv_dists[i] * angles[i] / sum_inv_dist / sum_angle;
sum_weight += p_neighbors[i].weight;
}
for (auto i : Range(inv_dists))
p_neighbors[i].weight /= sum_weight;
} }
for (auto i : Range(inv_dists))
p_neighbors[i].weight /= sum_weight;
}
return neighbors; return neighbors;
} }
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) if (seg.edgenr > new_max_edge_nr)
new_max_edge_nr = seg.edgenr; 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) if (seg.edgenr > new_max_edge_nr)
new_max_edge_nr = seg.edgenr; 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)
growth_vector_map[pi] = {&growthvectors[pi], total_height}; growth_vector_map[pi] = {&growthvectors[pi], total_height};
auto [gw, height] = growth_vector_map[pi]; auto [gw, height] = growth_vector_map[pi];
return height * (*gw); return height * (*gw);
}; };
auto addGW = [&](PointIndex pi, Vec<3> vec) { auto addGW = [&] (PointIndex pi, Vec<3> vec) {
if (growth_vector_map.count(pi) == 0) if (growth_vector_map.count(pi) == 0)
growth_vector_map[pi] = {&growthvectors[pi], total_height}; growth_vector_map[pi] = {&growthvectors[pi], total_height};
auto [gw, height] = growth_vector_map[pi]; auto [gw, height] = growth_vector_map[pi];
@ -87,129 +97,139 @@ void BoundaryLayerTool ::InterpolateGrowthVectors() {
if (max_edge_nr >= new_max_edge_nr) if (max_edge_nr >= new_max_edge_nr)
return; return;
auto edgenr2seg = ngcore::CreateSortedTable<Segment *, int>( auto edgenr2seg = ngcore::CreateSortedTable<Segment*, int>(
Range(segments.Size() + new_segments.Size()), Range(segments.Size() + new_segments.Size()),
[&](auto &table, size_t segi) { [&] (auto& table, size_t segi) {
auto &seg = segi < segments.Size() auto& seg = segi < segments.Size()
? segments[segi] ? segments[segi]
: new_segments[segi - segments.Size()]; : new_segments[segi - segments.Size()];
table.Add(seg.edgenr, &seg); table.Add(seg.edgenr, &seg);
}, },
new_max_edge_nr + 1); new_max_edge_nr + 1);
auto point2seg = ngcore::CreateSortedTable<Segment *, PointIndex>( auto point2seg = ngcore::CreateSortedTable<Segment*, PointIndex>(
Range(segments.Size() + new_segments.Size()), Range(segments.Size() + new_segments.Size()),
[&](auto &table, size_t segi) { [&] (auto& table, size_t segi) {
auto &seg = segi < segments.Size() auto& seg = segi < segments.Size()
? segments[segi] ? segments[segi]
: new_segments[segi - segments.Size()]; : new_segments[segi - segments.Size()];
table.Add(seg[0], &seg); table.Add(seg[0], &seg);
table.Add(seg[1], &seg); table.Add(seg[1], &seg);
}, },
mesh.GetNP()); mesh.GetNP());
for (auto edgenr : Range(max_edge_nr + 1, new_max_edge_nr + 1)) { for (auto edgenr : Range(max_edge_nr + 1, new_max_edge_nr + 1))
double edge_len = 0.; {
double edge_len = 0.;
auto is_end_point = [&](PointIndex pi) { auto is_end_point = [&] (PointIndex pi) {
auto segs = point2seg[pi]; auto segs = point2seg[pi];
if (segs.Size() == 1) if (segs.Size() == 1)
return true;
auto first_edgenr = (*segs[0]).edgenr;
for (auto *p_seg : segs)
if (p_seg->edgenr != first_edgenr)
return true; return true;
return false; auto first_edgenr = (*segs[0]).edgenr;
}; for (auto* p_seg : segs)
if (p_seg->edgenr != first_edgenr)
return true;
return false;
};
bool any_grows = false; bool any_grows = false;
Array<PointIndex> points; Array<PointIndex> points;
for (auto *p_seg : edgenr2seg[edgenr]) { for (auto* p_seg : edgenr2seg[edgenr])
auto &seg = *p_seg; {
auto& seg = *p_seg;
if (getGW(seg[0]).Length2() != 0 || getGW(seg[1]).Length2() != 0) if (getGW(seg[0]).Length2() != 0 || getGW(seg[1]).Length2() != 0)
any_grows = true; any_grows = true;
if (points.Size() == 0) if (points.Size() == 0)
for (auto i : Range(2)) for (auto i : Range(2))
if (is_end_point(seg[i])) { if (is_end_point(seg[i]))
points.Append(seg[i]); {
points.Append(seg[1 - i]); points.Append(seg[i]);
edge_len += (mesh[seg[1]] - mesh[seg[0]]).Length(); points.Append(seg[1 - i]);
break; edge_len += (mesh[seg[1]] - mesh[seg[0]]).Length();
} break;
} }
if (!any_grows) {
PrintMessage(1, "BLayer: skip interpolating growth vectors at edge ",
edgenr + 1);
continue;
}
if (!points.Size()) {
cerr << "Could not find startpoint for edge " << edgenr << endl;
continue;
}
std::set<PointIndex> points_set;
points_set.insert(points[0]);
points_set.insert(points[1]);
bool point_found = true;
while (point_found) {
if (is_end_point(points.Last()))
break;
point_found = false;
for (auto *p_seg : point2seg[points.Last()]) {
const auto &seg = *p_seg;
if (seg.edgenr != edgenr)
continue;
auto plast = points.Last();
if (plast != seg[0] && plast != seg[1])
continue;
auto pnew = plast == seg[0] ? seg[1] : seg[0];
if (pnew == points[0] && points.Size() > 1) {
} }
if (points_set.count(pnew) > 0 &&
(pnew != points[0] || points.Size() == 2)) if (!any_grows)
{
PrintMessage(1, "BLayer: skip interpolating growth vectors at edge ", edgenr + 1);
continue; continue;
edge_len += (mesh[points.Last()] - mesh[pnew]).Length(); }
points.Append(pnew);
points_set.insert(pnew);
point_found = true;
break;
}
}
if (!point_found) {
cerr << "Could not find connected list of line segments for edge "
<< edgenr << endl;
cerr << "current points: " << endl << points << endl;
continue;
}
if (getGW(points[0]).Length2() == 0 && getGW(points.Last()).Length2() == 0) if (!points.Size())
continue; {
cerr << "Could not find startpoint for edge " << edgenr << endl;
continue;
}
// tangential part of growth vectors std::set<PointIndex> points_set;
auto t1 = (mesh[points[1]] - mesh[points[0]]).Normalize(); points_set.insert(points[0]);
auto gt1 = getGW(points[0]) * t1 * t1; points_set.insert(points[1]);
auto t2 =
bool point_found = true;
while (point_found)
{
if (is_end_point(points.Last()))
break;
point_found = false;
for (auto* p_seg : point2seg[points.Last()])
{
const auto& seg = *p_seg;
if (seg.edgenr != edgenr)
continue;
auto plast = points.Last();
if (plast != seg[0] && plast != seg[1])
continue;
auto pnew = plast == seg[0] ? seg[1] : seg[0];
if (pnew == points[0] && points.Size() > 1)
{
}
if (points_set.count(pnew) > 0 && (pnew != points[0] || points.Size() == 2))
continue;
edge_len += (mesh[points.Last()] - mesh[pnew]).Length();
points.Append(pnew);
points_set.insert(pnew);
point_found = true;
break;
}
}
if (!point_found)
{
cerr << "Could not find connected list of line segments for edge "
<< edgenr << endl;
cerr << "current points: " << endl
<< points << endl;
continue;
}
if (getGW(points[0]).Length2() == 0 && getGW(points.Last()).Length2() == 0)
continue;
// 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(); (mesh[points.Last()] - mesh[points[points.Size() - 2]]).Normalize();
auto gt2 = getGW(points.Last()) * t2 * t2; auto gt2 = getGW(points.Last()) * t2 * t2;
double len = 0.; double len = 0.;
for (auto i : IntRange(1, points.Size() - 1)) { for (auto i : IntRange(1, points.Size() - 1))
auto pi = points[i]; {
len += (mesh[pi] - mesh[points[i - 1]]).Length(); auto pi = points[i];
auto t = getEdgeTangent(pi, edgenr, point2seg[pi]); len += (mesh[pi] - mesh[points[i - 1]]).Length();
auto lam = len / edge_len; auto t = getEdgeTangent(pi, edgenr, point2seg[pi]);
auto interpol = (1 - lam) * (gt1 * t) * t + lam * (gt2 * t) * t; auto lam = len / edge_len;
addGW(pi, interpol); auto interpol = (1 - lam) * (gt1 * t) * t + lam * (gt2 * t) * t;
addGW(pi, interpol);
}
} }
}
} }
void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() { void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors()
{
static Timer tall("InterpolateSurfaceGrowthVectors"); static Timer tall("InterpolateSurfaceGrowthVectors");
RegionTimer rtall(tall); RegionTimer rtall(tall);
static Timer tsmooth("InterpolateSurfaceGrowthVectors-Smoothing"); static Timer tsmooth("InterpolateSurfaceGrowthVectors-Smoothing");
@ -225,7 +245,7 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
// } // }
// cout << __FILE__ << ":" << __LINE__ << endl; // cout << __FILE__ << ":" << __LINE__ << endl;
auto getGW = [&](PointIndex pi) -> Vec<3> { auto getGW = [&] (PointIndex pi) -> Vec<3> {
return growthvectors[pi]; return growthvectors[pi];
// if (growth_vector_map.count(pi) == 0) { // if (growth_vector_map.count(pi) == 0) {
// non_bl_growth_vectors[pi] = .0; // non_bl_growth_vectors[pi] = .0;
@ -234,7 +254,7 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
// auto [gw, height] = growth_vector_map[pi]; // auto [gw, height] = growth_vector_map[pi];
// return height * (*gw); // return height * (*gw);
}; };
auto addGW = [&](PointIndex pi, Vec<3> vec) { auto addGW = [&] (PointIndex pi, Vec<3> vec) {
growthvectors[pi] += vec; growthvectors[pi] += vec;
// // cout << "add gw " << pi << "\t" << vec << endl; // // cout << "add gw " << pi << "\t" << vec << endl;
// if (growth_vector_map.count(pi) == 0) { // if (growth_vector_map.count(pi) == 0) {
@ -247,18 +267,18 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
// *gw += 1.0 / height * vec; // *gw += 1.0 / height * vec;
}; };
auto hasMoved = [&](PointIndex pi) { auto hasMoved = [&] (PointIndex pi) {
return (pi - PointIndex::BASE >= np_old) || mapto[pi].Size() > 0 || return (pi - PointIndex::BASE >= np_old) || mapto[pi].Size() > 0 || special_boundary_points.count(pi);
special_boundary_points.count(pi);
}; };
// cout << __FILE__ << ":" << __LINE__ << endl; // cout << __FILE__ << ":" << __LINE__ << endl;
std::set<PointIndex> points_set; std::set<PointIndex> points_set;
for (const auto &sel : mesh.SurfaceElements()) { for (const auto& sel : mesh.SurfaceElements())
for (auto pi : sel.PNums()) {
if (mesh[pi].Type() == SURFACEPOINT && hasMoved(pi)) for (auto pi : sel.PNums())
points_set.insert(pi); if (mesh[pi].Type() == SURFACEPOINT && hasMoved(pi))
} points_set.insert(pi);
}
Array<PointIndex> points; Array<PointIndex> points;
for (auto pi : points_set) for (auto pi : points_set)
@ -277,70 +297,76 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
BitArray interpolate_tangent(mesh.GetNP() + 1); BitArray interpolate_tangent(mesh.GetNP() + 1);
interpolate_tangent = false; interpolate_tangent = false;
for (auto pi : points) { for (auto pi : points)
for (auto sei : p2sel[pi]) {
if (is_boundary_moved[mesh[sei].GetIndex()]) for (auto sei : p2sel[pi])
interpolate_tangent.SetBit(pi); if (is_boundary_moved[mesh[sei].GetIndex()])
} interpolate_tangent.SetBit(pi);
}
constexpr int N_STEPS = 64; constexpr int N_STEPS = 64;
for ([[maybe_unused]] auto i : Range(N_STEPS)) { for ([[maybe_unused]] auto i : Range(N_STEPS))
for (auto i : points.Range()) { {
auto pi = points[i]; for (auto i : points.Range())
// cout << "AVERAGE " << pi << endl; {
auto &p_neighbors = neighbors[i]; auto pi = points[i];
// cout << "AVERAGE " << pi << endl;
auto& p_neighbors = neighbors[i];
ArrayMem<Vec<3>, 20> g_vectors; ArrayMem<Vec<3>, 20> g_vectors;
double max_len = 0.0; double max_len = 0.0;
double sum_len = 0.0; double sum_len = 0.0;
// average only tangent component on new bl points, average whole growth // average only tangent component on new bl points, average whole growth
// vector otherwise // vector otherwise
bool do_average_tangent = true; bool do_average_tangent = true;
for (const auto &s : p_neighbors) { for (const auto& s : p_neighbors)
auto gw_other = getGW(s.pi) + corrections[s.pi]; {
if (do_average_tangent) { auto gw_other = getGW(s.pi) + corrections[s.pi];
auto n = surf_normals[s.sei]; if (do_average_tangent)
gw_other = gw_other - (gw_other * n) * n; {
auto n = surf_normals[s.sei];
gw_other = gw_other - (gw_other * n) * n;
}
auto v = gw_other;
auto len = v.Length2();
sum_len += len;
max_len = max(max_len, len);
g_vectors.Append(v);
}
if (max_len == 0.0)
continue;
double lambda = 0;
if (i > N_STEPS / 4.)
lambda = 2.0 * (i - N_STEPS / 4.) / (N_STEPS / 2.);
lambda = min(1.0, lambda);
auto& correction = corrections[pi];
correction = 0.0;
for (const auto i : p_neighbors.Range())
{
auto v = g_vectors[i];
double weight = lambda * p_neighbors[i].weight + (1.0 - lambda) * v.Length2() / sum_len;
// if(pi == 19911) cout << "pi " << pi << "\tneighbor " <<
// p_neighbors[i].pi << "\tweight " << weight << endl;
correction += weight * v;
}
if (!do_average_tangent)
correction -= getGW(pi);
// if(pi == 19911) cout << "pi " << pi << "\tcorrection " << correction <<
// endl;
} }
auto v = gw_other;
auto len = v.Length2();
sum_len += len;
max_len = max(max_len, len);
g_vectors.Append(v);
}
if (max_len == 0.0)
continue;
double lambda = 0;
if (i > N_STEPS / 4.)
lambda = 2.0 * (i - N_STEPS / 4.) / (N_STEPS / 2.);
lambda = min(1.0, lambda);
auto &correction = corrections[pi];
correction = 0.0;
for (const auto i : p_neighbors.Range()) {
auto v = g_vectors[i];
double weight = lambda * p_neighbors[i].weight +
(1.0 - lambda) * v.Length2() / sum_len;
// if(pi == 19911) cout << "pi " << pi << "\tneighbor " <<
// p_neighbors[i].pi << "\tweight " << weight << endl;
correction += weight * v;
}
if (!do_average_tangent)
correction -= getGW(pi);
// if(pi == 19911) cout << "pi " << pi << "\tcorrection " << correction <<
// endl;
} }
}
for (auto pi : points) for (auto pi : points)
addGW(pi, corrections[pi]); addGW(pi, corrections[pi]);
} }
void BoundaryLayerTool ::FixSurfaceElements() { void BoundaryLayerTool ::FixSurfaceElements()
{
static Timer tall("FixSurfaceElements"); static Timer tall("FixSurfaceElements");
RegionTimer rtall(tall); RegionTimer rtall(tall);
auto np_old = this->np; auto np_old = this->np;
@ -348,34 +374,37 @@ void BoundaryLayerTool ::FixSurfaceElements() {
non_bl_growth_vectors.clear(); non_bl_growth_vectors.clear();
auto getGW = [&](PointIndex pi) -> Vec<3> { auto getGW = [&] (PointIndex pi) -> Vec<3> {
// return growthvectors[pi]; // return growthvectors[pi];
if (growth_vector_map.count(pi) == 0) { if (growth_vector_map.count(pi) == 0)
non_bl_growth_vectors[pi] = .0; {
growth_vector_map[pi] = {&non_bl_growth_vectors[pi], 1.0}; non_bl_growth_vectors[pi] = .0;
} growth_vector_map[pi] = {&non_bl_growth_vectors[pi], 1.0};
}
auto [gw, height] = growth_vector_map[pi]; auto [gw, height] = growth_vector_map[pi];
return height * (*gw); return height * (*gw);
}; };
auto addGW = [&](PointIndex pi, Vec<3> vec) { auto addGW = [&] (PointIndex pi, Vec<3> vec) {
if (growth_vector_map.count(pi) == 0) { if (growth_vector_map.count(pi) == 0)
non_bl_growth_vectors[pi] = .0; {
growth_vector_map[pi] = {&non_bl_growth_vectors[pi], 1.0}; non_bl_growth_vectors[pi] = .0;
} growth_vector_map[pi] = {&non_bl_growth_vectors[pi], 1.0};
}
auto [gw, height] = growth_vector_map[pi]; auto [gw, height] = growth_vector_map[pi];
*gw += 1.0 / height * vec; *gw += 1.0 / height * vec;
}; };
std::set<PointIndex> points_set; std::set<PointIndex> points_set;
// only smooth over old surface elements // only smooth over old surface elements
for (SurfaceElementIndex sei : Range(nse)) { for (SurfaceElementIndex sei : Range(nse))
const auto &sel = mesh[sei]; {
if (sel.GetNP() == 3 && is_boundary_moved[sel.GetIndex()]) const auto& sel = mesh[sei];
for (auto pi : sel.PNums()) if (sel.GetNP() == 3 && is_boundary_moved[sel.GetIndex()])
if (mesh[pi].Type() == SURFACEPOINT) for (auto pi : sel.PNums())
points_set.insert(pi); if (mesh[pi].Type() == SURFACEPOINT)
} points_set.insert(pi);
}
Array<PointIndex> points; Array<PointIndex> points;
for (auto pi : points_set) for (auto pi : points_set)
@ -388,41 +417,44 @@ void BoundaryLayerTool ::FixSurfaceElements() {
auto neighbors = BuildNeighbors(points, mesh); auto neighbors = BuildNeighbors(points, mesh);
constexpr int N_STEPS = 32; constexpr int N_STEPS = 32;
for ([[maybe_unused]] auto i : Range(N_STEPS)) { for ([[maybe_unused]] auto i : Range(N_STEPS))
for (auto i : points.Range()) { {
auto pi = points[i]; for (auto i : points.Range())
auto &p_neighbors = neighbors[i]; {
auto pi = points[i];
auto& p_neighbors = neighbors[i];
ArrayMem<Vec<3>, 20> g_vectors; ArrayMem<Vec<3>, 20> g_vectors;
double max_len = 0.0; double max_len = 0.0;
double sum_len = 0.0; double sum_len = 0.0;
for (const auto &s : p_neighbors) { for (const auto& s : p_neighbors)
auto v = getGW(s.pi) + corrections[s.pi]; {
auto len = v.Length2(); auto v = getGW(s.pi) + corrections[s.pi];
sum_len += len; auto len = v.Length2();
max_len = max(max_len, len); sum_len += len;
g_vectors.Append(v); max_len = max(max_len, len);
} g_vectors.Append(v);
}
if (max_len == 0.0) if (max_len == 0.0)
continue; continue;
double lambda = 0; double lambda = 0;
if (i > N_STEPS / 4.) if (i > N_STEPS / 4.)
lambda = 2.0 * (i - N_STEPS / 4.) / (N_STEPS / 2.); lambda = 2.0 * (i - N_STEPS / 4.) / (N_STEPS / 2.);
lambda = min(1.0, lambda); lambda = min(1.0, lambda);
auto &correction = corrections[pi]; auto& correction = corrections[pi];
correction = 0.0; correction = 0.0;
for (const auto i : p_neighbors.Range()) { for (const auto i : p_neighbors.Range())
auto v = g_vectors[i]; {
double weight = lambda * p_neighbors[i].weight + auto v = g_vectors[i];
(1.0 - lambda) * v.Length2() / sum_len; double weight = lambda * p_neighbors[i].weight + (1.0 - lambda) * v.Length2() / sum_len;
correction += weight * v; correction += weight * v;
} }
}
} }
}
for (auto pi : points) for (auto pi : points)
addGW(pi, corrections[pi]); addGW(pi, corrections[pi]);

736
libsrc/meshing/boundarylayer_limiter.hpp
View File

@ -1,9 +1,11 @@
#include "boundarylayer.hpp" #include "boundarylayer.hpp"
#include <core/array.hpp> #include <core/array.hpp>
namespace netgen { namespace netgen
{
struct Intersection_ { struct Intersection_
{
bool is_intersecting = false; bool is_intersecting = false;
double lam0 = -1, lam1 = -1; double lam0 = -1, lam1 = -1;
Point<3> p; Point<3> p;
@ -11,13 +13,14 @@ struct Intersection_ {
operator bool() const { return is_intersecting; } operator bool() const { return is_intersecting; }
}; };
struct GrowthVectorLimiter { struct GrowthVectorLimiter
{
typedef std::array<Point<3>, 2> Seg; typedef std::array<Point<3>, 2> Seg;
typedef std::array<Point<3>, 3> Trig; typedef std::array<Point<3>, 3> Trig;
BoundaryLayerTool &tool; BoundaryLayerTool& tool;
const BoundaryLayerParameters &params; const BoundaryLayerParameters& params;
Mesh &mesh; Mesh& mesh;
double height; double height;
Array<double, PointIndex> limits; Array<double, PointIndex> limits;
FlatArray<Vec<3>, PointIndex> growthvectors; FlatArray<Vec<3>, PointIndex> growthvectors;
@ -26,94 +29,102 @@ struct GrowthVectorLimiter {
Array<PointIndex, PointIndex> map_from; Array<PointIndex, PointIndex> map_from;
Table<SurfaceElementIndex, PointIndex> p2sel; Table<SurfaceElementIndex, PointIndex> p2sel;
GrowthVectorLimiter(BoundaryLayerTool &tool_) GrowthVectorLimiter(BoundaryLayerTool& tool_)
: tool(tool_), params(tool_.params), mesh(tool_.mesh), : tool(tool_), params(tool_.params), mesh(tool_.mesh), height(tool_.total_height), growthvectors(tool_.growthvectors), map_from(mesh.Points().Size())
height(tool_.total_height), growthvectors(tool_.growthvectors), {
map_from(mesh.Points().Size()) {
changed_domains = tool.domains; changed_domains = tool.domains;
if (!params.outside) if (!params.outside)
changed_domains.Invert(); changed_domains.Invert();
map_from = tool.mapfrom; map_from = tool.mapfrom;
p2sel = ngcore::CreateSortedTable<SurfaceElementIndex, PointIndex>( p2sel = ngcore::CreateSortedTable<SurfaceElementIndex, PointIndex>(
tool.new_sels.Range(), tool.new_sels.Range(),
[&](auto &table, SurfaceElementIndex ei) { [&] (auto& table, SurfaceElementIndex ei) {
for (PointIndex pi : tool.new_sels[ei].PNums()) for (PointIndex pi : tool.new_sels[ei].PNums())
table.Add(pi, ei); table.Add(pi, ei);
}, },
mesh.GetNP()); mesh.GetNP());
} }
auto SurfaceElementsRange() { return Range(tool.nse + tool.new_sels.Size()); } auto SurfaceElementsRange () { return Range(tool.nse + tool.new_sels.Size()); }
const auto &Get(SurfaceElementIndex sei) { const auto& Get (SurfaceElementIndex sei)
{
if (sei < tool.nse) if (sei < tool.nse)
return mesh[sei]; return mesh[sei];
return tool.new_sels[sei - tool.nse]; return tool.new_sels[sei - tool.nse];
} }
std::pair<double, double> GetMinMaxLimit(SurfaceElementIndex sei) { std::pair<double, double> GetMinMaxLimit (SurfaceElementIndex sei)
const auto &sel = Get(sei); {
const auto& sel = Get(sei);
double min_limit = GetLimit(sel[0]); double min_limit = GetLimit(sel[0]);
double max_limit = min_limit; double max_limit = min_limit;
for (auto i : IntRange(1, sel.GetNP())) { for (auto i : IntRange(1, sel.GetNP()))
auto limit = GetLimit(sel[i]); {
min_limit = min(min_limit, limit); auto limit = GetLimit(sel[i]);
max_limit = max(max_limit, limit); min_limit = min(min_limit, limit);
} max_limit = max(max_limit, limit);
}
return {min_limit, max_limit}; return {min_limit, max_limit};
} }
double GetLimit(PointIndex pi) { double GetLimit (PointIndex pi)
{
if (pi <= tool.np) if (pi <= tool.np)
return limits[pi]; return limits[pi];
return limits[map_from[pi]]; return limits[map_from[pi]];
} }
bool SetLimit(PointIndex pi, double new_limit) { bool SetLimit (PointIndex pi, double new_limit)
double &limit = (pi <= tool.np) ? limits[pi] : limits[map_from[pi]]; {
double& limit = (pi <= tool.np) ? limits[pi] : limits[map_from[pi]];
if (limit <= new_limit) if (limit <= new_limit)
return false; return false;
limit = new_limit; limit = new_limit;
return true; return true;
} }
bool ScaleLimit(PointIndex pi, double factor) { bool ScaleLimit (PointIndex pi, double factor)
double &limit = (pi <= tool.np) ? limits[pi] : limits[map_from[pi]]; {
double& limit = (pi <= tool.np) ? limits[pi] : limits[map_from[pi]];
return SetLimit(pi, limit * factor); return SetLimit(pi, limit * factor);
} }
Vec<3> GetVector(PointIndex pi_to, double shift = 1., Vec<3> GetVector (PointIndex pi_to, double shift = 1., bool apply_limit = false)
bool apply_limit = false) { {
auto [gw, height] = tool.growth_vector_map[pi_to]; auto [gw, height] = tool.growth_vector_map[pi_to];
if (apply_limit) if (apply_limit)
shift *= GetLimit(pi_to); shift *= GetLimit(pi_to);
return shift * height * (*gw); return shift * height * (*gw);
} }
Point<3> GetPoint(PointIndex pi_to, double shift = 1., Point<3> GetPoint (PointIndex pi_to, double shift = 1., bool apply_limit = false)
bool apply_limit = false) { {
if (pi_to <= tool.np || tool.growth_vector_map.count(pi_to) == 0) if (pi_to <= tool.np || tool.growth_vector_map.count(pi_to) == 0)
return mesh[pi_to]; return mesh[pi_to];
return mesh[pi_to] + GetVector(pi_to, shift, apply_limit); return mesh[pi_to] + GetVector(pi_to, shift, apply_limit);
} }
Point<3> GetMappedPoint(PointIndex pi_from, double shift = 1.) { Point<3> GetMappedPoint (PointIndex pi_from, double shift = 1.)
{
auto pi_to = tool.mapto[pi_from].Last(); auto pi_to = tool.mapto[pi_from].Last();
return GetPoint(pi_to, shift); return GetPoint(pi_to, shift);
} }
Seg GetMappedSeg(PointIndex pi_from, double shift = 1.) { Seg GetMappedSeg (PointIndex pi_from, double shift = 1.)
{
return {mesh[pi_from], GetMappedPoint(pi_from, shift)}; return {mesh[pi_from], GetMappedPoint(pi_from, shift)};
} }
Seg GetSeg(PointIndex pi_to, double shift = 1., bool apply_limit = false) { Seg GetSeg (PointIndex pi_to, double shift = 1., bool apply_limit = false)
{
return {GetPoint(pi_to, 0), GetPoint(pi_to, shift, apply_limit)}; return {GetPoint(pi_to, 0), GetPoint(pi_to, shift, apply_limit)};
} }
Trig GetTrig(SurfaceElementIndex sei, double shift = 0.0, Trig GetTrig (SurfaceElementIndex sei, double shift = 0.0, bool apply_limit = false)
bool apply_limit = false) { {
auto sel = Get(sei); auto sel = Get(sei);
Trig trig; Trig trig;
for (auto i : Range(3)) for (auto i : Range(3))
@ -121,7 +132,8 @@ struct GrowthVectorLimiter {
return trig; return trig;
} }
Trig GetMappedTrig(SurfaceElementIndex sei, double shift = 0.0) { Trig GetMappedTrig (SurfaceElementIndex sei, double shift = 0.0)
{
auto sel = Get(sei); auto sel = Get(sei);
Trig trig; Trig trig;
for (auto i : Range(3)) for (auto i : Range(3))
@ -129,8 +141,8 @@ struct GrowthVectorLimiter {
return trig; return trig;
} }
Trig GetSideTrig(SurfaceElementIndex sei, int index, double shift = 0.0, Trig GetSideTrig (SurfaceElementIndex sei, int index, double shift = 0.0, bool grow_first_vertex = true)
bool grow_first_vertex = true) { {
auto trig = GetMappedTrig(sei, 0.0); auto trig = GetMappedTrig(sei, 0.0);
auto sel = Get(sei); auto sel = Get(sei);
auto index1 = (index + 1) % 3; auto index1 = (index + 1) % 3;
@ -141,116 +153,126 @@ struct GrowthVectorLimiter {
} }
static constexpr double INTERSECTION_SAFETY = .9; static constexpr double INTERSECTION_SAFETY = .9;
bool LimitGrowthVector(PointIndex pi_to, SurfaceElementIndex sei, bool LimitGrowthVector (PointIndex pi_to, SurfaceElementIndex sei, double trig_shift, double seg_shift, bool check_prism_sides = false)
double trig_shift, double seg_shift, {
bool check_prism_sides = false) {
auto pi_from = map_from[pi_to]; auto pi_from = map_from[pi_to];
if (!pi_from.IsValid()) if (!pi_from.IsValid())
return false; return false;
auto seg = GetSeg(pi_to, seg_shift, true); auto seg = GetSeg(pi_to, seg_shift, true);
for (auto pi : Get(sei).PNums()) { for (auto pi : Get(sei).PNums())
if (pi == pi_from) {
return false; if (pi == pi_from)
if (map_from[pi] == pi_from) return false;
return false; if (map_from[pi] == pi_from)
} return false;
if (check_prism_sides || trig_shift > .0) {
auto [trig_min_limit, trig_max_limit] = GetMinMaxLimit(sei);
if (GetLimit(pi_to) < trig_min_limit)
return false;
auto getTrigs = [&](double scaling = 1.0) -> ArrayMem<Trig, 3> {
ArrayMem<Trig, 3> trigs;
if (check_prism_sides)
for (auto i : Range(3))
trigs.Append(GetSideTrig(sei, i, scaling * trig_shift, true));
else
trigs.Append(GetTrig(sei, scaling * trig_shift, true));
return trigs;
};
double scaling = 1.0;
while (true) {
bool have_intersection = false;
auto seg = GetSeg(pi_to, scaling * seg_shift, true);
for (auto trig : getTrigs(scaling))
have_intersection |= isIntersectingTrig(seg, trig);
if (!have_intersection)
break;
scaling *= 0.9;
} }
if (scaling == 1.0)
return false;
double new_limit = scaling * max(GetLimit(pi_to), trig_max_limit); if (check_prism_sides || trig_shift > .0)
SetLimit(pi_to, new_limit); {
for (auto pi : Get(sei).PNums()) auto [trig_min_limit, trig_max_limit] = GetMinMaxLimit(sei);
SetLimit(pi, new_limit); if (GetLimit(pi_to) < trig_min_limit)
return true; return false;
} else {
auto trig = GetTrig(sei, 0.0); auto getTrigs = [&] (double scaling = 1.0) -> ArrayMem<Trig, 3> {
auto intersection = isIntersectingTrig(seg, trig); ArrayMem<Trig, 3> trigs;
// checking with original surface elements -> allow only half the distance if (check_prism_sides)
auto new_seg_limit = 0.40 * intersection.lam0 * seg_shift; for (auto i : Range(3))
if (intersection && new_seg_limit < GetLimit(pi_from)) { trigs.Append(GetSideTrig(sei, i, scaling * trig_shift, true));
auto p0 = seg[0]; else
auto p1 = seg[1]; trigs.Append(GetTrig(sei, scaling * trig_shift, true));
auto d = Dist(p0, p1); return trigs;
auto [gw, height] = tool.growth_vector_map[pi_to]; };
return SetLimit(pi_from, new_seg_limit);
double scaling = 1.0;
while (true)
{
bool have_intersection = false;
auto seg = GetSeg(pi_to, scaling * seg_shift, true);
for (auto trig : getTrigs(scaling))
have_intersection |= isIntersectingTrig(seg, trig);
if (!have_intersection)
break;
scaling *= 0.9;
}
if (scaling == 1.0)
return false;
double new_limit = scaling * max(GetLimit(pi_to), trig_max_limit);
SetLimit(pi_to, new_limit);
for (auto pi : Get(sei).PNums())
SetLimit(pi, new_limit);
return true;
}
else
{
auto trig = GetTrig(sei, 0.0);
auto intersection = isIntersectingTrig(seg, trig);
// checking with original surface elements -> allow only half the distance
auto new_seg_limit = 0.40 * intersection.lam0 * seg_shift;
if (intersection && new_seg_limit < GetLimit(pi_from))
{
auto p0 = seg[0];
auto p1 = seg[1];
auto d = Dist(p0, p1);
auto [gw, height] = tool.growth_vector_map[pi_to];
return SetLimit(pi_from, new_seg_limit);
}
return false;
} }
return false;
}
} }
void EqualizeLimits(double factor = .5) { void EqualizeLimits (double factor = .5)
{
static Timer t("GrowthVectorLimiter::EqualizeLimits"); static Timer t("GrowthVectorLimiter::EqualizeLimits");
PrintMessage(5, "GrowthVectorLimiter - equalize limits"); PrintMessage(5, "GrowthVectorLimiter - equalize limits");
RegionTimer reg(t); RegionTimer reg(t);
if (factor == 0.0) if (factor == 0.0)
return; return;
for (PointIndex pi : IntRange(tool.np, mesh.GetNP())) { for (PointIndex pi : IntRange(tool.np, mesh.GetNP()))
auto pi_from = map_from[pi]; {
std::set<PointIndex> pis; auto pi_from = map_from[pi];
for (auto sei : p2sel[pi]) std::set<PointIndex> pis;
for (auto pi_ : tool.new_sels[sei].PNums()) for (auto sei : p2sel[pi])
pis.insert(pi_); for (auto pi_ : tool.new_sels[sei].PNums())
ArrayMem<double, 20> limits; pis.insert(pi_);
for (auto pi1 : pis) { ArrayMem<double, 20> limits;
auto limit = GetLimit(pi1); for (auto pi1 : pis)
if (limit > 0.0) {
limits.Append(GetLimit(pi1)); auto limit = GetLimit(pi1);
if (limit > 0.0)
limits.Append(GetLimit(pi1));
}
if (limits.Size() == 0)
continue;
QuickSort(limits);
double mean_limit = limits[limits.Size() / 2];
// if mean limit is the maximum limit, take the average of second-highest
// and highest value
if (mean_limit > limits[0] && mean_limit == limits.Last())
{
auto i = limits.Size() - 1;
while (limits[i] == limits.Last())
i--;
mean_limit = 0.5 * (limits[i] + limits.Last());
}
if (limits.Size() % 2 == 0)
mean_limit = 0.5 * (mean_limit + limits[(limits.Size() - 1) / 2]);
SetLimit(pi, factor * mean_limit + (1.0 - factor) * GetLimit(pi));
} }
if (limits.Size() == 0)
continue;
QuickSort(limits);
double mean_limit = limits[limits.Size() / 2];
// if mean limit is the maximum limit, take the average of second-highest
// and highest value
if (mean_limit > limits[0] && mean_limit == limits.Last()) {
auto i = limits.Size() - 1;
while (limits[i] == limits.Last())
i--;
mean_limit = 0.5 * (limits[i] + limits.Last());
}
if (limits.Size() % 2 == 0)
mean_limit = 0.5 * (mean_limit + limits[(limits.Size() - 1) / 2]);
SetLimit(pi, factor * mean_limit + (1.0 - factor) * GetLimit(pi));
}
} }
void LimitSelfIntersection(double safety = 1.4) { void LimitSelfIntersection (double safety = 1.4)
{
static Timer t("GrowthVectorLimiter::LimitSelfIntersection"); static Timer t("GrowthVectorLimiter::LimitSelfIntersection");
PrintMessage(5, "GrowthVectorLimiter - self intersection"); PrintMessage(5, "GrowthVectorLimiter - self intersection");
RegionTimer reg(t); RegionTimer reg(t);
// check for self-intersection within new elements (prisms/hexes) // check for self-intersection within new elements (prisms/hexes)
auto isIntersecting = [&](SurfaceElementIndex sei, double shift) { auto isIntersecting = [&] (SurfaceElementIndex sei, double shift) {
// checks if surface element is self intersecting when growing with factor // checks if surface element is self intersecting when growing with factor
// shift // shift
@ -260,52 +282,59 @@ struct GrowthVectorLimiter {
return false; return false;
const auto sel = Get(sei); const auto sel = Get(sei);
auto np = sel.GetNP(); auto np = sel.GetNP();
for (auto i : Range(np)) { for (auto i : Range(np))
if (sel[i] > tool.np) {
return false; if (sel[i] > tool.np)
if (tool.mapto[sel[i]].Size() == 0) return false;
return false; if (tool.mapto[sel[i]].Size() == 0)
} return false;
for (auto i : Range(np)) { }
auto seg = GetMappedSeg(sel[i], shift * limits[sel[i]]); for (auto i : Range(np))
for (auto fi : Range(np - 2)) { {
for (auto side : {true, false}) { auto seg = GetMappedSeg(sel[i], shift * limits[sel[i]]);
auto trig = GetSideTrig(sei, i + fi, 1.0, side); for (auto fi : Range(np - 2))
if (isIntersectingPlane(seg, trig)) {
return true; for (auto side : {true, false})
} {
auto trig = GetSideTrig(sei, i + fi, 1.0, side);
if (isIntersectingPlane(seg, trig))
return true;
}
}
} }
}
return false; return false;
}; };
for (SurfaceElementIndex sei : mesh.SurfaceElements().Range()) { for (SurfaceElementIndex sei : mesh.SurfaceElements().Range())
auto sel = mesh[sei]; {
const auto &fd = mesh.GetFaceDescriptor(sel.GetIndex()); auto sel = mesh[sei];
if (sel.GetNP() == 4) const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex());
continue; if (sel.GetNP() == 4)
continue;
auto np = sel.GetNP(); auto np = sel.GetNP();
double shift = 1.0; double shift = 1.0;
const double step_factor = 0.9; const double step_factor = 0.9;
while (isIntersecting(sei, shift * safety)) { while (isIntersecting(sei, shift * safety))
shift *= step_factor; {
double max_limit = 0; shift *= step_factor;
for (auto i : Range(np)) double max_limit = 0;
max_limit = max(max_limit, limits[sel[i]]); for (auto i : Range(np))
for (auto i : Range(np)) max_limit = max(max_limit, limits[sel[i]]);
if (max_limit == limits[sel[i]]) for (auto i : Range(np))
ScaleLimit(sel[i], step_factor); if (max_limit == limits[sel[i]])
// if (max_limit < 0.01) break; ScaleLimit(sel[i], step_factor);
// if (max_limit < 0.01) break;
}
} }
}
} }
// checks if a segment is intersecting a plane, spanned by three points, lam // checks if a segment is intersecting a plane, spanned by three points, lam
// will be set s.t. p_intersect = seg[0] + lam * (seg[1]-seg[0]) // will be set s.t. p_intersect = seg[0] + lam * (seg[1]-seg[0])
Intersection_ isIntersectingPlane(const Seg & seg, Intersection_ isIntersectingPlane (const Seg& seg,
const Trig & trig) { const Trig& trig)
{
auto t1 = trig[1] - trig[0]; auto t1 = trig[1] - trig[0];
auto t2 = trig[2] - trig[0]; auto t2 = trig[2] - trig[0];
auto n = Cross(t1, t2); auto n = Cross(t1, t2);
@ -315,14 +344,13 @@ struct GrowthVectorLimiter {
Intersection_ intersection; Intersection_ intersection;
intersection.lam0 = -v0n / (v1n - v0n); intersection.lam0 = -v0n / (v1n - v0n);
intersection.p = seg[0] + intersection.lam0 * (seg[1] - seg[0]); intersection.p = seg[0] + intersection.lam0 * (seg[1] - seg[0]);
intersection.is_intersecting = (v0n * v1n < 0) && intersection.is_intersecting = (v0n * v1n < 0) && (intersection.lam0 > -1e-8) && (intersection.lam0 < 1 + 1e-8);
(intersection.lam0 > -1e-8) &&
(intersection.lam0 < 1 + 1e-8);
return intersection; return intersection;
} }
Intersection_ isIntersectingTrig(const Seg &seg, const Trig &trig) { Intersection_ isIntersectingTrig (const Seg& seg, const Trig& trig)
{
auto intersection = isIntersectingPlane(seg, trig); auto intersection = isIntersectingPlane(seg, trig);
if (!intersection) if (!intersection)
return intersection; return intersection;
@ -338,174 +366,186 @@ struct GrowthVectorLimiter {
intersection.lam1 = 0; intersection.lam1 = 0;
double eps = 0.1; double eps = 0.1;
if (bary.X() >= -eps && bary.Y() >= -eps && if (bary.X() >= -eps && bary.Y() >= -eps && bary.X() + bary.Y() <= 1 + eps)
bary.X() + bary.Y() <= 1 + eps) { {
intersection.bary[0] = bary.X(); intersection.bary[0] = bary.X();
intersection.bary[1] = bary.Y(); intersection.bary[1] = bary.Y();
intersection.bary[2] = 1.0 - bary.X() - bary.Y(); intersection.bary[2] = 1.0 - bary.X() - bary.Y();
} else }
else
intersection.is_intersecting = false; intersection.is_intersecting = false;
return intersection; return intersection;
} }
Intersection_ isIntersectingTrig(PointIndex pi_from, PointIndex pi_to, Intersection_ isIntersectingTrig (PointIndex pi_from, PointIndex pi_to, SurfaceElementIndex sei, double shift = 0.0)
SurfaceElementIndex sei, {
double shift = 0.0) {
return isIntersectingTrig(GetSeg(pi_from, pi_to), GetTrig(sei, shift)); return isIntersectingTrig(GetSeg(pi_from, pi_to), GetTrig(sei, shift));
} }
void BuildSearchTree(double trig_shift) { void BuildSearchTree (double trig_shift)
{
static Timer t("BuildSearchTree"); static Timer t("BuildSearchTree");
RegionTimer rt(t); RegionTimer rt(t);
Box<3> bbox(Box<3>::EMPTY_BOX); Box<3> bbox(Box<3>::EMPTY_BOX);
for (PointIndex pi : mesh.Points().Range()) { for (PointIndex pi : mesh.Points().Range())
bbox.Add(mesh[pi]); {
bbox.Add(GetPoint(pi, 1.1)); bbox.Add(mesh[pi]);
} bbox.Add(GetPoint(pi, 1.1));
}
tree = make_unique<BoxTree<3>>(bbox); tree = make_unique<BoxTree<3>>(bbox);
for (auto sei : SurfaceElementsRange()) { for (auto sei : SurfaceElementsRange())
const auto &sel = Get(sei); {
auto sel_index = sel.GetIndex(); const auto& sel = Get(sei);
auto sel_index = sel.GetIndex();
Box<3> box(Box<3>::EMPTY_BOX); Box<3> box(Box<3>::EMPTY_BOX);
for (auto pi : sel.PNums()) { for (auto pi : sel.PNums())
box.Add(GetPoint(pi, 0.)); {
box.Add(GetPoint(pi, trig_shift * GetLimit(pi))); box.Add(GetPoint(pi, 0.));
box.Add(GetPoint(pi, trig_shift * GetLimit(pi)));
}
tree->Insert(box, sei);
} }
tree->Insert(box, sei);
}
} }
template <typename TFunc> template <typename TFunc>
void FindTreeIntersections(double trig_shift, double seg_shift, TFunc f, void FindTreeIntersections (double trig_shift, double seg_shift, TFunc f, BitArray* relevant_points = nullptr)
BitArray *relevant_points = nullptr) { {
static Timer t("GrowthVectorLimiter::FindTreeIntersections"); static Timer t("GrowthVectorLimiter::FindTreeIntersections");
RegionTimer rt(t); RegionTimer rt(t);
BuildSearchTree(trig_shift); BuildSearchTree(trig_shift);
auto np_new = mesh.Points().Size(); auto np_new = mesh.Points().Size();
int counter = 0; int counter = 0;
for (auto i : IntRange(tool.np, np_new)) { for (auto i : IntRange(tool.np, np_new))
PointIndex pi_to = i + PointIndex::BASE; {
PointIndex pi_from = map_from[pi_to]; PointIndex pi_to = i + PointIndex::BASE;
if (!pi_from.IsValid()) PointIndex pi_from = map_from[pi_to];
throw Exception("Point not mapped"); if (!pi_from.IsValid())
throw Exception("Point not mapped");
if (relevant_points && !relevant_points->Test(pi_to) && if (relevant_points && !relevant_points->Test(pi_to) && !relevant_points->Test(pi_from))
!relevant_points->Test(pi_from)) continue;
continue;
Box<3> box(Box<3>::EMPTY_BOX); Box<3> box(Box<3>::EMPTY_BOX);
auto seg = GetSeg(pi_to, seg_shift); auto seg = GetSeg(pi_to, seg_shift);
box.Add(GetPoint(pi_to, 0)); box.Add(GetPoint(pi_to, 0));
box.Add(GetPoint(pi_to, GetLimit(pi_from))); box.Add(GetPoint(pi_to, GetLimit(pi_from)));
tree->GetFirstIntersecting(box.PMin(), box.PMax(), tree->GetFirstIntersecting(box.PMin(), box.PMax(), [&] (SurfaceElementIndex sei) {
[&](SurfaceElementIndex sei) { const auto& sel = Get(sei);
const auto &sel = Get(sei); if (sel.PNums().Contains(pi_from))
if (sel.PNums().Contains(pi_from)) return false;
return false; if (sel.PNums().Contains(pi_to))
if (sel.PNums().Contains(pi_to)) return false;
return false; counter++;
counter++; f(pi_to, sei);
f(pi_to, sei); return false;
return false; });
}); }
}
} }
void FixIntersectingSurfaceTrigs() { void FixIntersectingSurfaceTrigs ()
{
static Timer t("GrowthVectorLimiter::FixIntersectingSurfaceTrigs"); static Timer t("GrowthVectorLimiter::FixIntersectingSurfaceTrigs");
RegionTimer reg(t); RegionTimer reg(t);
// check if surface trigs are intersecting each other // check if surface trigs are intersecting each other
bool changed = true; bool changed = true;
while (changed) { while (changed)
changed = false; {
Point3d pmin, pmax; changed = false;
mesh.GetBox(pmin, pmax); Point3d pmin, pmax;
BoxTree<3, SurfaceElementIndex> setree(pmin, pmax); mesh.GetBox(pmin, pmax);
BoxTree<3, SurfaceElementIndex> setree(pmin, pmax);
for (auto sei : SurfaceElementsRange()) { for (auto sei : SurfaceElementsRange())
const Element2d &tri = Get(sei); {
const Element2d& tri = Get(sei);
Box<3> box(Box<3>::EMPTY_BOX); Box<3> box(Box<3>::EMPTY_BOX);
for (PointIndex pi : tri.PNums()) for (PointIndex pi : tri.PNums())
box.Add(GetPoint(pi, 1.0, true)); box.Add(GetPoint(pi, 1.0, true));
box.Increase(1e-3 * box.Diam()); box.Increase(1e-3 * box.Diam());
setree.Insert(box, sei); setree.Insert(box, sei);
}
for (auto sei : SurfaceElementsRange()) {
const Element2d &tri = Get(sei);
Box<3> box(Box<3>::EMPTY_BOX);
for (PointIndex pi : tri.PNums())
box.Add(GetPoint(pi, 1.0, true));
setree.GetFirstIntersecting(box.PMin(), box.PMax(), [&](size_t sej) {
const Element2d &tri2 = Get(sej);
if (mesh[tri[0]].GetLayer() != mesh[tri2[0]].GetLayer())
return false;
netgen::Point<3> tri1_points[3], tri2_points[3];
const netgen::Point<3> *trip1[3], *trip2[3];
for (int k = 0; k < 3; k++) {
trip1[k] = &tri1_points[k];
trip2[k] = &tri2_points[k];
} }
auto set_points = [&]() {
for (int k = 0; k < 3; k++) {
tri1_points[k] = GetPoint(tri[k], 1.0, true);
tri2_points[k] = GetPoint(tri2[k], 1.0, true);
}
};
set_points(); for (auto sei : SurfaceElementsRange())
{
const Element2d& tri = Get(sei);
int counter = 0; Box<3> box(Box<3>::EMPTY_BOX);
while (IntersectTriangleTriangle(&trip1[0], &trip2[0])) { for (PointIndex pi : tri.PNums())
changed = true; box.Add(GetPoint(pi, 1.0, true));
PointIndex pi_max_limit = PointIndex::INVALID;
for (PointIndex pi :
{tri[0], tri[1], tri[2], tri2[0], tri2[1], tri2[2]})
if (pi > tool.np && (!pi_max_limit.IsValid() ||
GetLimit(pi) > GetLimit(pi_max_limit)))
pi_max_limit = map_from[pi];
if (!pi_max_limit.IsValid()) setree.GetFirstIntersecting(box.PMin(), box.PMax(), [&] (size_t sej) {
break; const Element2d& tri2 = Get(sej);
ScaleLimit(pi_max_limit, 0.9); if (mesh[tri[0]].GetLayer() != mesh[tri2[0]].GetLayer())
set_points(); return false;
counter++;
if (counter > 20) { netgen::Point<3> tri1_points[3], tri2_points[3];
cerr << "Limit intersecting surface elements: too many " const netgen::Point<3>*trip1[3], *trip2[3];
"limitation steps, sels: " for (int k = 0; k < 3; k++)
<< Get(sei) << '\t' << Get(sej) << endl; {
for (auto si : {sei, sej}) { trip1[k] = &tri1_points[k];
auto sel = Get(si); trip2[k] = &tri2_points[k];
cerr << "Limits: "; }
for (auto pi : sel.PNums()) auto set_points = [&] () {
cerr << GetLimit(pi) << ",\t"; for (int k = 0; k < 3; k++)
cerr << endl; {
for (auto pi : sel.PNums()) tri1_points[k] = GetPoint(tri[k], 1.0, true);
cerr << GetPoint(pi, 1.0, true) << "\t"; tri2_points[k] = GetPoint(tri2[k], 1.0, true);
cerr << endl; }
} };
cerr << "pi_max_limit " << pi_max_limit << endl;
break; set_points();
}
int counter = 0;
while (IntersectTriangleTriangle(&trip1[0], &trip2[0]))
{
changed = true;
PointIndex pi_max_limit = PointIndex::INVALID;
for (PointIndex pi :
{tri[0], tri[1], tri[2], tri2[0], tri2[1], tri2[2]})
if (pi > tool.np && (!pi_max_limit.IsValid() || GetLimit(pi) > GetLimit(pi_max_limit)))
pi_max_limit = map_from[pi];
if (!pi_max_limit.IsValid())
break;
ScaleLimit(pi_max_limit, 0.9);
set_points();
counter++;
if (counter > 20)
{
cerr << "Limit intersecting surface elements: too many "
"limitation steps, sels: "
<< Get(sei) << '\t' << Get(sej) << endl;
for (auto si : {sei, sej})
{
auto sel = Get(si);
cerr << "Limits: ";
for (auto pi : sel.PNums())
cerr << GetLimit(pi) << ",\t";
cerr << endl;
for (auto pi : sel.PNums())
cerr << GetPoint(pi, 1.0, true) << "\t";
cerr << endl;
}
cerr << "pi_max_limit " << pi_max_limit << endl;
break;
}
}
return false;
});
} }
return false;
});
} }
}
} }
void LimitOriginalSurface(double safety) { void LimitOriginalSurface (double safety)
{
static Timer t("GrowthVectorLimiter::LimitOriginalSurface"); static Timer t("GrowthVectorLimiter::LimitOriginalSurface");
RegionTimer reg(t); RegionTimer reg(t);
PrintMessage(5, "GrowthVectorLimiter - original surface"); PrintMessage(5, "GrowthVectorLimiter - original surface");
@ -513,14 +553,15 @@ struct GrowthVectorLimiter {
double trig_shift = 0; double trig_shift = 0;
double seg_shift = safety; double seg_shift = safety;
FindTreeIntersections( FindTreeIntersections(
trig_shift, seg_shift, [&](PointIndex pi_to, SurfaceElementIndex sei) { trig_shift, seg_shift, [&] (PointIndex pi_to, SurfaceElementIndex sei) {
if (sei >= tool.nse) if (sei >= tool.nse)
return; // ignore new surface elements in first pass return; // ignore new surface elements in first pass
LimitGrowthVector(pi_to, sei, trig_shift, seg_shift); LimitGrowthVector(pi_to, sei, trig_shift, seg_shift);
}); });
} }
void LimitBoundaryLayer(double safety = 1.1) { void LimitBoundaryLayer (double safety = 1.1)
{
static Timer t("GrowthVectorLimiter::LimitBoundaryLayer"); static Timer t("GrowthVectorLimiter::LimitBoundaryLayer");
PrintMessage(5, "GrowthVectorLimiter - boundary layer"); PrintMessage(5, "GrowthVectorLimiter - boundary layer");
// now limit again with shifted surface elements // now limit again with shifted surface elements
@ -533,45 +574,49 @@ struct GrowthVectorLimiter {
relevant_points_next.SetSize(mesh.Points().Size() + 1); relevant_points_next.SetSize(mesh.Points().Size() + 1);
relevant_points.Set(); relevant_points.Set();
while (limit_counter) { while (limit_counter)
RegionTimer reg(t); {
size_t find_counter = 0; RegionTimer reg(t);
limit_counter = 0; size_t find_counter = 0;
relevant_points_next.Clear(); limit_counter = 0;
FindTreeIntersections( relevant_points_next.Clear();
trig_shift, seg_shift, FindTreeIntersections(
[&](PointIndex pi_to, SurfaceElementIndex sei) { trig_shift, seg_shift, [&] (PointIndex pi_to, SurfaceElementIndex sei) {
find_counter++; find_counter++;
auto sel = Get(sei); auto sel = Get(sei);
if (LimitGrowthVector(pi_to, sei, trig_shift, seg_shift)) { if (LimitGrowthVector(pi_to, sei, trig_shift, seg_shift))
limit_counter++; {
relevant_points_next.SetBit(pi_to); limit_counter++;
relevant_points_next.SetBit(map_from[pi_to]); relevant_points_next.SetBit(pi_to);
for (auto pi : sel.PNums()) { relevant_points_next.SetBit(map_from[pi_to]);
relevant_points_next.SetBit(pi); for (auto pi : sel.PNums())
if (pi >= tool.np) {
relevant_points_next.SetBit(map_from[pi]); relevant_points_next.SetBit(pi);
else if (pi >= tool.np)
relevant_points_next.SetBit(map_from[pi]); relevant_points_next.SetBit(map_from[pi]);
else
relevant_points_next.SetBit(map_from[pi]);
}
} }
}
for (auto pi : sel.PNums()) { for (auto pi : sel.PNums())
if (pi >= tool.np) {
return; if (pi >= tool.np)
if (tool.mapto[pi].Size() == 0) return;
return; if (tool.mapto[pi].Size() == 0)
} return;
}
if (LimitGrowthVector(pi_to, sei, trig_shift, seg_shift, true)) if (LimitGrowthVector(pi_to, sei, trig_shift, seg_shift, true))
limit_counter++; limit_counter++;
}, },
&relevant_points); &relevant_points);
relevant_points = relevant_points_next; relevant_points = relevant_points_next;
} }
} }
void Perform() { void Perform ()
{
limits.SetSize(mesh.Points().Size()); limits.SetSize(mesh.Points().Size());
limits = 1.0; limits = 1.0;
@ -580,31 +625,34 @@ struct GrowthVectorLimiter {
// No smoothing in the last pass, to avoid generating new intersections // No smoothing in the last pass, to avoid generating new intersections
std::array smoothing_factors = {0.8, 0.7, 0.5, 0.0}; std::array smoothing_factors = {0.8, 0.7, 0.5, 0.0};
for (auto i_pass : Range(safeties.size())) { for (auto i_pass : Range(safeties.size()))
PrintMessage(4, "GrowthVectorLimiter pass ", i_pass); {
double safety = safeties[i_pass]; PrintMessage(4, "GrowthVectorLimiter pass ", i_pass);
// intersect segment with original surface elements double safety = safeties[i_pass];
LimitOriginalSurface(2.1); // intersect segment with original surface elements
// intersect prisms with themself LimitOriginalSurface(2.1);
LimitSelfIntersection(1.3 * safety); // intersect prisms with themself
// intesect segment with prism LimitSelfIntersection(1.3 * safety);
LimitBoundaryLayer(safety); // intesect segment with prism
LimitBoundaryLayer(safety);
for (auto i : Range(3)) for (auto i : Range(3))
EqualizeLimits(smoothing_factors[i_pass]); EqualizeLimits(smoothing_factors[i_pass]);
if (i_pass == safeties.size() - 1) if (i_pass == safeties.size() - 1)
FixIntersectingSurfaceTrigs(); FixIntersectingSurfaceTrigs();
} }
for (auto i : Range(growthvectors)) for (auto i : Range(growthvectors))
growthvectors[i] *= limits[i]; growthvectors[i] *= limits[i];
for (auto &[special_pi, special_point] : tool.special_boundary_points) { for (auto& [special_pi, special_point] : tool.special_boundary_points)
for (auto &group : special_point.growth_groups) { {
group.growth_vector *= limits[special_pi]; for (auto& group : special_point.growth_groups)
{
group.growth_vector *= limits[special_pi];
}
} }
}
} }
}; };