ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-24 21:10:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

boundarylayer fixes

Browse Source
This commit is contained in:
Christopher Lackner 2020-07-14 21:30:26 +02:00
parent fb13152004
commit ec3d7c3ec9
2 changed files with 88 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
libsrc/meshing/boundarylayer.cpp
View File

@ -373,9 +373,6 @@ namespace netgen
throw Exception("Couldn't find element on other side for " + ToString(segj[0]) + " to " + ToString(segj[1])); throw Exception("Couldn't find element on other side for " + ToString(segj[0]) + " to " + ToString(segj[1]));
const auto& commsel = mesh[pnt_commelem]; const auto& commsel = mesh[pnt_commelem];
auto surfelem_vect = GetSurfaceNormal(mesh, commsel);
if(blp.outside)
surfelem_vect *= -1;
Element2d sel(QUAD); Element2d sel(QUAD);
auto seg_p1 = segi[0]; auto seg_p1 = segi[0];
auto seg_p2 = segi[1]; auto seg_p2 = segi[1];
@ -428,38 +425,78 @@ namespace netgen
mesh.AddSegment(seg_2); mesh.AddSegment(seg_2);
} }
// in first layer insert new segments adjacent to
// new face
if(layer == 1 && !blp.surfid.Contains(segj.si))
{
Segment s3 = segj;
s3.si = map_surface_index(segj.si)-1;
s3[0] = mapto[s3[0]];
s3[1] = mapto[s3[1]];
mesh.AddSegment(s3);
}
// in last layer insert new segments // in last layer insert new segments
if(layer == blp.heights.Size()) if(layer == blp.heights.Size())
{ {
max_edge_nr++;
if(!blp.surfid.Contains(segj.si))
{
Segment s3 = segj;
s3.si = map_surface_index(segj.si)-1;
Swap(s3[0], s3[1]);
if(blp.outside)
{
s3[0] = mapto[s3[0]];
s3[1] = mapto[s3[1]];
}
else
s3.edgenr = max_edge_nr;
mesh.AddSegment(s3);
}
Segment s1 = segi; Segment s1 = segi;
Segment s2 = segj; Segment s2 = segj;
s1.edgenr = ++max_edge_nr; s1.edgenr = max_edge_nr;
s2.edgenr = max_edge_nr; s2.edgenr = max_edge_nr;
auto side_surf = domains_to_surf_index[make_tuple(s2.si, blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(s2.si).DomainOut())];
if(blp.surfid.Contains(segj.si)) if(blp.surfid.Contains(segj.si))
s2.si = map_surface_index(segj.si); s2.si = map_surface_index(segj.si);
else else
{ {
auto side_surf = domains_to_surf_index[make_tuple(s2.si, blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(s2.si).DomainOut())];
if(blp.outside) if(blp.outside)
s2.si = side_surf; {
s2.si = side_surf;
}
else else
segj.si = side_surf; mesh[sej].si = side_surf;
} }
s1.si = map_surface_index(s1.si); s1.si = map_surface_index(s1.si);
s1.surfnr1 = s1.surfnr2 = s2.surfnr1 = s2.surfnr2 = -1; s1.surfnr1 = s1.surfnr2 = s2.surfnr1 = s2.surfnr2 = -1;
mesh.AddSegment(s1); mesh.AddSegment(s1);
mesh.AddSegment(s2); mesh.AddSegment(s2);
} }
segmap.SetSize(mesh.LineSegments().Size());
for(auto sei2 : segmap[sei])
{
auto& s = mesh[sei2];
if(blp.outside && layer == blp.heights.Size())
{
if(blp.surfid.Contains(s.si))
s.si = map_surface_index(s.si);
s.edgenr = max_edge_nr;
}
else
{
s[0] = mapto[s[0]];
s[1] = mapto[s[1]];
}
}
for(auto sej2 : segmap[sej])
{
auto& s = mesh[sej2];
if(blp.outside && layer == blp.heights.Size())
{
if(blp.surfid.Contains(s.si))
s.si = map_surface_index(s.si);
s.edgenr = max_edge_nr;
}
else
{
s[0] = mapto[s[0]];
s[1] = mapto[s[1]];
}
}
// do not use segi (not even with reference, since // do not use segi (not even with reference, since
// mesh.AddSegment will resize segment array and // mesh.AddSegment will resize segment array and
// invalidate reference), this is why we copy it!!! // invalidate reference), this is why we copy it!!!
@ -477,6 +514,8 @@ namespace netgen
// if necessary map them // if necessary map them
for(SegmentIndex sej = 0; sej<nseg; sej++) for(SegmentIndex sej = 0; sej<nseg; sej++)
{ {
if(segsel.Test(sej))
{
if(mesh[sej][0] == mesh[sei][1] && if(mesh[sej][0] == mesh[sei][1] &&
mesh[sej][1] == mesh[sei][0]) mesh[sej][1] == mesh[sei][0])
{ {
@ -506,8 +545,7 @@ namespace netgen
mesh[sej][1] = mapto[mesh[sej][1]]; mesh[sej][1] = mapto[mesh[sej][1]];
} }
} }
else }
continue;
} }
} }
} }

