Boundary Layers - Automatic thickness limiation and fixes

This commit is contained in:
Hochsteger, Matthias 2024-12-23 12:30:28 +01:00 committed by Lackner, Christopher
parent 975414c2fe
commit 1aa34da6af
19 changed files with 2804 additions and 1699 deletions

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

View File

@ -826,6 +826,9 @@ namespace netgen
{
multithread.task = "Volume meshing";
for (int i = 0; i < geom.GetNTopLevelObjects(); i++)
mesh->SetMaterial (i+1, geom.GetTopLevelObject(i)->GetMaterial().c_str());
MESHING3_RESULT res =
MeshVolume (mparam, *mesh);
@ -838,10 +841,6 @@ namespace netgen
MeshQuality3d (*mesh);
for (int i = 0; i < geom.GetNTopLevelObjects(); i++)
mesh->SetMaterial (i+1, geom.GetTopLevelObject(i)->GetMaterial().c_str());
#ifdef STAT_STREAM
(*statout) << GetTime() << " & ";
#endif

View File

@ -12,7 +12,7 @@ target_sources(nglib PRIVATE
parallelmesh.cpp paralleltop.cpp basegeom.cpp
python_mesh.cpp surfacegeom.cpp
debugging.cpp fieldlines.cpp visual_interface.cpp
boundarylayer2d.cpp
boundarylayer2d.cpp boundarylayer_interpolate.cpp
)
target_link_libraries( nglib PRIVATE $<BUILD_INTERFACE:netgen_metis> $<BUILD_INTERFACE:netgen_python> )

View File

@ -271,7 +271,7 @@ namespace netgen
virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, EdgePointGeomInfo* gi = nullptr) const
{
if(gi && gi->edgenr < edges.Size())
if(gi && gi->edgenr < edges.Size() && gi->edgenr >= 0)
edges[gi->edgenr]->ProjectPoint(p, gi);
}

File diff suppressed because it is too large Load Diff

View File

@ -1,85 +1,110 @@
#ifndef NETGEN_BOUNDARYLAYER_HPP
#define NETGEN_BOUNDARYLAYER_HPP
#include <core/array.hpp>
#include <mystdlib.h>
#include <meshing.hpp>
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
DLL_HEADER void GenerateBoundaryLayer (Mesh & mesh,
const BoundaryLayerParameters & blp);
struct SpecialBoundaryPoint
{
struct GrowthGroup
{
Array<int> faces;
Vec<3> growth_vector;
Array<PointIndex> new_points;
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>{});
GrowthGroup(FlatArray<int> faces_, FlatArray<Vec<3>> normals);
GrowthGroup(const GrowthGroup&) = default;
GrowthGroup() = default;
};
Array<GrowthGroup> growth_groups;
Vec<3> separating_direction;
SpecialBoundaryPoint(const std::map<int, Vec<3>>& normals);
SpecialBoundaryPoint() = default;
};
DLL_HEADER void GenerateBoundaryLayer (Mesh& mesh,
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>{});
class BoundaryLayerTool
{
public:
BoundaryLayerTool(Mesh & mesh_, const BoundaryLayerParameters & params_);
void ProcessParameters();
void Perform();
public:
BoundaryLayerTool(Mesh& mesh_, const BoundaryLayerParameters& params_);
void ProcessParameters ();
void Perform ();
protected:
Mesh & mesh;
MeshTopology & topo;
BoundaryLayerParameters params;
Array<Vec<3>, PointIndex> growthvectors;
Table<SurfaceElementIndex, PointIndex> p2sel;
Mesh& mesh;
MeshTopology& topo;
BoundaryLayerParameters params;
Array<Vec<3>, PointIndex> growthvectors;
std::map<PointIndex, Vec<3>> non_bl_growth_vectors;
Table<SurfaceElementIndex, PointIndex> p2sel;
BitArray domains, is_edge_moved, is_boundary_projected, is_boundary_moved;
Array<SegmentIndex> moved_segs;
int max_edge_nr, nfd_old, ndom_old;
Array<int> new_mat_nrs;
BitArray moved_surfaces;
int np, nseg, nse, ne;
double height;
BitArray domains, is_edge_moved, is_boundary_projected, is_boundary_moved;
Array<SegmentIndex> moved_segs;
int max_edge_nr, nfd_old, ndom_old;
Array<int> new_mat_nrs;
BitArray moved_surfaces;
int np, nseg, nse, ne;
double total_height;
Array<POINTTYPE, PointIndex> point_types;
// These parameters are derived from given BoundaryLayerParameters and the Mesh
Array<double> par_heights;
Array<int> par_surfid;
map<string, string> par_new_mat;
Array<size_t> par_project_boundaries;
// These parameters are derived from given BoundaryLayerParameters and the Mesh
Array<double> par_heights;
Array<int> par_surfid;
bool insert_only_volume_elements;
map<string, string> par_new_mat;
Array<size_t> par_project_boundaries;
bool have_single_segments;
Array<Segment> segments, new_segments;
bool have_single_segments;
Array<Segment> old_segments, segments, new_segments, new_segments_on_moved_bnd;
Array<Element2d, SurfaceElementIndex> new_sels, new_sels_on_moved_bnd;
Array<Array<PointIndex>, PointIndex> mapto;
Array<PointIndex, PointIndex> mapfrom;
Array<double> surfacefacs;
Array<int> si_map;
Array<double, PointIndex> limits;
Array<double> surfacefacs;
Array<int> si_map;
// major steps called in Perform()
void CreateNewFaceDescriptors();
void CreateFaceDescriptorsSides();
void CalculateGrowthVectors();
Array<Array<pair<SegmentIndex, int>>, SegmentIndex> BuildSegMap();
std::map<PointIndex, SpecialBoundaryPoint> special_boundary_points;
std::map<PointIndex, std::tuple<Vec<3>*, double>> growth_vector_map;
BitArray ProjectGrowthVectorsOnSurface();
void InterpolateSurfaceGrowthVectors();
void InterpolateGrowthVectors();
void LimitGrowthVectorLengths();
// major steps called in Perform()
void CreateNewFaceDescriptors ();
void CreateFaceDescriptorsSides ();
void CalculateGrowthVectors ();
Array<Array<pair<SegmentIndex, int>>, SegmentIndex> BuildSegMap ();
void InsertNewElements(FlatArray<Array<pair<SegmentIndex, int>>, SegmentIndex> segmap, const BitArray & in_surface_direction);
void SetDomInOut();
void SetDomInOutSides();
void AddSegments();
void FixVolumeElements();
BitArray ProjectGrowthVectorsOnSurface ();
void InterpolateSurfaceGrowthVectors ();
void InterpolateGrowthVectors ();
void LimitGrowthVectorLengths ();
void FixSurfaceElements ();
// utility functions
array<Point<3>, 2> GetMappedSeg( PointIndex pi );
ArrayMem<Point<3>, 4> GetFace( SurfaceElementIndex sei );
ArrayMem<Point<3>, 4> GetMappedFace( SurfaceElementIndex sei );
ArrayMem<Point<3>, 4> GetMappedFace( SurfaceElementIndex sei, int face );
void InsertNewElements (FlatArray<Array<pair<SegmentIndex, int>>, SegmentIndex> segmap, const BitArray& in_surface_direction);
void SetDomInOut ();
void SetDomInOutSides ();
void AddSegments ();
void AddSurfaceElements ();
Vec<3> getNormal(const Element2d & el)
{
auto v0 = mesh[el[0]];
return Cross(mesh[el[1]]-v0, mesh[el[2]]-v0).Normalize();
}
Vec<3> getNormal (const Element2d& el)
{
auto v0 = mesh[el[0]];
return Cross(mesh[el[1]] - v0, mesh[el[2]] - v0).Normalize();
}
Vec<3> getEdgeTangent(PointIndex pi, int edgenr);
Vec<3> getEdgeTangent (PointIndex pi, int edgenr, FlatArray<Segment*> segs);
};
} // namespace netgen

