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

Cleanup, fix growthvector interpolation

Browse Source
This commit is contained in:
Matthias Hochsteger 2024-04-10 17:34:39 +02:00
parent cb2a8ddbb1
commit 7cd448c958
2 changed files with 52 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

90
libsrc/meshing/boundarylayer.cpp
View File

@ -564,27 +564,31 @@ void MergeAndAddSegments(Mesh& mesh, FlatArray<Segment> segments,
for (const auto& seg : new_segments) addSegment(seg); for (const auto& seg : new_segments) addSegment(seg);
} }
// TODO: Hack, move this to the header or restructure the whole growth_vectors storage
static std::map<PointIndex, Vec<3>> non_bl_growth_vectors;
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");
auto np_old = this->np; auto np_old = this->np;
auto np = mesh.GetNP(); auto np = mesh.GetNP();
BitArray is_point_on_bl_surface(np + 1);
is_point_on_bl_surface.Clear(); non_bl_growth_vectors.clear();
BitArray is_point_on_other_surface(np + 1);
is_point_on_other_surface.Clear();
auto getGW = [&](PointIndex pi) -> Vec<3> { auto getGW = [&](PointIndex pi) -> Vec<3> {
if (pi - PointIndex::BASE < np_old && mapto[pi].Size()==0) return {.0}; if (growth_vector_map.count(pi) == 0) {
if (growth_vector_map.count(pi) == 0) non_bl_growth_vectors[pi] = .0;
growth_vector_map[pi] = {&growthvectors[pi], total_height}; 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) {
growth_vector_map[pi] = {&growthvectors[pi], total_height}; 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;
}; };
@ -600,14 +604,31 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
}; };
std::set<PointIndex> points_set; std::set<PointIndex> points_set;
ParallelForRange(Range((size_t)nse, mesh.SurfaceElements().Size()), [&](auto myrange) { ParallelForRange(mesh.SurfaceElements().Range(), [&](auto myrange) {
for (SurfaceElementIndex sei : myrange) { for (SurfaceElementIndex sei : myrange) {
for (auto pi : mesh[sei].PNums()) for (auto pi : mesh[sei].PNums()) {
if (auto& p = mesh[mapfrom[pi]]; p.Type() == SURFACEPOINT) 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))
points_set.insert(pi); points_set.insert(pi);
}
} }
}); });
Array<bool> has_moved_points(max_edge_nr + 1);
has_moved_points = false;
std::set<PointIndex> moved_edge_points;
for (auto seg : segments) {
if (hasMoved(seg[0]) != hasMoved(seg[1]))
has_moved_points[seg.edgenr] = true;
}
for (auto seg : segments)
if (has_moved_points[seg.edgenr])
for (auto pi : seg.PNums())
if (mesh[pi].Type() == EDGEPOINT) points_set.insert(pi);
Array<PointIndex> points; Array<PointIndex> points;
for (auto pi : points_set) for (auto pi : points_set)
points.Append(pi); points.Append(pi);
@ -622,25 +643,31 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
for (auto pi : points) { for (auto pi : points) {
auto sels = p2sel[pi]; auto sels = p2sel[pi];
auto & correction = corrections[pi]; auto & correction = corrections[pi];
correction = 0.0;
std::set<PointIndex> suround; std::set<PointIndex> suround;
suround.insert(pi); suround.insert(pi);
double total_weight = 0;
// average only tangent component on new bl points, average whole growth vector otherwise
bool do_average_tangent = mapfrom[pi].IsValid();
correction = 0.0;
for (auto sei : sels) { for (auto sei : sels) {
const auto& sel = mesh[sei]; const auto& sel = mesh[sei];
for (auto pi1 : sel.PNums()) for (auto pi1 : sel.PNums()) {
if (suround.count(pi1) == 0) { if (suround.count(pi1)) continue;
suround.insert(pi1); suround.insert(pi1);
auto gw_other = getGW(pi1)+corrections[pi1]; auto gw_other = getGW(pi1)+corrections[pi1];
if(do_average_tangent) {
auto normal_other = getNormal(mesh[sei]); auto normal_other = getNormal(mesh[sei]);
auto tangent_part = auto tangent_part = gw_other - (gw_other * normal_other) * normal_other;
gw_other - (gw_other * normal_other) * normal_other; correction += tangent_part;
double weight = 1.0;
correction += weight * tangent_part;
total_weight += weight;
} }
else {
correction += gw_other;
}
}
} }
correction *= 1.0 / total_weight; correction *= 1.0 / suround.size();
if(!do_average_tangent)
correction -= getGW(pi);
} }
} }
for(auto pi: points) for(auto pi: points)
@ -1050,9 +1077,6 @@ void BoundaryLayerTool ::InterpolateGrowthVectors() {
// gt2 = 0.; // gt2 = 0.;
// } // }
// cout << "edgenr " << edgenr << endl;
// cout << "points " << endl << points << endl;
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]; auto pi = points[i];
@ -1060,13 +1084,6 @@ void BoundaryLayerTool ::InterpolateGrowthVectors() {
auto t = getEdgeTangent(pi, edgenr); auto t = getEdgeTangent(pi, edgenr);
auto lam = len / edge_len; auto lam = len / edge_len;
auto interpol = (1 - lam) * (gt1 * t) * t + lam * (gt2 * t) * t; auto interpol = (1 - lam) * (gt1 * t) * t + lam * (gt2 * t) * t;
// if(pi==89) {
// cout << "points " << points << endl;
// cout << "INTERPOL" << len << ',' << t << ',' << lam << ',' <<
// interpol << endl; cout << gt1 << endl; cout << gt2 << endl; cout <<
// getGW(pi) << endl;
// }
addGW(pi, interpol); addGW(pi, interpol);
} }
} }
@ -1214,7 +1231,6 @@ void BoundaryLayerTool ::InsertNewElements(
s0.edgenr = segj.edgenr; s0.edgenr = segj.edgenr;
s0.si = segj.si; s0.si = segj.si;
new_segments.Append(s0); new_segments.Append(s0);
// cout << __LINE__ <<"\t" << s0 << endl;
for (auto i : Range(params.heights)) { for (auto i : Range(params.heights)) {
Element2d sel(QUAD); Element2d sel(QUAD);
@ -1240,7 +1256,6 @@ void BoundaryLayerTool ::InsertNewElements(
s1.edgenr = getEdgeNr(segj.edgenr); s1.edgenr = getEdgeNr(segj.edgenr);
s1.si = segj.si; s1.si = segj.si;
// new_segments.Append(s1); // new_segments.Append(s1);
// cout << __LINE__ <<"\t" << s1 << endl;
Segment s2; Segment s2;
s2[0] = p4; s2[0] = p4;
s2[1] = p1; s2[1] = p1;
@ -1249,7 +1264,6 @@ void BoundaryLayerTool ::InsertNewElements(
s2.edgenr = getEdgeNr(segj.edgenr); s2.edgenr = getEdgeNr(segj.edgenr);
s2.si = segj.si; s2.si = segj.si;
// new_segments.Append(s2); // new_segments.Append(s2);
// cout << __LINE__ <<"\t" << s2 << endl;
p1 = p4; p1 = p4;
p2 = p3; p2 = p3;
} }
@ -1261,7 +1275,6 @@ void BoundaryLayerTool ::InsertNewElements(
s3.edgenr = getEdgeNr(segj.edgenr); s3.edgenr = getEdgeNr(segj.edgenr);
s3.si = segj.si; s3.si = segj.si;
new_segments.Append(s3); new_segments.Append(s3);
// cout << __LINE__ << "\t" << s3 << endl;
} }
} }
} }
@ -1661,6 +1674,7 @@ void BoundaryLayerTool ::Perform() {
InsertNewElements(segmap, in_surface_direction); InsertNewElements(segmap, in_surface_direction);
mapfrom.SetSize(mesh.GetNP()); mapfrom.SetSize(mesh.GetNP());
mapfrom = PointIndex::INVALID;
for (auto pi : mapto.Range()) for (auto pi : mapto.Range())
for (auto pi_to : mapto[pi]) mapfrom[pi_to] = pi; for (auto pi_to : mapto[pi]) mapfrom[pi_to] = pi;

6
libsrc/meshing/meshfunc.cpp
View File

@ -50,7 +50,6 @@ namespace netgen
ret[0].mp = mp; ret[0].mp = mp;
return ret; return ret;
} }
cout << "divide mesh" << endl;
ret.SetSize(num_domains); ret.SetSize(num_domains);
Array<Array<PointIndex, PointIndex>> ipmap; Array<Array<PointIndex, PointIndex>> ipmap;
@ -72,7 +71,6 @@ namespace netgen
m.SetLocalH(mesh.GetLocalH()); m.SetLocalH(mesh.GetLocalH());
cout << "set imap size " << i << " to " << num_points << endl;
ipmap[i].SetSize(num_points); ipmap[i].SetSize(num_points);
ipmap[i] = PointIndex::INVALID; ipmap[i] = PointIndex::INVALID;
m.SetDimension( mesh.GetDimension() ); m.SetDimension( mesh.GetDimension() );
@ -157,8 +155,6 @@ namespace netgen
auto & imap = ipmap[i]; auto & imap = ipmap[i];
auto nmax = identifications.GetMaxNr (); auto nmax = identifications.GetMaxNr ();
auto & m_ident = m.GetIdentifications(); auto & m_ident = m.GetIdentifications();
cout << "imap " << imap << endl;
cout << imap.Size() << endl;
for (auto & sel : m.SurfaceElements()) for (auto & sel : m.SurfaceElements())
for(auto & pi : sel.PNums()) for(auto & pi : sel.PNums())
@ -175,7 +171,6 @@ namespace netgen
for(auto pair : pairs) for(auto pair : pairs)
{ {
// cout << "get pair " << pair[0] << ',' << pair[1] << endl;
auto pi0 = imap[pair[0]]; auto pi0 = imap[pair[0]];
auto pi1 = imap[pair[1]]; auto pi1 = imap[pair[1]];
if(!pi0.IsValid() || !pi1.IsValid()) if(!pi0.IsValid() || !pi1.IsValid())
@ -307,7 +302,6 @@ namespace netgen
rulep = hexrules; rulep = hexrules;
break; break;
case 1: case 1:
cout << "Apply prismrules " << endl;
rulep = prismrules2; rulep = prismrules2;
break; break;
case 2: // connect pyramid to triangle case 2: // connect pyramid to triangle