ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-11-11 16:49:16 +05:00
Code Issues Packages Projects Releases Wiki Activity

Boundarylayer grows pyramids if created on interior bnd

Browse Source
This commit is contained in:
Christopher Lackner 2020-11-23 23:48:49 +01:00
parent e7b9baa93b
commit 7e78056ade
2 changed files with 116 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

107
libsrc/meshing/boundarylayer.cpp
View File

@ -408,6 +408,10 @@ namespace netgen
}
}
BitArray fixed_points(np+1);
fixed_points.Clear();
BitArray moveboundarypoint(np+1);
moveboundarypoint.Clear();
for(SurfaceElementIndex si = 0; si < nse; si++)
{
// copy because surfaceels array will be resized!
@ -436,10 +440,29 @@ namespace netgen
newel.SetIndex(si_map[sel.GetIndex()]);
mesh.AddSurfaceElement(newel);
}
else
{
bool has_moved = false;
for(auto p : sel.PNums())
if(mapto[p].Size())
has_moved = true;
if(has_moved)
for(auto p : sel.PNums())
{
if(!mapto[p].Size())
{
fixed_points.SetBit(p);
if(move_boundaries.Test(sel.GetIndex()))
moveboundarypoint.SetBit(p);
}
}
}
if(move_boundaries.Test(sel.GetIndex()))
for(auto& p : mesh[si].PNums())
if(mapto[p].Size())
p = mapto[p].Last();
{
for(auto& p : mesh[si].PNums())
if(mapto[p].Size())
p = mapto[p].Last();
}
}
for(SegmentIndex sei = 0; sei < nseg; sei++)
@ -453,13 +476,85 @@ namespace netgen
for(ElementIndex ei = 0; ei < ne; ei++)
{
auto& el = mesh[ei];
if(!domains[el.GetIndex()])
auto el = mesh[ei];
ArrayMem<PointIndex,4> fixed;
ArrayMem<PointIndex,4> moved;
bool moved_bnd = false;
for(const auto& p : el.PNums())
{
for(auto& p : el.PNums())
if(fixed_points.Test(p))
fixed.Append(p);
if(mapto[p].Size())
moved.Append(p);
if(moveboundarypoint.Test(p))
moved_bnd = true;
}
bool do_move, do_insert;
if(domains.Test(el.GetIndex()))
{
do_move = fixed.Size() && moved_bnd;
do_insert = do_move;
}
else
{
do_move = !fixed.Size() || moved_bnd;
do_insert = !do_move;
}
if(do_move)
{
for(auto& p : mesh[ei].PNums())
if(mapto[p].Size())
p = mapto[p].Last();
}
if(do_insert)
{
if(el.GetType() != TET)
throw Exception("Boundarylayer only implemented for tets outside yet!");
if(moved.Size() == 2)
{
if(fixed.Size() == 2)
throw Exception("This should not be possible!");
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)
{
PointIndex p1 = moved[0];
for(auto i : Range(blp.heights))
{
Element nel = el;
PointIndex p2 = mapto[moved[0]][i];
for(auto& p : nel.PNums())
{
if(p == moved[0])
p = p2;
if(p == fixed[0])
p = p1;
}
p1 = p2;
mesh.AddVolumeElement(nel);
}
}
}
}
for(auto i : Range(1, fd_old+1))

15
tests/pytest/test_boundarylayer.py
View File

@ -92,3 +92,18 @@ def test_splitted_surface():
mesh = ngs.Mesh(mesh)
assert ngs.Integrate(1, mesh) == pytest.approx(1)
assert ngs.Integrate(1, mesh.Materials("slot")) == pytest.approx(0.4)
@pytest.mark.parametrize("outside", [True, False])
def test_pyramids(outside):
geo = CSGeometry()
box = OrthoBrick((0,0,0), (1,1,1))
plate = OrthoBrick((0.3,0.3,0.4),(0.7,0.7,1)) * Plane((0,0,0.6), (0,0,1)).bc("top")
geo.Add((box-plate).mat("air"))
geo.Add(plate.mat("plate"))
mesh = geo.GenerateMesh()
mesh.BoundaryLayer("top", [0.01], "layer", "plate", outside=outside)
ngs = pytest.importorskip("ngsolve")
mesh = ngs.Mesh(mesh)
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)