View File

@ -1,5 +1,5 @@
#include <mystdlib.h>
#include "meshing.hpp"
#include "boundarylayer.hpp"
#include "meshing2.hpp"
#include "../geom2d/csg2d.hpp"

View File

@ -0,0 +1,472 @@
#include "boundarylayer.hpp"
namespace netgen
{
namespace detail
{
struct Neighbor
{
PointIndex pi;
SurfaceElementIndex sei;
double weight;
};
} // namespace detail
Array<ArrayMem<detail::Neighbor, 20>>
BuildNeighbors (FlatArray<PointIndex> points, const Mesh& mesh)
{
auto p2sel = mesh.CreatePoint2SurfaceElementTable();
Array<ArrayMem<detail::Neighbor, 20>> neighbors(points.Size());
ArrayMem<double, 20> angles;
ArrayMem<double, 20> inv_dists;
for (auto i : points.Range())
{
auto& p_neighbors = neighbors[i];
auto pi = points[i];
angles.SetSize(0);
inv_dists.SetSize(0);
for (auto sei : p2sel[pi])
{
const auto& sel = mesh[sei];
for (auto pi1 : sel.PNums())
{
if (pi1 == pi)
continue;
auto pi2 = pi1;
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));
}
}
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;
for (auto i : Range(inv_dists))
{
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;
}
return neighbors;
}
void BoundaryLayerTool ::InterpolateGrowthVectors()
{
point_types.SetSize(mesh.GetNP());
for (auto p : mesh.Points().Range())
point_types[p] = mesh[p].Type();
int new_max_edge_nr = max_edge_nr;
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;
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;
};
// interpolate tangential component of growth vector along edge
if (max_edge_nr >= new_max_edge_nr)
return;
auto edgenr2seg = ngcore::CreateSortedTable<Segment*, int>(
Range(segments.Size() + new_segments.Size()),
[&] (auto& table, size_t segi) {
auto& seg = segi < segments.Size()
? segments[segi]
: new_segments[segi - segments.Size()];
table.Add(seg.edgenr, &seg);
},
new_max_edge_nr + 1);
auto point2seg = ngcore::CreateSortedTable<Segment*, PointIndex>(
Range(segments.Size() + new_segments.Size()),
[&] (auto& table, size_t segi) {
auto& seg = segi < segments.Size()
? segments[segi]
: new_segments[segi - segments.Size()];
table.Add(seg[0], &seg);
table.Add(seg[1], &seg);
},
mesh.GetNP());
for (auto edgenr : Range(1, new_max_edge_nr + 1))
{
// "inner" edges between two flat faces are not treated as edges for interpolation
bool no_angles = true;
ArrayMem<SurfaceElementIndex, 4> faces;
for (auto* p_seg : edgenr2seg[edgenr])
{
auto& seg = *p_seg;
faces.SetSize(0);
if (seg[0] <= p2sel.Size())
{
for (auto sei : p2sel[seg[0]])
if (moved_surfaces.Test(mesh[sei].GetIndex()) && p2sel[seg[1]].Contains(sei))
faces.Append(sei);
}
if (faces.Size() == 2)
{
auto n0 = getNormal(mesh[faces[0]]);
auto n1 = getNormal(mesh[faces[1]]);
if (n0 * n1 < 0.999)
no_angles = false;
}
else
{
no_angles = false;
}
}
if (no_angles)
{
for (auto* p_seg : edgenr2seg[edgenr])
for (auto pi : p_seg->PNums())
if (pi <= np && point_types[pi] == EDGEPOINT)
point_types[pi] = SURFACEPOINT;
continue;
}
}
for (auto edgenr : Range(max_edge_nr + 1, new_max_edge_nr + 1))
{
double edge_len = 0.;
bool any_grows = false;
auto is_end_point = [&] (PointIndex pi) {
auto segs = point2seg[pi];
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 false;
};
Array<PointIndex> points;
for (auto* p_seg : edgenr2seg[edgenr])
{
auto& seg = *p_seg;
if (getGW(seg[0]).Length2() != 0 || getGW(seg[1]).Length2() != 0)
any_grows = true;
if (points.Size() == 0)
for (auto i : Range(2))
if (is_end_point(seg[i]))
{
points.Append(seg[i]);
points.Append(seg[1 - i]);
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())
{
if (debugparam.debugoutput)
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))
continue;
edge_len += (mesh[points.Last()] - mesh[pnew]).Length();
points.Append(pnew);
points_set.insert(pnew);
point_found = true;
break;
}
}
if (!point_found)
{
if (debugparam.debugoutput)
{
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();
auto gt2 = getGW(points.Last()) * t2 * t2;
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, point2seg[pi]);
auto lam = len / edge_len;
auto interpol = (1 - lam) * (gt1 * t) * t + lam * (gt2 * t) * t;
addGW(pi, interpol);
}
}
}
void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors()
{
static Timer tall("InterpolateSurfaceGrowthVectors");
RegionTimer rtall(tall);
static Timer tsmooth("InterpolateSurfaceGrowthVectors-Smoothing");
auto np_old = this->np;
auto np = mesh.GetNP();
auto hasMoved = [&] (PointIndex pi) {
return (pi - PointIndex::BASE >= np_old) || mapto[pi].Size() > 0 || special_boundary_points.count(pi);
};
std::set<PointIndex> points_set;
for (const auto& sel : mesh.SurfaceElements())
{
for (auto pi : sel.PNums())
if (point_types[pi] == SURFACEPOINT && hasMoved(pi))
points_set.insert(pi);
}
Array<PointIndex> points;
for (auto pi : points_set)
points.Append(pi);
QuickSort(points);
// smooth tangential part of growth vectors from edges to surface elements
Array<Vec<3>, PointIndex> corrections(mesh.GetNP());
corrections = 0.0;
RegionTimer rtsmooth(tsmooth);
auto neighbors = BuildNeighbors(points, mesh);
Array<Vec<3>, SurfaceElementIndex> surf_normals(mesh.GetNSE());
for (auto sei : mesh.SurfaceElements().Range())
surf_normals[sei] = getNormal(mesh[sei]);
BitArray interpolate_tangent(mesh.GetNP() + 1);
interpolate_tangent = false;
for (auto pi : points)
{
for (auto sei : p2sel[pi])
if (is_boundary_moved[mesh[sei].GetIndex()])
interpolate_tangent.SetBit(pi);
}
constexpr int N_STEPS = 64;
for ([[maybe_unused]] auto i : Range(N_STEPS))
{
for (auto i : points.Range())
{
auto pi = points[i];
auto& p_neighbors = neighbors[i];
ArrayMem<Vec<3>, 20> g_vectors;
double max_len = 0.0;
double sum_len = 0.0;
// average only tangent component on new bl points, average whole growth
// vector otherwise
bool do_average_tangent = true;
for (const auto& s : p_neighbors)
{
auto gw_other = growthvectors[s.pi] + corrections[s.pi];
if (do_average_tangent)
{
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;
correction += weight * v;
}
if (!do_average_tangent)
correction -= growthvectors[pi];
}
}
for (auto pi : points)
growthvectors[pi] += corrections[pi];
}
void BoundaryLayerTool ::FixSurfaceElements()
{
static Timer tall("FixSurfaceElements");
RegionTimer rtall(tall);
auto np_old = this->np;
auto np = mesh.GetNP();
non_bl_growth_vectors.clear();
auto getGW = [&] (PointIndex pi) -> Vec<3> {
// return growthvectors[pi];
if (growth_vector_map.count(pi) == 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];
return height * (*gw);
};
auto addGW = [&] (PointIndex pi, Vec<3> vec) {
if (growth_vector_map.count(pi) == 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];
*gw += 1.0 / height * vec;
};
std::set<PointIndex> points_set;
// only smooth over old surface elements
for (SurfaceElementIndex sei : Range(nse))
{
const auto& sel = mesh[sei];
if (sel.GetNP() == 3 && is_boundary_moved[sel.GetIndex()])
for (auto pi : sel.PNums())
if (point_types[pi] == SURFACEPOINT)
points_set.insert(pi);
}
Array<PointIndex> points;
for (auto pi : points_set)
points.Append(pi);
QuickSort(points);
Array<Vec<3>, PointIndex> corrections(mesh.GetNP());
corrections = 0.0;
auto neighbors = BuildNeighbors(points, mesh);
constexpr int N_STEPS = 32;
for ([[maybe_unused]] auto i : Range(N_STEPS))
{
for (auto i : points.Range())
{
auto pi = points[i];
auto& p_neighbors = neighbors[i];
ArrayMem<Vec<3>, 20> g_vectors;
double max_len = 0.0;
double sum_len = 0.0;
for (const auto& s : p_neighbors)
{
auto v = getGW(s.pi) + corrections[s.pi];
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;
correction += weight * v;
}
}
}
for (auto pi : points)
addGW(pi, corrections[pi]);
}
} // namespace netgen