38
tests/pytest/test_boundarylayer.py
View File

@ -2,11 +2,10 @@
import pytest import pytest
from netgen.csg import * from netgen.csg import *
def GetNSurfaceElements(mesh, boundary): def GetNSurfaceElements(mesh, boundaries):
nse_in_layer = 0 nse_in_layer = 0
for el in mesh.Elements2D(): for el in mesh.Elements2D():
print(el.index) if mesh.GetBCName(el.index-1) in boundaries:
if mesh.GetBCName(el.index-1) == boundary:
nse_in_layer += 1 nse_in_layer += 1
return nse_in_layer return nse_in_layer
@ -14,18 +13,43 @@ def GetNSurfaceElements(mesh, boundary):
def test_boundarylayer(outside, capfd): def test_boundarylayer(outside, capfd):
mesh = unit_cube.GenerateMesh(maxh=0.3) mesh = unit_cube.GenerateMesh(maxh=0.3)
ne_before = mesh.ne ne_before = mesh.ne
layer_surfacenames = ["right", "top"] layer_surfacenames = ["right", "top", "left", "back", "bottom"]
mesh.BoundaryLayer("|".join(layer_surfacenames), [0.01, 0.02], "layer", outside=outside, grow_edges=True) mesh.BoundaryLayer("|".join(layer_surfacenames), [0.01, 0.02], "layer", outside=outside, grow_edges=True)
should_ne = ne_before + 2 * sum([GetNSurfaceElements(mesh, surf) for surf in layer_surfacenames]) should_ne = ne_before + 2 * GetNSurfaceElements(mesh, layer_surfacenames)
assert mesh.ne == should_ne assert mesh.ne == should_ne
capture = capfd.readouterr() capture = capfd.readouterr()
assert not "elements are not matching" in capture.out assert not "elements are not matching" in capture.out
for side in ["front"]: for side in ["front"]:
mesh.BoundaryLayer(side, [0.001], "layer", outside=outside, grow_edges=True) mesh.BoundaryLayer(side, [0.001, 0.002], "layer", outside=outside, grow_edges=True)
should_ne += GetNSurfaceElements(mesh, side) should_ne += 2 * GetNSurfaceElements(mesh, [side])
assert mesh.ne == should_ne assert mesh.ne == should_ne
capture = capfd.readouterr() capture = capfd.readouterr()
assert not "elements are not matching" in capture.out assert not "elements are not matching" in capture.out
@pytest.mark.parametrize("outside", [True, False])
@pytest.mark.parametrize("version", [1]) # version 2 not working yet
def test_boundarylayer2(outside, version, capfd):
geo = CSGeometry()
top = Plane(Pnt(0,0,0.5), Vec(0,0,1))
bot = Plane(Pnt(0,0,0.1), Vec(0,0,-1))
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"))
else:
geo.Add((outer-top).mat("rest"))
geo.Add((outer-bot).mat("rest"))
geo.Add((outer*top*bot-bigpart).mat("rest"))
geo.CloseSurfaces(top, bot, [])
mesh = geo.GenerateMesh()
should_ne = mesh.ne + 2 * GetNSurfaceElements(mesh, ["default"])
layersize = 0.05
mesh.BoundaryLayer("default", [0.5 * layersize, layersize], "layer", domains="part", outside=outside, grow_edges=True)
assert mesh.ne == should_ne
assert not "elements are not matching" in capfd.readouterr().out