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

more stable boundarylayer, also cut prisms at outside

Browse Source
This commit is contained in:
Christopher Lackner 2021-03-29 14:02:00 +02:00
parent c27fa6899b
commit d2dc84b02c
2 changed files with 163 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

182
libsrc/meshing/boundarylayer.cpp
View File

@ -414,14 +414,14 @@ namespace netgen
} }
} }
BitArray fixed_points(np+1); // 1 if moved away from point, 2 if moved towards point
fixed_points.Clear(); Array<Array<tuple<PointIndex, bool>>> fixed_points(np+1);
BitArray moveboundarypoint(np+1);
moveboundarypoint.Clear();
for(SurfaceElementIndex si = 0; si < nse; si++) for(SurfaceElementIndex si = 0; si < nse; si++)
{ {
// copy because surfaceels array will be resized! // copy because surfaceels array will be resized!
auto sel = mesh[si]; auto sel = mesh[si];
// force move this sel
bool move_sel = false;
if(si_map[sel.GetIndex()] != -1) if(si_map[sel.GetIndex()] != -1)
{ {
Array<PointIndex> points(sel.PNums()); Array<PointIndex> points(sel.PNums());
@ -448,22 +448,69 @@ namespace netgen
} }
else else
{ {
bool has_moved = false; ArrayMem<tuple<PointIndex, Vec<3>>, 4> moved_direction;
for(auto p : sel.PNums()) for(auto p : sel.PNums())
if(mapto[p].Size()) if(mapto[p].Size())
has_moved = true; {
if(has_moved) Vec<3> dir = growthvectors[p];
moved_direction.Append({ p, dir.Normalize() });
}
const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex());
bool in_inside = blp.domains.Test(fd.DomainIn());
bool sel_between_in_and_out = in_inside != blp.domains.Test(fd.DomainOut());
if(moved_direction.Size() && sel_between_in_and_out)
for(auto p : sel.PNums()) for(auto p : sel.PNums())
{ {
if(!mapto[p].Size()) if(!mapto[p].Size())
{ {
fixed_points.SetBit(p); tuple<PointIndex, double> max = { PointIndex::INVALID, 0 };
if(move_boundaries.Test(sel.GetIndex())) for(auto& [po, dir] : moved_direction)
moveboundarypoint.SetBit(p); {
Vec<3> mydir = mesh[p] - mesh[po];
mydir.Normalize();
double ang = mydir * dir;
if(fabs(ang) >= fabs(get<1>(max)))
max = { po, ang };
}
fixed_points[p].Append({ get<0>(max), get<1>(max) > 0 });
} }
} }
// if only one point is moved and moved into sel
// then sel needs to be divided
if(!sel_between_in_and_out && !in_inside && moved_direction.Size() == 1)
{
for(const auto& p : sel.PNums())
{
if(!mapto[p].Size())
{
const auto& [po, dir] = moved_direction[0];
Vec<3> mydir = mesh[p] - mesh[po];
mydir.Normalize();
if(mydir * dir > 1 - 1e-6)
{
PointIndex p1 = po;
for(auto i : Range(blp.heights))
{
PointIndex p2 = mapto[po][i];
Element2d newel = sel;
for(auto& pi : newel.PNums())
{
if(pi == po)
pi = p1;
if(pi == p)
pi = p2;
}
p1 = p2;
newel.SetIndex(sel.GetIndex());
mesh.AddSurfaceElement(newel);
move_sel = true;
}
}
}
}
}
} }
if(move_boundaries.Test(sel.GetIndex())) if(move_sel || move_boundaries.Test(sel.GetIndex()))
{ {
for(auto& p : mesh[si].PNums()) for(auto& p : mesh[si].PNums())
if(mapto[p].Size()) if(mapto[p].Size())
@ -484,28 +531,65 @@ namespace netgen
{ {
auto el = mesh[ei]; auto el = mesh[ei];
ArrayMem<PointIndex,4> fixed; ArrayMem<PointIndex,4> fixed;
ArrayMem<PointIndex,4> fixed_other;
ArrayMem<PointIndex,4> moved; ArrayMem<PointIndex,4> moved;
bool moved_bnd = false; bool moved_bnd = false;
for(const auto& p : el.PNums()) for(const auto& p : el.PNums())
{ {
if(fixed_points.Test(p)) if(domains.Test(el.GetIndex()))
fixed.Append(p); {
if(fixed_points[p].Size())
{
bool found = false;
bool found_other = false;
for(const auto& tup : fixed_points[p])
if(get<1>(tup))
{
found_other = true;
if(el.PNums().Contains(get<0>(tup)))
found = true;
}
if(found_other)
fixed_other.Append(p);
if(found)
fixed.Append(p);
}
}
if(!domains.Test(el.GetIndex()))
{
if(fixed_points[p].Size())
{
bool found = false;
bool found_other = false;
for(const auto& tup : fixed_points[p])
{
if(!get<1>(tup))
{
found_other = true;
if(el.PNums().Contains(get<0>(tup)))
found = true;
}
}
if(found)
fixed.Append(p);
if(found_other)
fixed_other.Append(p);
}
}
if(mapto[p].Size()) if(mapto[p].Size())
moved.Append(p); moved.Append(p);
if(moveboundarypoint.Test(p))
moved_bnd = true;
} }
bool do_move, do_insert; bool do_move, do_insert;
if(domains.Test(el.GetIndex())) if(domains.Test(el.GetIndex()))
{ {
do_move = fixed.Size() && moved_bnd; do_move = fixed.Size() && moved.Size();
do_insert = do_move; do_insert = do_move;
} }
else else
{ {
do_move = !fixed.Size() || moved_bnd; do_move = !fixed.Size() && moved.Size();
do_insert = !do_move; do_insert = !do_move && moved.Size();
} }
if(do_move) if(do_move)
@ -576,9 +660,9 @@ namespace netgen
for(auto& p : nel.PNums()) for(auto& p : nel.PNums())
{ {
if(p == moved[0]) if(p == moved[0])
p = p1; p = domains.Test(el.GetIndex()) ? p1 : p2;
else if(p == fixed[0]) else if(p == fixed[0])
p = p2; p = domains.Test(el.GetIndex()) ? p2 : p1;
} }
p1 = p2; p1 = p2;
mesh.AddVolumeElement(nel); mesh.AddVolumeElement(nel);
@ -614,6 +698,64 @@ namespace netgen
else if(moved.Size() == 1) else if(moved.Size() == 1)
throw Exception("This case is not implemented yet!"); throw Exception("This case is not implemented yet!");
} }
else if(el.GetType() == PRISM)
{
if(moved.Size() == 2)
{
auto pnums = el.PNums();
bool cut_vertical = (moved[0] == pnums[0] && moved[1] == pnums[3]) || (moved[0] == pnums[1] && moved[1] == pnums[4]) || (moved[0] == pnums[2] && moved[1] == pnums[5]);
if(cut_vertical)
{
if(fixed.Size() != 2)
{
cout << "in domain = " << domains.Test(el.GetIndex()) << endl;
cout << "moved = " << moved << endl;
cout << "fixed = " << fixed << endl;
cout << "pnums = " << el.PNums() << endl;
throw Exception("Cut prism vert only implemented for fixed size == 2!");
}
auto bot_trig = pnums.Range(3);
if(bot_trig.Contains(moved[1]))
Swap(moved[0], moved[1]);
if(bot_trig.Contains(fixed[1]))
Swap(fixed[0], fixed[1]);
PointIndex tip_bot = pnums[0] + pnums[1] + pnums[2] - moved[0] - fixed[0];
PointIndex tip_top = pnums[3] + pnums[4] + pnums[5] - moved[1] - fixed[1];
int index_moved = bot_trig.Pos(moved[0]);
int index_tip = bot_trig.Pos(tip_bot);
int index_fixed = 3 - index_moved - index_tip;
PointIndex p1_bot = moved[0];
PointIndex p1_top = moved[1];
for(auto i : Range(blp.heights))
{
PointIndex p2_bot = mapto[moved[0]][i];
PointIndex p2_top = mapto[moved[1]][i];
Element nel(PRISM);
nel[index_moved] = p1_bot;
nel[index_tip] = tip_bot;
nel[index_fixed] = p2_bot;
nel[3+index_moved] = p1_top;
nel[3+index_tip] = tip_top;
nel[3+index_fixed] = p2_top;
nel.SetIndex(el.GetIndex());
mesh.AddVolumeElement(nel);
p1_bot = p2_bot;
p1_top = p2_top;
}
}
else
throw Exception("Boundarylayer horizontally threw PRISM not implemented yet!");
}
else
{
cout << "in domain = " << domains.Test(el.GetIndex()) << endl;
cout << "moved = " << moved << endl;
cout << "fixed = " << fixed << endl;
cout << "pnums = " << el.PNums() << endl;
throw Exception("Prism with moved points != 2 not implemented yet!");
}
}
else else
throw Exception("Boundarylayer only implemented for tets and pyramids outside yet!"); throw Exception("Boundarylayer only implemented for tets and pyramids outside yet!");
} }

2
tests/pytest/test_boundarylayer.py
View File

@ -33,7 +33,7 @@ def test_boundarylayer(outside, capfd):
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("outside", [True, False])
@pytest.mark.parametrize("version", [1, 2]) # version 2 not working yet @pytest.mark.parametrize("version", [1, 2])
def test_boundarylayer2(outside, version, capfd): def test_boundarylayer2(outside, version, capfd):
geo = CSGeometry() geo = CSGeometry()
top = Plane(Pnt(0,0,0.5), Vec(0,0,1)) top = Plane(Pnt(0,0,0.5), Vec(0,0,1))