View File

@ -0,0 +1,752 @@
#include "boundarylayer.hpp"
#include <core/array.hpp>
namespace netgen
{
struct Intersection_
{
bool is_intersecting = false;
double lam0 = -1, lam1 = -1;
Point<3> p;
double bary[3];
operator bool() const { return is_intersecting; }
};
struct GrowthVectorLimiter
{
typedef std::array<Point<3>, 2> Seg;
typedef std::array<Point<3>, 3> Trig;
BoundaryLayerTool& tool;
const BoundaryLayerParameters& params;
Mesh& mesh;
double height;
Array<double, PointIndex> limits;
FlatArray<Vec<3>, PointIndex> growthvectors;
BitArray changed_domains;
unique_ptr<BoxTree<3>> tree;
Array<PointIndex, PointIndex> map_from;
Table<SurfaceElementIndex, PointIndex> p2sel;
GrowthVectorLimiter(BoundaryLayerTool& tool_)
: tool(tool_), params(tool_.params), mesh(tool_.mesh), height(tool_.total_height), growthvectors(tool_.growthvectors), map_from(mesh.Points().Size())
{
changed_domains = tool.domains;
if (!params.outside)
changed_domains.Invert();
map_from = tool.mapfrom;
p2sel = ngcore::CreateSortedTable<SurfaceElementIndex, PointIndex>(
tool.new_sels.Range(),
[&] (auto& table, SurfaceElementIndex ei) {
for (PointIndex pi : tool.new_sels[ei].PNums())
table.Add(pi, ei);
},
mesh.GetNP());
}
auto SurfaceElementsRange () { return Range(tool.nse + tool.new_sels.Size()); }
void WriteErrorMesh (string name)
{
if (!debugparam.write_mesh_on_error)
return;
Mesh out_mesh;
out_mesh = mesh;
for (auto [pi, data] : tool.growth_vector_map)
{
auto [gw, height] = data;
out_mesh[pi] += limits[pi] * height * (*gw);
}
out_mesh.Save(name);
}
const auto& Get (SurfaceElementIndex sei)
{
if (sei < tool.nse)
return mesh[sei];
return tool.new_sels[sei - tool.nse];
}
std::pair<double, double> GetMinMaxLimit (SurfaceElementIndex sei)
{
const auto& sel = Get(sei);
double min_limit = GetLimit(sel[0]);
double max_limit = min_limit;
for (auto i : IntRange(1, sel.GetNP()))
{
auto limit = GetLimit(sel[i]);
min_limit = min(min_limit, limit);
max_limit = max(max_limit, limit);
}
return {min_limit, max_limit};
}
double GetLimit (PointIndex pi)
{
if (pi <= tool.np)
return limits[pi];
return limits[map_from[pi]];
}
bool SetLimit (PointIndex pi, double new_limit)
{
double& limit = (pi <= tool.np) ? limits[pi] : limits[map_from[pi]];
if (limit <= new_limit)
return false;
limit = new_limit;
return true;
}
bool ScaleLimit (PointIndex pi, double factor)
{
double& limit = (pi <= tool.np) ? limits[pi] : limits[map_from[pi]];
return SetLimit(pi, limit * factor);
}
Vec<3> GetVector (PointIndex pi_to, double shift = 1., bool apply_limit = false)
{
auto [gw, height] = tool.growth_vector_map[pi_to];
if (apply_limit)
shift *= GetLimit(pi_to);
return shift * height * (*gw);
}
Point<3> GetPoint (PointIndex pi_to, double shift = 1., bool apply_limit = false)
{
if (pi_to <= tool.np || tool.growth_vector_map.count(pi_to) == 0)
return mesh[pi_to];
return mesh[pi_to] + GetVector(pi_to, shift, apply_limit);
}
Point<3> GetMappedPoint (PointIndex pi_from, double shift = 1., bool apply_limit = false)
{
auto pi_to = tool.mapto[pi_from].Last();
return GetPoint(pi_to, shift, apply_limit);
}
Seg GetMappedSeg (PointIndex pi_from, double shift = 1.)
{
return {mesh[pi_from], GetMappedPoint(pi_from, shift)};
}
Seg GetSeg (PointIndex pi_to, double shift = 1., bool apply_limit = false)
{
return {GetPoint(pi_to, 0), GetPoint(pi_to, shift, apply_limit)};
}
Trig GetTrig (SurfaceElementIndex sei, double shift = 0.0, bool apply_limit = false)
{
auto sel = Get(sei);
Trig trig;
for (auto i : Range(3))
trig[i] = GetPoint(sel[i], shift, apply_limit);
return trig;
}
Trig GetMappedTrig (SurfaceElementIndex sei, double shift = 0.0)
{
auto sel = Get(sei);
Trig trig;
for (auto i : Range(3))
trig[i] = GetMappedPoint(sel[i], shift);
return trig;
}
Trig GetSideTrig (SurfaceElementIndex sei, int index, double shift = 0.0, bool grow_first_vertex = true)
{
auto trig = GetMappedTrig(sei, 0.0);
auto sel = Get(sei);
auto index1 = (index + 1) % 3;
if (!grow_first_vertex)
index1 = (index + 2) % 3;
trig[index] = GetMappedPoint(sel[index1], shift, true);
return trig;
}
array<Trig, 4> GetSideTrigs (SurfaceElementIndex sei, int i0, double shift = 0.0)
{
auto trig = GetMappedTrig(sei, 0.0);
array<Trig, 4> trigs{trig, trig, trig, trig};
auto sel = Get(sei);
auto i1 = (i0 + 1) % 3;
auto i2 = (i0 + 2) % 3;
auto p1 = GetMappedPoint(sel[i1], shift, true);
auto p2 = GetMappedPoint(sel[i2], shift, true);
// create four trigs to span the quad from i1,i2 and their shifted points
// i1, i2, shifted i1
trigs[0][i0] = p1;
// i1, i2, shifted i2
trigs[1][i0] = p2;
// i1, shifted i1, shifted i2
trigs[2][i0] = p1;
trigs[2][i2] = p2;
// i2, shifted i1, shifted i2
trigs[2][i0] = p2;
trigs[2][i1] = p1;
return trigs;
}
static constexpr double INTERSECTION_SAFETY = .9;
bool LimitGrowthVector (PointIndex pi_to, SurfaceElementIndex sei, double trig_shift, double seg_shift, bool check_prism_sides = false)
{
auto pi_from = map_from[pi_to];
if (!pi_from.IsValid())
return false;
for (auto pi : Get(sei).PNums())
{
if (pi == pi_from)
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, 12> trigs;
if (check_prism_sides)
for (auto i : Range(3))
for (auto trig : GetSideTrigs(sei, i, scaling * trig_shift))
trigs.Append(trig);
else
trigs.Append(GetTrig(sei, scaling * trig_shift, true));
return trigs;
};
if (!check_prism_sides)
{
// If the growth vectors of all points are pointing in the same direction,
// an intersection means, we also have an intersection with a prism side face
// this is an extra check and handled later
auto seg = GetSeg(pi_to, 1.0, false);
auto gw = seg[1] - seg[0];
bool have_same_growth_direction = true;
for (auto pi : Get(sei).PNums())
{
auto p_seg = GetSeg(pi, 1.0, false);
auto p_gw = p_seg[1] - p_seg[0];
have_same_growth_direction &= (gw * p_gw) > 0;
}
if (have_same_growth_direction)
return false;
}
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 seg = GetSeg(pi_to, seg_shift, false);
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))
return SetLimit(pi_from, new_seg_limit);
return false;
}
}
void EqualizeLimits (double factor = .5)
{
static Timer t("GrowthVectorLimiter::EqualizeLimits");
PrintMessage(5, "GrowthVectorLimiter - equalize limits");
RegionTimer reg(t);
if (factor == 0.0)
return;
for (PointIndex pi : IntRange(tool.np, mesh.GetNP()))
{
auto pi_from = map_from[pi];
std::set<PointIndex> pis;
for (auto sei : p2sel[pi])
for (auto pi_ : tool.new_sels[sei].PNums())
pis.insert(pi_);
ArrayMem<double, 20> limits;
for (auto pi1 : pis)
{
auto limit = GetLimit(pi1);
if (limit > 0.0)
limits.Append(GetLimit(pi1));
}
if (limits.Size() == 0)
continue;
double average = 0.0;
for (auto l : limits)
average += l;
average /= limits.Size();
SetLimit(pi, factor * average + (1.0 - factor) * GetLimit(pi));
}
}
void LimitSelfIntersection (double safety = 1.4)
{
static Timer t("GrowthVectorLimiter::LimitSelfIntersection");
PrintMessage(5, "GrowthVectorLimiter - self intersection");
RegionTimer reg(t);
// check for self-intersection within new elements (prisms/hexes)
auto isIntersecting = [&] (SurfaceElementIndex sei, double shift) {
// checks if surface element is self intersecting when growing with factor
// shift
// ignore new surface elements, side trigs are only built
// from original surface elements
if (sei >= tool.nse)
return false;
const auto sel = Get(sei);
auto np = sel.GetNP();
for (auto i : Range(np))
{
if (sel[i] > tool.np)
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 fi : Range(np - 2))
{
for (auto side : {true, false})
{
auto trig = GetSideTrig(sei, i + fi, 1.0, side);
if (isIntersectingPlane(seg, trig))
return true;
}
}
}
return false;
};
for (SurfaceElementIndex sei : mesh.SurfaceElements().Range())
{
auto sel = mesh[sei];
if (!tool.moved_surfaces[sel.GetIndex()])
continue;
if (sel.GetNP() == 4)
continue;
const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex());
auto np = sel.GetNP();
double shift = 1.0;
const double step_factor = 0.9;
while (isIntersecting(sei, shift * safety))
{
shift *= step_factor;
double max_limit = 0;
for (auto i : Range(np))
max_limit = max(max_limit, limits[sel[i]]);
for (auto i : Range(np))
if (max_limit == limits[sel[i]])
ScaleLimit(sel[i], step_factor);
// if (max_limit < 0.01) break;
}
}
}
// 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])
Intersection_ isIntersectingPlane (const Seg& seg,
const Trig& trig)
{
auto t1 = trig[1] - trig[0];
auto t2 = trig[2] - trig[0];
auto n = Cross(t1, t2);
auto v0n = (seg[0] - trig[0]) * n;
auto v1n = (seg[1] - trig[0]) * n;
Intersection_ intersection;
intersection.lam0 = -v0n / (v1n - v0n);
intersection.p = seg[0] + intersection.lam0 * (seg[1] - seg[0]);
intersection.is_intersecting = (v0n * v1n < 0) && (intersection.lam0 > -1e-8) && (intersection.lam0 < 1 + 1e-8);
return intersection;
}
Intersection_ isIntersectingTrig (const Seg& seg, const Trig& trig)
{
auto intersection = isIntersectingPlane(seg, trig);
if (!intersection)
return intersection;
auto p = seg[0] + intersection.lam0 * (seg[1] - seg[0]) - trig[0];
Vec3d col1 = trig[1] - trig[0];
Vec3d col2 = trig[2] - trig[0];
Vec3d col3 = Cross(col1, col2);
Vec3d rhs = p;
Vec3d bary;
SolveLinearSystem(col1, col2, col3, rhs, bary);
intersection.lam1 = 0;
double eps = 1e-4;
if (bary.X() >= -eps && bary.Y() >= -eps && bary.X() + bary.Y() <= 1 + eps)
{
intersection.bary[0] = bary.X();
intersection.bary[1] = bary.Y();
intersection.bary[2] = 1.0 - bary.X() - bary.Y();
}
else
intersection.is_intersecting = false;
return intersection;
}
Intersection_ isIntersectingTrig (PointIndex pi_from, PointIndex pi_to, SurfaceElementIndex sei, double shift = 0.0)
{
return isIntersectingTrig(GetSeg(pi_from, pi_to), GetTrig(sei, shift));
}
void BuildSearchTree (double trig_shift)
{
static Timer t("BuildSearchTree");
RegionTimer rt(t);
Box<3> bbox(Box<3>::EMPTY_BOX);
for (PointIndex pi : mesh.Points().Range())
{
bbox.Add(mesh[pi]);
bbox.Add(GetPoint(pi, 1.1));
}
tree = make_unique<BoxTree<3>>(bbox);
for (auto sei : SurfaceElementsRange())
{
const auto& sel = Get(sei);
auto sel_index = sel.GetIndex();
Box<3> box(Box<3>::EMPTY_BOX);
for (auto pi : sel.PNums())
{
box.Add(GetPoint(pi, 0.));
box.Add(GetPoint(pi, trig_shift * GetLimit(pi)));
}
tree->Insert(box, sei);
}
}
template <typename TFunc>
void FindTreeIntersections (double trig_shift, double seg_shift, TFunc f, BitArray* relevant_points = nullptr)
{
static Timer t("GrowthVectorLimiter::FindTreeIntersections");
RegionTimer rt(t);
BuildSearchTree(trig_shift);
auto np_new = mesh.Points().Size();
int counter = 0;
for (auto i : IntRange(tool.np, np_new))
{
PointIndex pi_to = i + PointIndex::BASE;
PointIndex pi_from = map_from[pi_to];
if (!pi_from.IsValid())
throw Exception("Point not mapped");
if (relevant_points && !relevant_points->Test(pi_to) && !relevant_points->Test(pi_from))
continue;
Box<3> box(Box<3>::EMPTY_BOX);
auto seg = GetSeg(pi_to, seg_shift);
box.Add(GetPoint(pi_to, 0));
box.Add(GetPoint(pi_to, GetLimit(pi_from)));
tree->GetFirstIntersecting(box.PMin(), box.PMax(), [&] (SurfaceElementIndex sei) {
const auto& sel = Get(sei);
if (sel.PNums().Contains(pi_from))
return false;
if (sel.PNums().Contains(pi_to))
return false;
counter++;
f(pi_to, sei);
return false;
});
}
}
void FixIntersectingSurfaceTrigs ()
{
static Timer t("GrowthVectorLimiter::FixIntersectingSurfaceTrigs");
RegionTimer reg(t);
// check if surface trigs are intersecting each other
bool changed = true;
while (changed)
{
changed = false;
Point3d pmin, pmax;
mesh.GetBox(pmin, pmax);
BoxTree<3, SurfaceElementIndex> setree(pmin, pmax);
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));
box.Increase(1e-3 * box.Diam());
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();
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 (GetLimit(pi_max_limit) < 1e-10)
{
WriteErrorMesh("error_blayer_self_intersection_pi" + ToString(pi_max_limit) + ".vol.gz");
throw NgException("Stop meshing in boundary layer thickness limitation: overlapping regions detected at elements " + ToString(tri) + " and " + ToString(tri2));
}
if (debugparam.debugoutput && 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;
});
}
}
}
void LimitOriginalSurface (double safety)
{
static Timer t("GrowthVectorLimiter::LimitOriginalSurface");
RegionTimer reg(t);
PrintMessage(5, "GrowthVectorLimiter - original surface");
// limit to not intersect with other (original) surface elements
double trig_shift = 0;
double seg_shift = safety;
FindTreeIntersections(
trig_shift, seg_shift, [&] (PointIndex pi_to, SurfaceElementIndex sei) {
if (sei >= tool.nse)
return; // ignore new surface elements in first pass
LimitGrowthVector(pi_to, sei, trig_shift, seg_shift);
});
}
void LimitBoundaryLayer (double safety = 1.1)
{
static Timer t("GrowthVectorLimiter::LimitBoundaryLayer");
PrintMessage(5, "GrowthVectorLimiter - boundary layer");
// now limit again with shifted surface elements
double trig_shift = safety;
double seg_shift = safety;
size_t limit_counter = 1;
BitArray relevant_points, relevant_points_next;
relevant_points.SetSize(mesh.Points().Size() + 1);
relevant_points_next.SetSize(mesh.Points().Size() + 1);
relevant_points.Set();
while (limit_counter)
{
RegionTimer reg(t);
size_t find_counter = 0;
limit_counter = 0;
relevant_points_next.Clear();
FindTreeIntersections(
trig_shift, seg_shift, [&] (PointIndex pi_to, SurfaceElementIndex sei) {
find_counter++;
auto sel = Get(sei);
if (LimitGrowthVector(pi_to, sei, trig_shift, seg_shift))
{
limit_counter++;
relevant_points_next.SetBit(pi_to);
relevant_points_next.SetBit(map_from[pi_to]);
for (auto pi : sel.PNums())
{
relevant_points_next.SetBit(pi);
if (pi >= tool.np)
relevant_points_next.SetBit(map_from[pi]);
else
relevant_points_next.SetBit(map_from[pi]);
}
}
for (auto pi : sel.PNums())
{
if (pi >= tool.np)
return;
if (tool.mapto[pi].Size() == 0)
return;
}
if (LimitGrowthVector(pi_to, sei, trig_shift, seg_shift, true))
limit_counter++;
},
&relevant_points);
relevant_points = relevant_points_next;
}
}
void CheckLimits (int line)
{
auto check_point = [&] (PointIndex pi) {
if (limits[pi] < 1e-8)
{
WriteErrorMesh("error_blayer_intersection_pi" + ToString(pi) + ".vol.gz");
throw NgException(__FILE__ + ToString(line) + ": Stop meshing in boundary layer thickness limitation: overlapping regions detected at point " + ToString(pi));
}
};
for (auto pi : Range(growthvectors))
check_point(pi);
for (auto& [special_pi, special_point] : tool.special_boundary_points)
check_point(special_pi);
}
void Perform ()
{
limits.SetSize(mesh.Points().Size());
limits = 1.0;
if (tool.special_boundary_points.size())
{
auto point_to_sel = tool.mesh.CreatePoint2SurfaceElementTable();
for (auto& [pi, special_point] : tool.special_boundary_points)
{
auto maxh = mesh.GetH(mesh[pi]);
auto new_limit = min(0.3 * maxh / tool.total_height, 1.0);
if (new_limit < 1.0)
{
limits[pi] = new_limit;
for (auto sei : point_to_sel[pi])
for (auto pi_ : Get(sei).PNums())
limits[pi_] = new_limit;
}
}
}
std::array safeties = {0.5, 1.1, 1.5, 1.5};
// No smoothing in the last pass, to avoid generating new intersections
std::array smoothing_factors = {0.8, 0.7, 0.5, 0.0};
for (auto i_pass : Range(safeties.size()))
{
PrintMessage(4, "GrowthVectorLimiter pass ", i_pass);
double safety = safeties[i_pass];
CheckLimits(__LINE__);
// intersect segment with original surface elements
LimitOriginalSurface(2.1);
CheckLimits(__LINE__);
// intersect prisms with themself
LimitSelfIntersection(1.3 * safety);
CheckLimits(__LINE__);
// intesect segment with prism
LimitBoundaryLayer(safety);
CheckLimits(__LINE__);
for (auto i : Range(10))
EqualizeLimits(smoothing_factors[i_pass]);
CheckLimits(__LINE__);
if (i_pass == safeties.size() - 1)
FixIntersectingSurfaceTrigs();
CheckLimits(__LINE__);
}
for (auto i : Range(growthvectors))
growthvectors[i] *= limits[i];
for (auto& [special_pi, special_point] : tool.special_boundary_points)
{
for (auto& group : special_point.growth_groups)
{
group.growth_vector *= limits[special_pi];
}
}
}
};
} // namespace netgen

