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

boundarylayers - inner corners at end of layer now possible too

Browse Source
This commit is contained in:
Christopher Lackner 2021-03-02 11:47:40 +01:00
parent 4c134567fd
commit f11cb4fcfb
2 changed files with 91 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

54
libsrc/meshing/boundarylayer.cpp
View File

@ -515,28 +515,52 @@ namespace netgen
{
if(el.GetType() == TET)
{
if(moved.Size() == 2)
if(moved.Size() == 3) // inner corner
{
if(fixed.Size() == 2)
throw Exception("This should not be possible!");
PointIndex p1 = moved[0];
PointIndex p2 = moved[1];
PointIndex p3 = moved[2];
auto v1 = mesh[p1];
auto n = Cross(mesh[p2]-v1, mesh[p3]-v1);
auto d = mesh[mapto[p1][0]] - v1;
if(n*d > 0)
Swap(p2,p3);
PointIndex p4 = p1;
PointIndex p5 = p2;
PointIndex p6 = p3;
for(auto i : Range(blp.heights))
{
PointIndex p3 = mapto[moved[1]][i];
PointIndex p4 = mapto[moved[0]][i];
Element nel(PYRAMID);
nel[0] = p1;
nel[1] = p2;
nel[2] = p3;
nel[3] = p4;
nel[4] = el[0] + el[1] + el[2] + el[3] - fixed[0] - moved[0] - moved[1];
if(Cross(mesh[p2]-mesh[p1], mesh[p4]-mesh[p1]) * (mesh[nel[4]]-mesh[nel[1]]) > 0)
Swap(nel[1], nel[3]);
Element nel(PRISM);
nel[0] = p4; nel[1] = p5; nel[2] = p6;
p4 = mapto[p1][i]; p5 = mapto[p2][i]; p6 = mapto[p3][i];
nel[3] = p4; nel[4] = p5; nel[5] = p6;
nel.SetIndex(el.GetIndex());
mesh.AddVolumeElement(nel);
p1 = p4;
p2 = p3;
}
}
if(moved.Size() == 2)
{
if(fixed.Size() == 1)
{
PointIndex p1 = moved[0];
PointIndex p2 = moved[1];
for(auto i : Range(blp.heights))
{
PointIndex p3 = mapto[moved[1]][i];
PointIndex p4 = mapto[moved[0]][i];
Element nel(PYRAMID);
nel[0] = p1;
nel[1] = p2;
nel[2] = p3;
nel[3] = p4;
nel[4] = el[0] + el[1] + el[2] + el[3] - fixed[0] - moved[0] - moved[1];
if(Cross(mesh[p2]-mesh[p1], mesh[p4]-mesh[p1]) * (mesh[nel[4]]-mesh[nel[1]]) > 0)
Swap(nel[1], nel[3]);
nel.SetIndex(el.GetIndex());
mesh.AddVolumeElement(nel);
p1 = p4;
p2 = p3;
}
}
}
if(moved.Size() == 1 && fixed.Size() == 1)

53
tests/pytest/test_boundarylayer.py
View File

@ -66,7 +66,7 @@ def test_boundarylayer2(outside, version, capfd):
@pytest.mark.parametrize("outside", [True, False])
def test_wrong_orientation(outside):
def test_wrong_orientation(outside, capfd):
geo = CSGeometry()
brick = OrthoBrick((-1,0,0),(1,1,1)) - Plane((0,0,0), (1,0,0))
geo.Add(brick.mat("air"))
@ -107,3 +107,54 @@ def test_pyramids(outside):
assert ngs.Integrate(1, mesh.Materials("plate")) == pytest.approx(0.032 if outside else 0.0304)
assert ngs.Integrate(1, mesh.Materials("layer")) == pytest.approx(0.0016)
assert ngs.Integrate(1, mesh.Materials("air")) == pytest.approx(0.9664 if outside else 0.968)
@pytest.mark.parametrize("outside", [True, False])
def test_with_inner_corner(outside, capfd):
geo = CSGeometry()
core_thickness = 0.1
limb_distance = 0.5
core_height = 0.5
coil_r1 = 0.08
coil_r2 = 0.16
coil_h = 0.2
domain_size = 1.2
domain_size_y = 0.7
def CreateCoil(x):
outer = Cylinder((x, 0, -1), (x, 0, 1), coil_r2)
inner = Cylinder((x, 0, -1), (x, 0, 1), coil_r1)
top = Plane((0,0,coil_h/2), (0,0,1))
bot = Plane((0,0,-coil_h/2), (0,0,-1))
return ((outer - inner) * top * bot, (outer, inner, top, bot))
core_front = Plane((0,-core_thickness/2, 0), (0,-1,0)).bc("core_front")
core_back = Plane((0,core_thickness/2, 0), (0,1,0)).bc("core_front")
core_limb1 = OrthoBrick((-limb_distance/2-core_thickness/2, -core_thickness, -core_height/2),(-limb_distance/2+core_thickness/2, core_thickness, core_height/2))
core_limb2 = OrthoBrick((limb_distance/2-core_thickness/2, -core_thickness, -core_height/2),(limb_distance/2+core_thickness/2, core_thickness, core_height/2))
core_top = OrthoBrick((-limb_distance/2-core_thickness/2, -core_thickness, core_height/2-core_thickness/2),(limb_distance/2+core_thickness/2, core_thickness, core_height/2+core_thickness/2))
core_bot = OrthoBrick((-limb_distance/2-core_thickness/2, -core_thickness, -core_height/2-core_thickness/2),(limb_distance/2+core_thickness/2, core_thickness, -core_height/2+core_thickness/2))
core = (core_limb1 + core_limb2 + core_top + core_bot).bc("core_rest")
core = core * core_front * core_back
core.maxh(core_thickness * 0.4)
coil1, (outer1, inner1, top1, bot1) = CreateCoil(-limb_distance/2)
coil1.mat("coil_1")
coil2, (outer2, inner2, top2, bot2) = CreateCoil(limb_distance/2)
coil2.mat("coil_2")
oil = OrthoBrick((-domain_size/2, -domain_size_y/2, -domain_size/2), (domain_size/2, domain_size_y/2, domain_size/2)).bc("tank") - core # - coil1 - coil2
geo.Add(core.mat("core"), col=(0.4,0.4,0.4))
geo.Add(coil1, col=(0.72, 0.45, 0.2))
geo.Add(coil2, col=(0.72, 0.45, 0.2))
geo.Add(oil.mat("oil"), transparent=True)
mesh = geo.GenerateMesh()
mesh.BoundaryLayer("core_front", [0.001, 0.002], "core", "core", outside=outside)
ngs = pytest.importorskip("ngsolve")
mesh = ngs.Mesh(mesh)
capture = capfd.readouterr()
assert not "elements are not matching" in capture.out
assert ngs.Integrate(1, mesh.Materials("core")) == pytest.approx(0.0212 if outside else 0.02)
assert ngs.Integrate(1, mesh.Materials("oil")) == pytest.approx(0.9868 if outside else 0.988)