some more

This commit is contained in:
Matthias Hochsteger 2024-10-04 20:38:58 +02:00
parent 9693a42a6a
commit dbffb79550
7 changed files with 179 additions and 79 deletions

View File

@ -961,10 +961,18 @@ namespace ngcore
f.precision(4);
f << R"CODE_(
<head>
<script src="https://d3js.org/d3.v5.min.js"></script>
<script src="https://cdn.jsdelivr.net/npm/d3@7"></script>
<script src="https://unpkg.com/sunburst-chart"></script>
<style>body { margin: 0 }</style>
<style>
body { margin: 0 }
.tooltip {
white-space: pre-line !important;
max-width: 800px !important;
word-wrap: break-word !important;
padding: 10px !important;
}
</style>
)CODE_";
if(!time_or_memory)
f << "<title>Maximum Memory Consumption</title>\n";

View File

@ -556,11 +556,12 @@ void BoundaryLayerTool ::InsertNewElements(
for (auto i : Range(par_heights)) {
height += par_heights[i];
auto pi_new = mesh.AddPoint(p);
// mesh.AddLockedPoint(pi_new);
mapfrom.Append(pi);
new_points.Append(pi_new);
growth_vector_map[pi_new] = {&growth_vector, height};
if (special_boundary_points.count(pi) > 0)
mesh.AddLockedPoint(pi_new);
// if (special_boundary_points.count(pi) > 0)
// mesh.AddLockedPoint(pi_new);
pi_last = pi_new;
}
};
@ -680,7 +681,8 @@ void BoundaryLayerTool ::InsertNewElements(
s0.edgenr = segj.edgenr;
s0.si = segj.si;
new_segments.Append(s0);
if(type==3) new_segments_on_moved_bnd.Append(s0);
if (type == 3)
new_segments_on_moved_bnd.Append(s0);
for (auto i : Range(par_heights)) {
Element2d sel(QUAD);
@ -726,25 +728,17 @@ void BoundaryLayerTool ::InsertNewElements(
s3.edgenr = getEdgeNr(segj.edgenr);
s3.si = segj.si;
new_segments.Append(s3);
if(type==3) new_segments_on_moved_bnd.Append(s0);
}
else if (type == 1) {
if (type == 3)
new_segments_on_moved_bnd.Append(s0);
} else if (type == 3) {
PointIndex pp1 = segj[1];
PointIndex pp2 = segj[0];
if (in_surface_direction.Test(segj.si)) {
if (!in_surface_direction.Test(segj.si)) {
Swap(pp1, pp2);
}
PointIndex p1 = pp1;
PointIndex p2 = pp2;
PointIndex p3, p4;
Segment s0;
s0[0] = p1;
s0[1] = p2;
s0[2] = PointIndex::INVALID;
s0.edgenr = segj.edgenr;
s0.si = segj.si;
new_segments.Append(s0);
if(type==3) new_segments_on_moved_bnd.Append(s0);
for (auto i : Range(par_heights)) {
Element2d sel(QUAD);
@ -761,6 +755,9 @@ void BoundaryLayerTool ::InsertNewElements(
sel.SetIndex(si_map[segj.si]);
new_sels.Append(sel);
new_sels_on_moved_bnd.Append(sel);
p1 = p4;
p2 = p3;
}
}
}
}
@ -841,16 +838,15 @@ void BoundaryLayerTool ::InsertNewElements(
}
}
// for (SegmentIndex sei = 0; sei < nseg; sei++) {
// auto &seg = segments[sei];
// if (is_boundary_moved.Test(seg.si)) {
// if(insert_only_volume_elements) throw
// Exception("insert_only_volume_elements and is_boundary_moved are
// incompatible"); for (auto &p : seg.PNums())
// if (hasMoved(p))
// p = newPoint(p);
// }
// }
for (SegmentIndex sei = 0; sei < nseg; sei++) {
auto &seg = segments[sei];
if (is_boundary_moved.Test(seg.si)) {
// cout << "moved setg " << seg << endl;
for (auto &p : seg.PNums())
if (hasMoved(p))
p = newPoint(p);
}
}
// fill holes in surface mesh at special boundary points (i.e. points with >=4
// adjacent boundary faces)
@ -987,8 +983,9 @@ void BoundaryLayerTool ::SetDomInOutSides() {
}
void BoundaryLayerTool ::AddSegments() {
auto & new_segs = insert_only_volume_elements ? new_segments_on_moved_bnd : new_segments;
cout << "add new segs " << endl << new_segs << endl;
auto &new_segs =
insert_only_volume_elements ? new_segments_on_moved_bnd : new_segments;
// cout << "add new segs " << endl << new_segs << endl;
if (have_single_segments)
MergeAndAddSegments(mesh, segments, new_segs);
else {
@ -1153,13 +1150,19 @@ void BoundaryLayerTool ::Perform() {
SetDomInOut();
AddSegments();
AddSurfaceElements();
mesh.CalcSurfacesOfNode();
topo.SetBuildVertex2Element(true);
mesh.UpdateTopology();
// cout << "GW 19911 " << growthvectors[19911] << endl;
InterpolateGrowthVectors();
// cout << "GW 19911 " << growthvectors[19911] << endl;
InterpolateSurfaceGrowthVectors();
// cout << "GW 19911 " << growthvectors[19911] << endl;
AddSurfaceElements();
if (params.limit_growth_vectors)
LimitGrowthVectorLengths();

View File

@ -158,22 +158,34 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
auto np = mesh.GetNP();
non_bl_growth_vectors.clear();
// for(const auto & sel : new_sels) {
// for(auto pi : sel.PNums()) {
// if(mesh[pi].Type() == INNERPOINT)
// mesh[pi].SetType(SURFACEPOINT);
// }
// }
// cout << __FILE__ << ":" << __LINE__ << endl;
auto getGW = [&](PointIndex pi) -> Vec<3> {
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);
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;
growthvectors[pi] += vec;
// // cout << "add gw " << pi << "\t" << vec << endl;
// if (growth_vector_map.count(pi) == 0) {
// // cout << "\t make new entry" << endl;
// non_bl_growth_vectors[pi] = .0;
// growth_vector_map[pi] = {&non_bl_growth_vectors[pi], 1.0};
// }
// auto [gw, height] = growth_vector_map[pi];
// // cout << "\tcurrent gw " << *gw << "\t" << height << endl;
// *gw += 1.0 / height * vec;
};
auto hasMoved = [&](PointIndex pi) {
@ -181,38 +193,47 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
special_boundary_points.count(pi);
};
// cout << __FILE__ << ":" << __LINE__ << endl;
std::set<PointIndex> points_set;
for (SurfaceElementIndex sei : mesh.SurfaceElements().Range()) {
for (auto pi : mesh[sei].PNums()) {
auto pi_from = mapfrom[pi];
if ((pi_from.IsValid() && mesh[pi_from].Type() == SURFACEPOINT) ||
(!pi_from.IsValid() && mapto[pi].Size() == 0 &&
mesh[pi].Type() == SURFACEPOINT))
for (const auto &sel: mesh.SurfaceElements())
for (auto pi : sel.PNums())
if(mesh[pi].Type() == SURFACEPOINT && hasMoved(pi))
points_set.insert(pi);
}
}
Array<bool> has_moved_points(max_edge_nr + 1);
has_moved_points = false;
std::set<PointIndex> moved_edge_points;
// Array<bool> has_moved_points(max_edge_nr + 1);
// has_moved_points = false;
// std::set<PointIndex> moved_edge_points;
for (auto seg : segments) {
if (hasMoved(seg[0]) != hasMoved(seg[1]))
has_moved_points[seg.edgenr] = true;
}
// for (auto seg : segments) {
// if (hasMoved(seg[0]) != hasMoved(seg[1]))
// has_moved_points[seg.edgenr] = true;
// }
// cout << __FILE__ << ":" << __LINE__ << endl;
for (auto seg : segments)
if (has_moved_points[seg.edgenr])
for (auto pi : seg.PNums())
if (mesh[pi].Type() == EDGEPOINT)
points_set.insert(pi);
// for (auto seg : segments)
// if (has_moved_points[seg.edgenr])
// for (auto pi : seg.PNums())
// if (mesh[pi].Type() == EDGEPOINT)
// points_set.insert(pi);
Array<PointIndex> points;
for (auto pi : points_set)
points.Append(pi);
QuickSort(points);
// cout << __FILE__ << ":" << __LINE__ << endl;
// cout << "points to interpolate " << endl << points << endl;
// cout << __FILE__ << ":" << __LINE__ << endl;
auto p2sel = mesh.CreatePoint2SurfaceElementTable();
// cout << __FILE__ << ":" << __LINE__ << endl;
// auto p2sel = ngcore::CreateSortedTable<SurfaceElementIndex, PointIndex>(
// new_sels.Range(),
// [&](auto &table, SurfaceElementIndex ei) {
// for (PointIndex pi : new_sels[ei].PNums())
// table.Add(pi, ei);
// },
// mesh.GetNP());
// smooth tangential part of growth vectors from edges to surface elements
Array<Vec<3>, PointIndex> corrections(mesh.GetNP());
corrections = 0.0;
@ -223,6 +244,7 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
double weight;
};
Array<ArrayMem<Neighbor, 20>> neighbors(points.Size());
// cout << __FILE__ << ":" << __LINE__ << endl;
ArrayMem<double, 20> angles;
ArrayMem<double, 20> inv_dists;
@ -244,6 +266,9 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
}
}
p_neighbors.Append({pi1, sei, 0.0});
// if((mesh[pi1]-mesh[pi]).Length() < 1e-10) {
// cout << "close points " << pi << "\t" << pi1 << endl;
// }
inv_dists.Append(1.0 / (mesh[pi1] - mesh[pi]).Length());
auto dot = (mesh[pi1] - mesh[pi]).Normalize() *
(mesh[pi2] - mesh[pi]).Normalize();
@ -257,6 +282,9 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
for (auto angle : angles)
sum_angle += angle;
// cout << "angles " << angles << endl;
// cout << "inv_dists " << inv_dists << endl;
double sum_weight = 0.0;
for (auto i : Range(inv_dists)) {
p_neighbors[i].weight =
@ -265,16 +293,38 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
}
for (auto i : Range(inv_dists))
p_neighbors[i].weight /= sum_weight;
// if(pi == 19911) {
// cout << "pi " << pi << endl;
// for(auto & nb : p_neighbors) {
// cout << "neighbor " << nb.pi << "\t" << nb.weight << endl;
// }
// cout << "inv_dist " << inv_dists << endl;
// cout << "angles " << angles << endl;
// }
}
// cout << __FILE__ << ":" << __LINE__ << endl;
Array<Vec<3>, SurfaceElementIndex> surf_normals(mesh.GetNSE());
for (auto sei : mesh.SurfaceElements().Range())
surf_normals[sei] = getNormal(mesh[sei]);
// cout << __FILE__ << ":" << __LINE__ << endl;
BitArray interpolate_tangent(mesh.GetNP()+1);
// cout << __FILE__ << ":" << __LINE__ << endl;
interpolate_tangent = false;
for(auto pi : points) {
for (auto sei : p2sel[pi])
if(is_boundary_moved[mesh[sei].GetIndex()])
interpolate_tangent.SetBit(pi);
}
// cout << __FILE__ << ":" << __LINE__ << endl;
constexpr int N_STEPS = 64;
for ([[maybe_unused]] auto i : Range(N_STEPS)) {
for (auto i : points.Range()) {
auto pi = points[i];
// cout << "AVERAGE " << pi << endl;
auto &p_neighbors = neighbors[i];
ArrayMem<Vec<3>, 20> g_vectors;
@ -283,17 +333,21 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
// average only tangent component on new bl points, average whole growth
// vector otherwise
bool do_average_tangent = mapfrom[pi].IsValid();
bool do_average_tangent = true;
for (const auto &s : p_neighbors) {
auto gw_other = getGW(s.pi) + corrections[s.pi];
// if(pi == 19911) cout << "neighbor " << s.pi << "\t" << s.weight << '\t' << gw_other << "\tdo avg: " << do_average_tangent << endl;
// if(pi == 19911) cout << "\tneighbor gw" << gw_other << endl;
if (do_average_tangent) {
auto n = surf_normals[s.sei];
gw_other = gw_other - (gw_other * n) * n;
}
// if(pi == 19911) cout << "\tneighbor gw" << gw_other << endl;
auto v = gw_other;
auto len = v.Length2();
sum_len += len;
max_len = max(max_len, len);
// if(pi == 19911) cout << "\tneighbor v" << v << endl;
g_vectors.Append(v);
}
@ -311,11 +365,13 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
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;
}
}