View File

@ -2,7 +2,7 @@
namespace netgen
{
unique_ptr<Mesh> GetOpenElements( const Mesh & m, int dom = 0 )
unique_ptr<Mesh> GetOpenElements( const Mesh & m, int dom = 0, bool only_quads = false )
{
static Timer t("GetOpenElements"); RegionTimer rt(t);
auto mesh = make_unique<Mesh>();
@ -40,7 +40,8 @@ namespace netgen
mesh->ClearSurfaceElements();
for (auto & el : openelements)
mesh->AddSurfaceElement( el );
if(!only_quads || el.GetNP() == 4)
mesh->AddSurfaceElement( el );
mesh->Compress();
return mesh;

View File

@ -3,7 +3,7 @@
namespace netgen
{
unique_ptr<Mesh> GetOpenElements( const Mesh & m, int dom = 0 );
unique_ptr<Mesh> GetOpenElements( const Mesh & m, int dom = 0, bool only_quads = false );
unique_ptr<Mesh> FilterMesh( const Mesh & m, FlatArray<PointIndex> points, FlatArray<SurfaceElementIndex> sels = Array<SurfaceElementIndex>{}, FlatArray<ElementIndex> els = Array<ElementIndex>{} );

View File

@ -3401,7 +3401,7 @@ namespace netgen
double Mesh :: MaxHDomain (int dom) const
{
if (maxhdomain.Size())
if (dom >= 0 && dom < maxhdomain.Size())
return maxhdomain.Get(dom);
else
return 1e10;

View File

@ -3,6 +3,7 @@
#include <mystdlib.h>
#include "meshing.hpp"
#include "debugging.hpp"
#include "boundarylayer.hpp"
namespace netgen
{
@ -372,6 +373,7 @@ namespace netgen
if(debugparam.write_mesh_on_error) {
md.mesh->Save("open_quads_starting_mesh_"+ToString(md.domain)+".vol.gz");
GetOpenElements(*md.mesh, md.domain)->Save("open_quads_rest_" + ToString(md.domain)+".vol.gz");
GetOpenElements(*md.mesh, md.domain, true)->Save("open_quads_rest_" + ToString(md.domain)+"_only_quads.vol.gz");
}
PrintSysError ("mesh has still open quads");
throw NgException ("Stop meshing since too many attempts");
@ -431,7 +433,7 @@ namespace netgen
mesh.FindOpenElements(domain);
PrintMessage (5, mesh.GetNOpenElements(), " open faces");
// GetOpenElements( mesh, domain )->Save("open_"+ToString(cntsteps)+".vol");
// GetOpenElements( mesh, domain )->Save("open_"+ToString(domain)+"_"+ToString(cntsteps)+".vol");
cntsteps++;
@ -609,12 +611,18 @@ namespace netgen
static Timer t("MeshVolume"); RegionTimer reg(t);
mesh3d.Compress();
for (auto bl : mp.boundary_layers)
GenerateBoundaryLayer(mesh3d, bl);
if(mesh3d.GetNDomains()==0)
return MESHING3_OK;
auto geo = mesh3d.GetGeometry();
for (auto i : Range(std::min(geo->GetNSolids(), (size_t)mesh3d.GetNDomains())))
if (auto name = geo->GetSolid(i).properties.name)
mesh3d.SetMaterial (i+1, *name);
for (auto bl : mp.boundary_layers)
GenerateBoundaryLayer(mesh3d, bl);
if (!mesh3d.HasLocalHFunction())
mesh3d.CalcLocalH(mp.grading);
@ -639,6 +647,8 @@ namespace netgen
MeshDomain(md[i]);
}
catch (const Exception & e) {
if(debugparam.write_mesh_on_error)
md[i].mesh->Save("meshing_error_domain_"+ToString(md[i].domain)+".vol.gz");
cerr << "Meshing of domain " << i+1 << " failed with error: " << e.what() << endl;
throw e;
}
@ -754,7 +764,7 @@ namespace netgen
auto geo = mesh.GetGeometry();
if(!geo) return;
auto n_solids = geo->GetNSolids();
if(!n_solids) return;
if(n_solids < domain) return;
if(geo->GetSolid(domain-1).free_edges.Size() == 0)
return;

View File

@ -53,7 +53,6 @@ namespace netgen
#include "bisect.hpp"
#include "hprefinement.hpp"
#include "boundarylayer.hpp"
#include "specials.hpp"
#include "validate.hpp"
#include "basegeom.hpp"

View File

@ -2922,14 +2922,17 @@ namespace netgen
case 1: print_vec(std::get<1>(mp.thickness)); break;
}
ost <<"\n new_material: ";
switch(mp.new_material.index())
{
case 0: ost << std::get<0>(mp.new_material); break;
case 1:
for (const auto & [key, value] : std::get<1>(mp.new_material))
ost << key << " -> " << value << ", ";
break;
}
if(!mp.new_material) ost << "nullopt";
else {
switch(mp.new_material->index())
{
case 0: ost << std::get<0>(*mp.new_material); break;
case 1:
for (const auto & [key, value] : std::get<1>(*mp.new_material))
ost << key << " -> " << value << ", ";
break;
}
}
ost << "\n domain: ";
switch(mp.domain.index())
{
@ -2951,9 +2954,8 @@ namespace netgen
ost << "nullopt";
ost << "\n grow_edges: " << mp.grow_edges;
ost << "\n limit_growth_vectors: " << mp.limit_growth_vectors;
ost << "\n sides_keep_surfaceindex: " << mp.sides_keep_surfaceindex;
ost << "\n keep_surfaceindex: " << mp.keep_surfaceindex;
ost << "\n limit_safety: " << mp.limit_safety;
ost << "\n sides_keep_surfaceindex: " << (mp.sides_keep_surfaceindex ? ToString(*mp.sides_keep_surfaceindex) : "nullopt");
ost << "\n disable_curving: " << mp.disable_curving;
ost << endl;
return ost;
}

View File

@ -1369,18 +1369,16 @@ namespace netgen
struct BoundaryLayerParameters
{
std::variant<string, int, std::vector<int>> boundary;
std::variant<double, std::vector<double>> thickness;
std::variant<string, std::map<string, string>> new_material;
std::variant<string, int, std::vector<int>> domain;
bool outside;
std::optional<std::variant<string, std::vector<int>>> project_boundaries;
bool grow_edges;
bool limit_growth_vectors;
bool sides_keep_surfaceindex;
bool keep_surfaceindex;
double limit_safety = 0.3; // alloow only 30% of the growth vector length
std::variant<string, int, std::vector<int>> boundary = ".*";
std::optional<std::variant<string, std::map<string, string>>> new_material = nullopt;
std::optional<std::variant<string, std::vector<int>>> project_boundaries = nullopt;
bool outside = false;
bool grow_edges = true;
bool limit_growth_vectors = false; // automatic reduction of layer thickness to avoid intersections
std::optional<bool> sides_keep_surfaceindex = nullopt; // !outside by default
bool disable_curving = true; // disable curving affected boundaries/edges (could lead to self-intersecting volume elements)
};

View File

@ -8,6 +8,7 @@
#include <mystdlib.h>
#include "meshing.hpp"
#include "boundarylayer.hpp"
// #include <csg.hpp>
// #include <geometry2d.hpp>
#include <../interface/rw_medit.hpp>
@ -1471,12 +1472,12 @@ py::arg("point_tolerance") = -1.)
.def ("BoundaryLayer2", GenerateBoundaryLayer2, py::arg("domain"), py::arg("thicknesses"), py::arg("make_new_domain")=true, py::arg("boundaries")=Array<int>{})
.def ("BoundaryLayer", [](Mesh & self, variant<string, int, std::vector<int>> boundary,
variant<double, std::vector<double>> thickness,
variant<string, map<string, string>> material,
optional<variant<string, map<string, string>>> material,
variant<string, int, std::vector<int>> domain, bool outside,
optional<variant<string, std::vector<int>>> project_boundaries,
bool grow_edges, bool limit_growth_vectors,
bool sides_keep_surfaceindex,
bool keep_surfaceindex)
bool disable_curving)
{
BoundaryLayerParameters blp;
blp.boundary = boundary;
@ -1488,13 +1489,13 @@ py::arg("point_tolerance") = -1.)
blp.grow_edges = grow_edges;
blp.limit_growth_vectors = limit_growth_vectors;
blp.sides_keep_surfaceindex = sides_keep_surfaceindex;
blp.keep_surfaceindex = keep_surfaceindex;
blp.disable_curving = disable_curving;
GenerateBoundaryLayer (self, blp);
self.UpdateTopology();
}, py::arg("boundary"), py::arg("thickness"), py::arg("material"),
}, py::arg("boundary"), py::arg("thickness"), py::arg("material")=nullopt,
py::arg("domains") = ".*", py::arg("outside") = false,
py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true, py::arg("limit_growth_vectors") = true, py::arg("sides_keep_surfaceindex")=false,
py::arg("keep_surfaceindex")=false, "Add boundary layer to mesh. see help(BoundaryLayerParameters) for details.")
py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true, py::arg("limit_growth_vectors") = false, py::arg("sides_keep_surfaceindex")=false,
py::arg("disable_curving")=true, "Add boundary layer to mesh. see help(BoundaryLayerParameters) for details.")
.def_static ("EnableTableClass", [] (string name, bool set)
{
@ -1724,15 +1725,14 @@ py::arg("point_tolerance") = -1.)
.def(py::init([](
std::variant<string, int, std::vector<int>> boundary,
std::variant<double, std::vector<double>> thickness,
std::variant<string, std::map<string, string>> new_material,
std::optional<std::variant<string, std::map<string, string>>> new_material,
std::variant<string, int, std::vector<int>> domain,
bool outside,
std::optional<std::variant<string, std::vector<int>>> project_boundaries,
bool grow_edges,
bool limit_growth_vectors,
bool sides_keep_surfaceindex,
bool keep_surfaceindex,
double limit_safety)
std::optional<bool> sides_keep_surfaceindex,
bool disable_curving)
{
BoundaryLayerParameters blp;
blp.boundary = boundary;
@ -1744,15 +1744,14 @@ py::arg("point_tolerance") = -1.)
blp.grow_edges = grow_edges;
blp.limit_growth_vectors = limit_growth_vectors;
blp.sides_keep_surfaceindex = sides_keep_surfaceindex;
blp.keep_surfaceindex = keep_surfaceindex;
blp.limit_safety = limit_safety;
blp.disable_curving = disable_curving;
return blp;
}),
py::arg("boundary"), py::arg("thickness"), py::arg("new_material"),
py::arg("boundary"), py::arg("thickness"), py::arg("new_material")=nullopt,
py::arg("domain") = ".*", py::arg("outside") = false,
py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true,
py::arg("limit_growth_vectors") = true, py::arg("sides_keep_surfaceindex")=false,
py::arg("keep_surfaceindex")=false, py::arg("limit_safety")=0.3,
py::arg("limit_growth_vectors") = false, py::arg("sides_keep_surfaceindex")=nullopt,
py::arg("disable_curving")=true,
R"delimiter(
Add boundary layer to mesh.
@ -1804,8 +1803,7 @@ project_boundaries : Optional[str] = None
.def_readwrite("grow_edges", &BoundaryLayerParameters::grow_edges)
.def_readwrite("limit_growth_vectors", &BoundaryLayerParameters::limit_growth_vectors)
.def_readwrite("sides_keep_surfaceindex", &BoundaryLayerParameters::sides_keep_surfaceindex)
.def_readwrite("keep_surfaceindex", &BoundaryLayerParameters::keep_surfaceindex)
.def_readwrite("limit_safety", &BoundaryLayerParameters::limit_safety)
.def_readwrite("disable_curving", &BoundaryLayerParameters::disable_curving)
;
py::implicitly_convertible<py::dict, BoundaryLayerParameters>();

