From f11cb4fcfb9314d842d52869fe8075b2dd3732d4 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Tue, 2 Mar 2021 11:47:40 +0100 Subject: [PATCH] boundarylayers - inner corners at end of layer now possible too --- libsrc/meshing/boundarylayer.cpp | 54 +++++++++++++++++++++--------- tests/pytest/test_boundarylayer.py | 53 ++++++++++++++++++++++++++++- 2 files changed, 91 insertions(+), 16 deletions(-) diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index d0bfe335..4539887c 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -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) diff --git a/tests/pytest/test_boundarylayer.py b/tests/pytest/test_boundarylayer.py index f9310955..798cb1cc 100644 --- a/tests/pytest/test_boundarylayer.py +++ b/tests/pytest/test_boundarylayer.py @@ -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)