View File

@ -1,3 +1,4 @@
#include <core/array.hpp>
#include "boundarylayer.hpp"
namespace netgen {
@ -207,33 +208,41 @@ struct GrowthVectorLimiter {
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];
// if(pi_from == 21110) cout << "equalize " << pi << "\tfactor " << factor << endl;
std::set<PointIndex> pis;
for (auto sei : p2sel[pi])
for (auto pi_ : Get(sei).PNums())
for (auto pi_ : tool.new_sels[sei].PNums())
pis.insert(pi_);
ArrayMem<double, 20> limits;
for (auto pi : pis) {
auto limit = GetLimit(pi);
for (auto pi1 : pis) {
auto limit = GetLimit(pi1);
if (limit > 0.0)
limits.Append(GetLimit(pi));
limits.Append(GetLimit(pi1));
// if(pi_from == 21110) cout << "\tneighbor point " << map_from[pi1] << " -> " << pi1 << " with limit " << limit << endl;
}
// if(pi_from == 21110) cout << "\town limit " << GetLimit(pi) << endl;
if (limits.Size() == 0)
continue;
QuickSort(limits);
double mean_limit = limits[limits.Size() / 2];
// if(pi_from == 21110) cout << "\tmean limit " << mean_limit << endl;
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(pi_from == 21110) cout << "\tnew limit " << GetLimit(pi) << endl;
}
}
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) {
@ -342,6 +351,7 @@ struct GrowthVectorLimiter {
}
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]);
@ -487,6 +497,7 @@ struct GrowthVectorLimiter {
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;
@ -500,6 +511,7 @@ struct GrowthVectorLimiter {
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;
@ -531,12 +543,13 @@ struct GrowthVectorLimiter {
limits.SetSize(mesh.Points().Size());
limits = 1.0;
std::array safeties = {0.5, 0.8, 1.1, 1.1};
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];
LimitOriginalSurface(2.1);
@ -550,6 +563,18 @@ struct GrowthVectorLimiter {
FixIntersectingSurfaceTrigs();
}
// cout << "limits " << endl << limits << endl;
// double avg_limit = 0;
// for(auto pi : limits.Range()) {
// avg_limit += limits[pi];
// }
// avg_limit = avg_limit / limits.Size();
// for(auto pi : limits.Range()) {
// if(limits[pi] < 0.01 * avg_limit)
// cout << "small limit " << pi << "\t" << limits[pi] << endl;
// }
for (auto i : Range(growthvectors))
growthvectors[i] *= limits[i];

View File

@ -1320,7 +1320,7 @@ void MeshOptimize3d :: SwapImprove (const NgBitArray * working_elements)
for (ElementIndex eli : myrange)
{
const auto & el = mesh[eli];
if(el.Flags().fixed)
if(el.Flags().fixed || el.GetType() != TET)
continue;
if(mp.only3D_domain_nr && mp.only3D_domain_nr != el.GetIndex())
@ -1381,12 +1381,20 @@ void MeshOptimize3d :: SwapImprove (const NgBitArray * working_elements)
break;
auto [pi0, pi1] = edges[i];
try {
double d_badness = SwapImproveEdge (working_elements, elementsonnode, faces, pi0, pi1, true);
if(d_badness<0.0)
{
int index = improvement_counter++;
candidate_edges[index] = make_tuple(d_badness, i);
}
}
catch(const ngcore::Exception &e) {
cerr << "Error in SwapImproveEdge " << pi0 << '-' << pi1 << endl;
if(debugparam.write_mesh_on_error)
mesh.Save("error_swapimproveedge_" +ToString(pi0) + "_" + ToString(pi1) +".vol");
throw;
}
}
}, TasksPerThread (4));
@ -2649,6 +2657,7 @@ double MeshOptimize3d :: SplitImprove2Element (
Element & elem = mesh[ei1];
if (elem.IsDeleted()) return false;
if (ei1 == ei) continue;
if (elem.GetType() != TET) return false;
if (elem[0] == pi3 || elem[1] == pi3 || elem[2] == pi3 || elem[3] == pi3 || (elem.GetNP()==5 && elem[4]==pi3))
if(!has_both_points1.Contains(ei1))

View File

@ -499,6 +499,7 @@ namespace netgen
PrintMessage (3, "Call remove problem");
RemoveProblem (mesh, domain);
mesh.FindOpenElements(domain);
GetOpenElements( mesh, domain )->Save("open_"+ToString(domain)+"_"+ToString(cntsteps)+".vol");
}
else
{
@ -546,14 +547,12 @@ namespace netgen
{
auto first_new_pi = m_.pmap.Range().Next();
auto & m = *m_.mesh;
Array<PointIndex, PointIndex> pmap(m.Points().Size());
for(auto pi : Range(PointIndex(PointIndex::BASE), first_new_pi))
pmap[pi] = m_.pmap[pi];
Array<PointIndex, PointIndex> pmap;
pmap = m_.pmap;
pmap.SetSize(mesh.GetNP() + m.GetNP() - first_new_pi);
for (auto pi : Range(first_new_pi, m.Points().Range().Next()))
pmap[pi] = mesh.AddPoint(m[pi]);
for ( auto el : m.VolumeElements() )
{
for (auto i : Range(el.GetNP()))
@ -590,7 +589,7 @@ namespace netgen
// extern double teterrpow;
MESHING3_RESULT MeshVolume (const MeshingParameters & mp, Mesh& mesh3d)
{
static Timer t("MeshVolume"); RegionTimer reg(t);
static Timer t("MeshVolume1"); RegionTimer reg(t);
mesh3d.Compress();
for (auto bl : mp.boundary_layers)

View File

@ -1472,7 +1472,7 @@ 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,
@ -1492,7 +1492,7 @@ py::arg("point_tolerance") = -1.)
blp.keep_surfaceindex = keep_surfaceindex;
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.")
@ -1725,7 +1725,7 @@ 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,
@ -1747,7 +1747,7 @@ py::arg("point_tolerance") = -1.)
blp.keep_surfaceindex = keep_surfaceindex;
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")=nullopt,