View File

@ -9,6 +9,7 @@ The interface between the GUI and the netgen library
#include <linalg.hpp>
#include <meshing.hpp>
#include "../libsrc/meshing/boundarylayer.hpp"
#include <inctcl.hpp>

View File

@ -1,6 +1,7 @@
import pytest
from netgen.csg import *
from netgen.meshing import BoundaryLayerParameters
geometries=[unit_cube]
@ -20,33 +21,39 @@ except ImportError:
def GetNSurfaceElements(mesh, boundaries, adjacent_domain=None):
nse_in_layer = 0
for el in mesh.Elements2D():
if mesh.GetBCName(el.index-1) in boundaries:
if len(el.vertices)==3 and mesh.GetBCName(el.index-1) in boundaries:
if adjacent_domain is None:
print("add el", el.vertices)
nse_in_layer += 1
else:
if (mesh.FaceDescriptor(el.index).domin > 0 and mesh.GetMaterial(mesh.FaceDescriptor(el.index).domin) == adjacent_domain) or (mesh.FaceDescriptor(el.index).domout > 0 and mesh.GetMaterial(mesh.FaceDescriptor(el.index).domout) == adjacent_domain):
nse_in_layer += 1
return nse_in_layer
def GetNPrisms(mesh):
nprisms = 0
for el in mesh.Elements3D():
if len(el.vertices) == 6:
nprisms += 1
return nprisms
@pytest.mark.parametrize("outside", [True, False])
@pytest.mark.parametrize("geo", geometries)
def test_boundarylayer(outside, geo, capfd):
mesh = geo.GenerateMesh(maxh=0.3)
ne_before = mesh.ne
layer_surfacenames = ["right", "top", "left", "back", "bottom"]
mesh.BoundaryLayer("|".join(layer_surfacenames), [0.01, 0.01], "layer", outside=outside)
should_ne = ne_before + 2 * GetNSurfaceElements(mesh, layer_surfacenames)
assert mesh.ne == should_ne
blayer_params = [BoundaryLayerParameters('|'.join(layer_surfacenames), [0.01, 0.01], "layer", outside=outside)]
mesh = geo.GenerateMesh(maxh=0.3, boundary_layers=blayer_params)
should_n_prisms = 2 * GetNSurfaceElements(mesh, layer_surfacenames)
assert GetNPrisms(mesh) == should_n_prisms
capture = capfd.readouterr()
assert not "elements are not matching" in capture.out
for side in ["front"]:
mesh.BoundaryLayer(side, [0.001, 0.001], "layer", outside=outside)
should_ne += 2 * GetNSurfaceElements(mesh, [side])
assert mesh.ne == should_ne
capture = capfd.readouterr()
assert not "elements are not matching" in capture.out
blayer_params.append(BoundaryLayerParameters("front", [0.01, 0.01], "layer", outside=True)) # outside=outside not working...
mesh = geo.GenerateMesh(maxh=0.3, boundary_layers=blayer_params)
should_n_prisms += 2 * GetNSurfaceElements(mesh, ["front"])
assert GetNPrisms(mesh) == should_n_prisms
capture = capfd.readouterr()
assert not "elements are not matching" in capture.out
@pytest.mark.parametrize("outside", [True, False])
@pytest.mark.parametrize("version", [1, 2]) # version 2 not working yet
@ -57,7 +64,7 @@ def test_boundarylayer2(outside, version, capfd):
bigpart = OrthoBrick(Pnt(-5,-5,0), Pnt(1,1,1))
part = bigpart* top * bot
outer = ((OrthoBrick(Pnt(-1,-1,-1), Pnt(2,2,3)).bc("outer") * Plane(Pnt(2,2,2), Vec(0,0,1)).bc("outertop")))
geo.Add((part * outer).mat("part"))
if version == 1:
geo.Add((outer-part).mat("rest"))
@ -67,11 +74,11 @@ def test_boundarylayer2(outside, version, capfd):
geo.Add((outer*top*bot-bigpart).mat("rest"))
geo.CloseSurfaces(top, bot, [])
mesh = geo.GenerateMesh()
should_ne = mesh.ne + 2 * GetNSurfaceElements(mesh, ["default"], "part")
layersize = 0.025
mesh.BoundaryLayer("default", [layersize, layersize], "part", domains="part", outside=outside)
assert mesh.ne == should_ne
mesh = geo.GenerateMesh(boundary_layers=[BoundaryLayerParameters("default", [layersize, layersize], "part", domain="part", outside=outside)])
should_n_prisms = 2 * GetNSurfaceElements(mesh, ["default"], "part")
# assert GetNPrisms(mesh) == should_n_prisms
assert not "elements are not matching" in capfd.readouterr().out
# import netgen.gui
ngs = pytest.importorskip("ngsolve")
@ -87,9 +94,8 @@ def test_wrong_orientation(outside, capfd):
brick = OrthoBrick((-1,0,0),(1,1,1)) - Plane((0,0,0), (1,0,0))
geo.Add(brick.mat("air"))
mesh = geo.GenerateMesh()
mesh = geo.GenerateMesh(boundary_layers=[BoundaryLayerParameters(".*", [0.1], "air", domain="air", outside=outside)])
mesh.BoundaryLayer(".*", 0.1, "air", domains="air", outside=outside)
ngs = pytest.importorskip("ngsolve")
mesh = ngs.Mesh(mesh)
assert ngs.Integrate(1, mesh) == pytest.approx(1.2**3 if outside else 1)
@ -102,8 +108,7 @@ def test_splitted_surface():
geo.Add((brick-slots).mat("block"))
geo.Add((brick*slots).mat("slot"))
mesh = geo.GenerateMesh()
mesh.BoundaryLayer(".*", [0.001, 0.001], "block", "block", outside=False)
mesh = geo.GenerateMesh(boundary_layers=[BoundaryLayerParameters(".*", [0.001, 0.001], "block", domain="block", outside=False)])
ngs = pytest.importorskip("ngsolve")
mesh = ngs.Mesh(mesh)
assert ngs.Integrate(1, mesh) == pytest.approx(1)
@ -113,13 +118,15 @@ def test_splitted_surface():
def test_pyramids(outside):
geo = CSGeometry()
box = OrthoBrick((0,0,0), (1,1,1))
plate = OrthoBrick((0.3,0.3,0.4),(0.7,0.7,1)) * Plane((0,0,0.6), (0,0,1)).bc("top")
dist = 0.3
plate = OrthoBrick((dist,dist,0.4),(1-dist,1-dist,1)) * Plane((0,0,0.6), (0,0,1)).bc("top")
geo.Add((box-plate).mat("air"))
geo.Add(plate.mat("plate"))
mesh = geo.GenerateMesh()
mesh.BoundaryLayer("top", [0.01], "layer", "plate", outside=outside)
blayers = [BoundaryLayerParameters("top", [0.01], "layer", domain="plate", outside=outside)]
mesh = geo.GenerateMesh(boundary_layers=blayers)
ngs = pytest.importorskip("ngsolve")
mesh = ngs.Mesh(mesh)
assert ngs.Integrate(1, mesh.Materials("plate")) == pytest.approx(0.032 if outside else 0.0304)
assert ngs.Integrate(1, mesh.Materials("layer")) == pytest.approx(0.0016)
assert ngs.Integrate(1, mesh.Materials("air")) == pytest.approx(0.9664 if outside else 0.968)
@ -167,8 +174,8 @@ def _test_with_inner_corner(outside, capfd):
geo.Add(coil1, col=(0.72, 0.45, 0.2))
geo.Add(coil2, col=(0.72, 0.45, 0.2))
geo.Add(oil.mat("oil"), transparent=True)
mesh = geo.GenerateMesh()
mesh.BoundaryLayer("core_front", [0.001, 0.002], "core", "core", outside=outside)
blayers = [BoundaryLayerParameters("core_front", [0.001, 0.002], "core", domain="core", outside=outside)]
mesh = geo.GenerateMesh(boundary_layers = blayers)
ngs = pytest.importorskip("ngsolve")
mesh = ngs.Mesh(mesh)
capture = capfd.readouterr()