Some fixes for boundary layers when adjacent faces are mapped to different new materials

This commit is contained in:
Matthias Hochsteger 2025-02-25 18:04:18 +01:00
parent 8b0b9e507f
commit 5ab7a4995c
6 changed files with 139 additions and 4 deletions

View File

@ -3,6 +3,7 @@
#include <regex>
#include <set>
#include <variant>
#include "debugging.hpp"
#include "global.hpp"
@ -1095,6 +1096,8 @@ void BoundaryLayerTool ::SetDomInOutSides()
if (done.Test(index))
continue;
done.SetBit(index);
if (index < nfd_old && moved_surfaces.Test(index))
continue;
auto& fd = mesh.GetFaceDescriptor(index);
if (fd.DomainIn() != -1)
continue;
@ -1110,7 +1113,7 @@ void BoundaryLayerTool ::SetDomInOutSides()
if (e2)
dom[1] = mesh.VolumeElement(e2).GetIndex();
auto& fd_old = mesh.GetFaceDescriptor(inv_si_map[index]);
const auto& fd_old = mesh.GetFaceDescriptor(inv_si_map[index]);
int dom_old[2] = {fd_old.DomainIn(), fd_old.DomainOut()};
for (auto i : Range(2))
@ -1125,7 +1128,9 @@ void BoundaryLayerTool ::SetDomInOutSides()
}
// Check if the old domain adjacent to this face gets a new boundary layer domain, if so, use that number
int dom_new = domains.Test(dom_old[i]) ? new_mat_nrs[dom_old[i]] : dom_old[i];
int dom_new = dom_old[i];
if (domains.Test(dom_old[i]) && new_mat_nrs[dom_old[i]] > 0)
dom_new = new_mat_nrs[dom_old[i]];
// This case is tested by test_boundarylayer.py::test_pyramids[False] -> look at the generated mesh to understand the text below :)
// Special case check here: when growing "outside" the new face could have the same domain on both sides (before adding blayer elements).
@ -1249,7 +1254,10 @@ void BoundaryLayerTool ::ProcessParameters()
if (string* mat = get_if<string>(&*params.new_material); mat)
par_new_mat = {{".*", *mat}};
else
par_new_mat = *get_if<map<string, string>>(&*params.new_material);
{
par_new_mat = *get_if<map<string, string>>(&*params.new_material);
have_material_map = true;
}
}
if (params.project_boundaries.has_value())
@ -1416,6 +1424,12 @@ void BoundaryLayerTool ::Perform()
mesh.SetNextMajorTimeStamp();
mesh.UpdateTopology();
SetDomInOutSides();
if (have_material_map)
{
AddFacesBetweenDomains(mesh);
mesh.SplitFacesByAdjacentDomains();
}
}
void GenerateBoundaryLayer (Mesh& mesh, const BoundaryLayerParameters& blp)

View File

@ -67,6 +67,7 @@ public:
Array<int> par_surfid;
bool insert_only_volume_elements;
map<string, string> par_new_mat;
bool have_material_map = false;
Array<size_t> par_project_boundaries;
bool have_single_segments;

View File

