fixes for boundarylayer edge tangent computation and some more

This commit is contained in:
Christopher Lackner 2022-06-07 14:51:41 +02:00 committed by Matthias Hochsteger
parent 6b12050fd1
commit a3408b537a
4 changed files with 81 additions and 23 deletions

View File

@ -194,7 +194,8 @@ namespace netgen
ArrayMem<Point<3>, 4> BoundaryLayerTool :: GetMappedFace( SurfaceElementIndex sei, int face ) ArrayMem<Point<3>, 4> BoundaryLayerTool :: GetMappedFace( SurfaceElementIndex sei, int face )
{ {
if(face == -1) return GetMappedFace(sei); if(face == -1) return GetFace(sei);
if(face == -2) return GetMappedFace(sei);
const auto & sel = mesh[sei]; const auto & sel = mesh[sei];
auto np = sel.GetNP(); auto np = sel.GetNP();
auto pi0 = sel[face % np]; auto pi0 = sel[face % np];
@ -210,23 +211,24 @@ namespace netgen
Vec<3> BoundaryLayerTool :: getEdgeTangent(PointIndex pi, int edgenr) Vec<3> BoundaryLayerTool :: getEdgeTangent(PointIndex pi, int edgenr)
{ {
Vec<3> tangent = 0.0; Vec<3> tangent = 0.0;
ArrayMem<PointIndex,2> pts;
for(auto segi : topo.GetVertexSegments(pi)) for(auto segi : topo.GetVertexSegments(pi))
{ {
auto & seg = mesh[segi]; auto & seg = mesh[segi];
if(seg.edgenr != edgenr) if(seg.edgenr != edgenr+1)
continue; continue;
PointIndex other = seg[0]+seg[1]-pi; PointIndex other = seg[0]+seg[1]-pi;
tangent += mesh[other] - mesh[pi]; pts.Append(other);
} }
if(pts.Size() != 2)
throw Exception("Something went wrong in getEdgeTangent!");
tangent = mesh[pts[1]] - mesh[pts[0]];
return tangent.Normalize(); return tangent.Normalize();
} }
void BoundaryLayerTool :: LimitGrowthVectorLengths() void BoundaryLayerTool :: LimitGrowthVectorLengths()
{ {
static Timer tall("BoundaryLayerTool::LimitGrowthVectorLengths"); RegionTimer rtall(tall); static Timer tall("BoundaryLayerTool::LimitGrowthVectorLengths"); RegionTimer rtall(tall);
height = 0.0;
for (auto h : params.heights)
height += h;
limits.SetSize(np); limits.SetSize(np);
limits = 1.0; limits = 1.0;
@ -320,7 +322,7 @@ namespace netgen
const auto & sel = mesh[sei]; const auto & sel = mesh[sei];
if(sel.PNums().Contains(pi)) if(sel.PNums().Contains(pi))
return false; return false;
auto face = GetFace(sei); auto face = GetMappedFace(sei, -2);
double lam_ = 999; double lam_ = 999;
bool is_bl_sel = params.surfid.Contains(sel.GetIndex()); bool is_bl_sel = params.surfid.Contains(sel.GetIndex());
@ -464,6 +466,8 @@ namespace netgen
auto np = mesh.GetNP(); auto np = mesh.GetNP();
BitArray is_point_on_bl_surface(np+1); BitArray is_point_on_bl_surface(np+1);
is_point_on_bl_surface.Clear(); is_point_on_bl_surface.Clear();
BitArray is_point_on_other_surface(np+1);
is_point_on_other_surface.Clear();
Array<Vec<3>, PointIndex> normals(np); Array<Vec<3>, PointIndex> normals(np);
for(auto pi : Range(growthvectors)) for(auto pi : Range(growthvectors))
normals[pi] = growthvectors[pi]; normals[pi] = growthvectors[pi];
@ -474,20 +478,33 @@ namespace netgen
{ {
auto facei = mesh[sei].GetIndex(); auto facei = mesh[sei].GetIndex();
if(facei < nfd_old && !params.surfid.Contains(facei)) if(facei < nfd_old && !params.surfid.Contains(facei))
continue; {
for(auto pi : mesh[sei].PNums()) for(auto pi : mesh[sei].PNums())
if(mesh[pi].Type() == SURFACEPOINT) if(mesh[pi].Type() == SURFACEPOINT)
is_point_on_other_surface.SetBitAtomic(pi);
}
else
{
for(auto pi : mesh[sei].PNums())
if(mesh[pi].Type() == SURFACEPOINT)
is_point_on_bl_surface.SetBitAtomic(pi); is_point_on_bl_surface.SetBitAtomic(pi);
}
} }
}); });
Array<PointIndex> points; Array<PointIndex> points;
for(PointIndex pi : mesh.Points().Range()) for(PointIndex pi : mesh.Points().Range())
{
if(is_point_on_bl_surface[pi]) if(is_point_on_bl_surface[pi])
{ {
points.Append(pi); points.Append(pi);
growthvectors[pi] = 0.0; growthvectors[pi] = 0.0;
} }
if(is_point_on_other_surface[pi])
{
points.Append(pi);
}
}
// smooth tangential part of growth vectors from edges to surface elements // smooth tangential part of growth vectors from edges to surface elements
RegionTimer rtsmooth(tsmooth); RegionTimer rtsmooth(tsmooth);
@ -511,7 +528,10 @@ namespace netgen
auto gw_other = growthvectors[pi1]; auto gw_other = growthvectors[pi1];
auto normal_other = getNormal(mesh[sei]); auto normal_other = getNormal(mesh[sei]);
auto tangent_part = gw_other - (gw_other*normal_other)*normal_other; auto tangent_part = gw_other - (gw_other*normal_other)*normal_other;
new_gw += tangent_part; if(is_point_on_bl_surface[pi])
new_gw += tangent_part;
else
new_gw += gw_other;
} }
} }
@ -533,6 +553,10 @@ namespace netgen
//for(auto & seg : mesh.LineSegments()) //for(auto & seg : mesh.LineSegments())
//seg.edgenr = seg.epgeominfo[1].edgenr; //seg.edgenr = seg.epgeominfo[1].edgenr;
height = 0.0;
for (auto h : params.heights)
height += h;
max_edge_nr = -1; max_edge_nr = -1;
for(const auto& seg : mesh.LineSegments()) for(const auto& seg : mesh.LineSegments())
if(seg.edgenr > max_edge_nr) if(seg.edgenr > max_edge_nr)
@ -777,7 +801,7 @@ namespace netgen
// interpolate tangential component of growth vector along edge // interpolate tangential component of growth vector along edge
for(auto edgenr : Range(max_edge_nr)) for(auto edgenr : Range(max_edge_nr))
{ {
if(!is_edge_moved[edgenr+1]) continue; // if(!is_edge_moved[edgenr+1]) continue;
// build sorted list of edge // build sorted list of edge
Array<PointIndex> points; Array<PointIndex> points;
@ -785,6 +809,9 @@ namespace netgen
double edge_len = 0.; double edge_len = 0.;
auto is_end_point = [&] (PointIndex pi) auto is_end_point = [&] (PointIndex pi)
{ {
// if(mesh[pi].Type() == FIXEDPOINT)
// return true;
// return false;
auto segs = topo.GetVertexSegments(pi); auto segs = topo.GetVertexSegments(pi);
auto first_edgenr = mesh[segs[0]].edgenr; auto first_edgenr = mesh[segs[0]].edgenr;
for(auto segi : segs) for(auto segi : segs)
@ -792,17 +819,31 @@ namespace netgen
return true; return true;
return false; return false;
}; };
bool any_grows = false;
for(const auto& seg : segments) for(const auto& seg : segments)
{ {
if(seg.edgenr-1 == edgenr && is_end_point(seg[0])) if(seg.edgenr-1 == edgenr)
{ {
points.Append(seg[0]); if(growthvectors[seg[0]].Length2() != 0 ||
points.Append(seg[1]); growthvectors[seg[1]].Length2() != 0)
edge_len += (mesh[seg[1]] - mesh[seg[0]]).Length(); any_grows = true;
break; if(points.Size() == 0 && is_end_point(seg[0]))
{
points.Append(seg[0]);
points.Append(seg[1]);
edge_len += (mesh[seg[1]] - mesh[seg[0]]).Length();
}
} }
} }
if(!any_grows)
continue;
if(!points.Size())
throw Exception("Could not find startpoint for edge " + ToString(edgenr));
while(true) while(true)
{ {
bool point_found = false; bool point_found = false;
@ -831,17 +872,28 @@ namespace netgen
break; break;
if(!point_found) if(!point_found)
{ {
cout << "points = " << points << endl;
throw Exception(string("Could not find connected list of line segments for edge ") + edgenr); throw Exception(string("Could not find connected list of line segments for edge ") + edgenr);
} }
} }
if(growthvectors[points[0]].Length2() == 0 &&
growthvectors[points.Last()].Length2() == 0)
continue;
// tangential part of growth vectors // tangential part of growth vectors
auto t1 = (mesh[points[1]]-mesh[points[0]]).Normalize(); auto t1 = (mesh[points[1]]-mesh[points[0]]).Normalize();
auto gt1 = growthvectors[points[0]] * t1 * t1; auto gt1 = growthvectors[points[0]] * t1 * t1;
auto t2 = (mesh[points.Last()]-mesh[points[points.Size()-2]]).Normalize(); auto t2 = (mesh[points.Last()]-mesh[points[points.Size()-2]]).Normalize();
auto gt2 = growthvectors[points.Last()] * t2 * t2; auto gt2 = growthvectors[points.Last()] * t2 * t2;
if(!is_edge_moved[edgenr+1])
{
if(growthvectors[points[0]] * (mesh[points[1]] - mesh[points[0]]) < 0)
gt1 = 0.;
if(growthvectors[points.Last()] * (mesh[points[points.Size()-2]] - mesh[points.Last()]) < 0)
gt2 = 0.;
}
double len = 0.; double len = 0.;
for(size_t i = 1; i < points.Size()-1; i++) for(size_t i = 1; i < points.Size()-1; i++)
{ {
@ -1262,7 +1314,9 @@ namespace netgen
auto in_surface_direction = ProjectGrowthVectorsOnSurface(); auto in_surface_direction = ProjectGrowthVectorsOnSurface();
InterpolateGrowthVectors(); InterpolateGrowthVectors();
LimitGrowthVectorLengths();
if(params.limit_growth_vectors)
LimitGrowthVectorLengths();
FixVolumeElements(); FixVolumeElements();
InsertNewElements(segmap, in_surface_direction); InsertNewElements(segmap, in_surface_direction);
SetDomInOut(); SetDomInOut();

View File

@ -18,6 +18,7 @@ public:
BitArray domains; BitArray domains;
bool outside = false; // set the boundary layer on the outside bool outside = false; // set the boundary layer on the outside
bool grow_edges = false; bool grow_edges = false;
bool limit_growth_vectors = true;
Array<size_t> project_boundaries; Array<size_t> project_boundaries;
}; };

View File

@ -250,6 +250,8 @@ namespace netgen
el.SetIndex(md.domain); el.SetIndex(md.domain);
mesh.AddVolumeElement(el); mesh.AddVolumeElement(el);
// TODO: Fix double hexes
return;
} }
} }
} }
@ -577,8 +579,8 @@ namespace netgen
if (md[i].mesh->CheckOverlappingBoundary()) if (md[i].mesh->CheckOverlappingBoundary())
throw NgException ("Stop meshing since boundary mesh is overlapping"); throw NgException ("Stop meshing since boundary mesh is overlapping");
if(md[i].mesh->GetGeometry()->GetGeomType() == Mesh::GEOM_OCC) // if(md[i].mesh->GetGeometry()->GetGeomType() == Mesh::GEOM_OCC)
FillCloseSurface( md[i] ); // FillCloseSurface( md[i] );
CloseOpenQuads( md[i] ); CloseOpenQuads( md[i] );
MeshDomain(md[i]); MeshDomain(md[i]);
}, md.Size()); }, md.Size());

View File

@ -1196,7 +1196,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
string material, string material,
variant<string, int> domain, bool outside, variant<string, int> domain, bool outside,
optional<string> project_boundaries, optional<string> project_boundaries,
bool grow_edges) bool grow_edges, bool limit_growth_vectors)
{ {
BoundaryLayerParameters blp; BoundaryLayerParameters blp;
if(int* bc = get_if<int>(&boundary); bc) if(int* bc = get_if<int>(&boundary); bc)
@ -1265,12 +1265,13 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
blp.outside = outside; blp.outside = outside;
blp.grow_edges = grow_edges; blp.grow_edges = grow_edges;
blp.limit_growth_vectors = limit_growth_vectors;
GenerateBoundaryLayer (self, blp); GenerateBoundaryLayer (self, blp);
self.UpdateTopology(); self.UpdateTopology();
}, py::arg("boundary"), py::arg("thickness"), py::arg("material"), }, py::arg("boundary"), py::arg("thickness"), py::arg("material"),
py::arg("domains") = ".*", py::arg("outside") = false, py::arg("domains") = ".*", py::arg("outside") = false,
py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true, py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true, py::arg("limit_growth_vectors") = true,
R"delimiter( R"delimiter(
Add boundary layer to mesh. Add boundary layer to mesh.