@ -151,6 +151,11 @@ void BoundaryLayerTool ::InterpolateGrowthVectors()
no_angles = false;
}
}
if (no_angles && faces.Size() == 2 && have_material_map)
if (par_new_mat[mesh.GetBCName(mesh[faces[0]].GetIndex() - 1)] != par_new_mat[mesh.GetBCName(mesh[faces[1]].GetIndex() - 1)])
no_angles = false;
if (no_angles)
{
for (auto* p_seg : edgenr2seg[edgenr])

View File

@ -1,7 +1,9 @@
#include <algorithm>
#include <mystdlib.h>
#include <atomic>
#include <regex>
#include <set>
#include "core/array.hpp"
#include "meshing.hpp"
#include "../general/gzstream.h"
@ -7292,15 +7294,23 @@ namespace netgen
int eli0, eli1;
GetTopology().GetSurface2VolumeElement(sei+1, eli0, eli1);
// auto [ei0,ei1] = GetTopology().GetSurface2VolumeElement(sei); // the way to go
if(eli0 == 0)
continue;
auto & sel = (*this)[sei];
int face = sel.GetIndex();
int domin = VolumeElement(eli0).GetIndex();
int domout = eli1 ? VolumeElement(eli1).GetIndex() : 0;
if(domin < domout)
swap(domin, domout);
auto key = std::make_tuple(face, domin, domout);
if(face_doms_2_new_face.find(key) == face_doms_2_new_face.end())
{
{
auto & fd = FaceDescriptors()[face-1];
if(domout == 0 && min(fd.DomainIn(), fd.DomainOut()) > 0)
continue;
}
if(!first_visit[face-1]) {
nfaces++;
FaceDescriptor new_fd = FaceDescriptors()[face-1];
@ -7817,4 +7827,108 @@ namespace netgen
return nm_;
}
void AddFacesBetweenDomains(Mesh & mesh)
{
static Timer timer("AddFacesBetweenDomains"); RegionTimer rt(timer);
auto & topo = mesh.GetTopology();
auto p2el = mesh.CreatePoint2ElementTable();
Array<size_t> els_per_domain(mesh.GetNDomains()+1);
els_per_domain = 0;
for(const auto & el : mesh.VolumeElements())
els_per_domain[el.GetIndex()]++;
std::map<tuple<int,int>, int> doms_2_new_face;
for(const auto & [facei, fd]: Enumerate(mesh.FaceDescriptors()))
{
auto dom0 = fd.DomainIn();
auto dom1 = fd.DomainOut();
if(dom0 > dom1)
swap(dom0, dom1);
doms_2_new_face[{dom0, dom1}] = facei+1;
}
for(auto dom : Range(1, 1+mesh.GetNDomains()))
{
if(els_per_domain[dom] == 0)
continue;
mesh.UpdateTopology();
mesh.FindOpenElements(dom);
for(const auto & openel : mesh.OpenElements())
{
std::set<ElementIndex> has_p1, has_p2, has_p3;
for (auto ei: topo.GetVertexElements(openel[0]))
has_p1.insert(ei);
for (auto ei: topo.GetVertexElements(openel[1]))
has_p2.insert(ei);
for (auto ei: topo.GetVertexElements(openel[2]))
has_p3.insert(ei);
std::set<ElementIndex> has_p12, has_all;
set_intersection(has_p1.begin(), has_p1.end(),
has_p2.begin(), has_p2.end(),
inserter(has_p12, has_p12.begin()));
set_intersection(has_p12.begin(), has_p12.end(),
has_p3.begin(), has_p3.end(),
inserter(has_all, has_all.begin()));
ArrayMem<ElementIndex, 5> els;
for(auto ei : has_all)
els.Append(ei);
if(els.Size() == 2 && mesh[els[0]].GetIndex() != mesh[els[1]].GetIndex())
{
auto dom0 = mesh[els[0]].GetIndex();
auto dom1 = mesh[els[1]].GetIndex();
ElementIndex ei0 = els[0];
if(dom0 > dom1)
{
Swap(dom0, dom1);
ei0 = els[1];
}
if(dom1 == dom)
continue;
if(doms_2_new_face.count({dom0, dom1}) == 0)
{
auto fd = FaceDescriptor(-1, dom0, dom1, -1);
auto new_si = mesh.GetNFD()+1;
fd.SetBCProperty(new_si);
auto new_face = mesh.AddFaceDescriptor(fd);
mesh.SetBCName(new_si - 1, "default");
doms_2_new_face[{dom0, dom1}] = new_face;
}
for(auto face : topo.GetFaces(ei0)) {
auto verts = topo.GetFaceVertices(face);
if(verts.Contains(openel[0]) && verts.Contains(openel[1]) && verts.Contains(openel[2])) {
Element2d sel(static_cast<int>(verts.Size()));
sel.SetIndex(doms_2_new_face[{dom0, dom1}]);
for(auto j : Range(verts.Size()))
sel[j] = verts[j];
auto normal = Cross(mesh[sel[1]]-mesh[sel[0]], mesh[sel[2]]-mesh[sel[0]]);
Vec<3> surf_center = Vec<3>(Center(mesh[sel[0]] , mesh[sel[1]] , mesh[sel[2]]));
Vec<3> center(0., 0., 0.);
for(auto pi : mesh[ei0].PNums())
center += Vec<3>(mesh[pi]);
center *= 1.0/mesh[ei0].GetNP();
if((normal * (center - surf_center)) < 0)
sel.Invert();
mesh.AddSurfaceElement(sel);
break;
}
}
}
}
}
}
}

View File

@ -1035,6 +1035,7 @@ namespace netgen
return FlatArray<T_FACE>(GetNFaces ( (*mesh)[elnr].GetType()), &faces[elnr][0]);
}
DLL_HEADER void AddFacesBetweenDomains(Mesh & mesh);
}

View File

@ -105,7 +105,7 @@ namespace netgen
for( auto dom : {dom_in, dom_out} )
{
if(dom==0)
if(dom<=0)
continue;
auto & sels = ret[dom-1].mesh->SurfaceElements();