mirror of
https://github.com/NGSolve/netgen.git
synced 2024-12-24 21:10:33 +05:00
improve functionality of boundarylayer
This commit is contained in:
parent
7fe4ca9c4f
commit
d752ada2bd
@ -103,7 +103,7 @@ namespace netgen
|
|||||||
function, in order to calculate the effective direction
|
function, in order to calculate the effective direction
|
||||||
in which the prismatic layer should grow
|
in which the prismatic layer should grow
|
||||||
*/
|
*/
|
||||||
Vec<3> GetSurfaceNormal(Mesh & mesh, const Element2d & el)
|
inline Vec<3> GetSurfaceNormal(Mesh & mesh, const Element2d & el)
|
||||||
{
|
{
|
||||||
auto v0 = mesh[el[0]];
|
auto v0 = mesh[el[0]];
|
||||||
auto v1 = mesh[el[1]];
|
auto v1 = mesh[el[1]];
|
||||||
@ -115,12 +115,9 @@ namespace netgen
|
|||||||
return normal;
|
return normal;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/*
|
/*
|
||||||
Philippose Rajan - 11 June 2009
|
Philippose Rajan - 11 June 2009
|
||||||
|
modified by Christopher Lackner Apr 2020
|
||||||
|
|
||||||
Added an initial experimental function for
|
Added an initial experimental function for
|
||||||
generating prismatic boundary layers on
|
generating prismatic boundary layers on
|
||||||
@ -149,10 +146,38 @@ namespace netgen
|
|||||||
PrintMessage(3, "Old NP: ", mesh.GetNP());
|
PrintMessage(3, "Old NP: ", mesh.GetNP());
|
||||||
PrintMessage(3, "Old NSE: ",mesh.GetNSE());
|
PrintMessage(3, "Old NSE: ",mesh.GetNSE());
|
||||||
|
|
||||||
|
map<tuple<int, int, int>, int> domains_to_surf_index;
|
||||||
|
map<tuple<PointIndex, PointIndex>, int> pi_to_edgenr;
|
||||||
|
|
||||||
|
map<int, int> last_layer_surface_index_map;
|
||||||
|
int max_surface_index = mesh.GetNFD();
|
||||||
|
|
||||||
|
int max_edge_nr = -1;
|
||||||
|
for(const auto& seg : mesh.LineSegments())
|
||||||
|
if(seg.edgenr > max_edge_nr)
|
||||||
|
max_edge_nr = seg.edgenr;
|
||||||
|
|
||||||
for(int layer = blp.heights.Size(); layer >= 1; layer--)
|
for(int layer = blp.heights.Size(); layer >= 1; layer--)
|
||||||
{
|
{
|
||||||
PrintMessage(3, "Generating layer: ", layer);
|
PrintMessage(3, "Generating layer: ", layer);
|
||||||
|
|
||||||
|
auto map_surface_index = [&](auto si)
|
||||||
|
{
|
||||||
|
if(last_layer_surface_index_map.find(si) == last_layer_surface_index_map.end())
|
||||||
|
{
|
||||||
|
last_layer_surface_index_map[si] = ++max_surface_index;
|
||||||
|
auto& old_fd = mesh.GetFaceDescriptor(si);
|
||||||
|
int domout = blp.outside ? old_fd.DomainOut() : blp.new_matnrs[layer-1];
|
||||||
|
int domin = blp.outside ? blp.new_matnrs[layer-1] : old_fd.DomainIn();
|
||||||
|
FaceDescriptor fd(max_surface_index-1,
|
||||||
|
domin, domout, -1);
|
||||||
|
fd.SetBCProperty(max_surface_index);
|
||||||
|
mesh.AddFaceDescriptor(fd);
|
||||||
|
return max_surface_index;
|
||||||
|
}
|
||||||
|
return last_layer_surface_index_map[si];
|
||||||
|
};
|
||||||
|
|
||||||
mesh.UpdateTopology();
|
mesh.UpdateTopology();
|
||||||
auto& meshtopo = mesh.GetTopology();
|
auto& meshtopo = mesh.GetTopology();
|
||||||
|
|
||||||
@ -213,7 +238,7 @@ namespace netgen
|
|||||||
}
|
}
|
||||||
|
|
||||||
if (!blp.grow_edges)
|
if (!blp.grow_edges)
|
||||||
for(const auto& sel : mesh.SurfaceElements())
|
for(const auto& sel : mesh.LineSegments())
|
||||||
{
|
{
|
||||||
bndnodes.Clear(sel[0]);
|
bndnodes.Clear(sel[0]);
|
||||||
bndnodes.Clear(sel[1]);
|
bndnodes.Clear(sel[1]);
|
||||||
@ -241,8 +266,6 @@ namespace netgen
|
|||||||
}
|
}
|
||||||
growthvectors[pi] += veci;
|
growthvectors[pi] += veci;
|
||||||
}
|
}
|
||||||
// growthvectors[pi].Normalize();
|
|
||||||
// growthvectors[pi] *= -1.0;
|
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
@ -251,7 +274,6 @@ namespace netgen
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
// Add quad surface elements at edges for surfaces which
|
// Add quad surface elements at edges for surfaces which
|
||||||
// don't have boundary layers
|
// don't have boundary layers
|
||||||
|
|
||||||
@ -271,13 +293,15 @@ namespace netgen
|
|||||||
|
|
||||||
// Only go in if the segment is still active, and if both its
|
// Only go in if the segment is still active, and if both its
|
||||||
// surface index is part of the "hit-list"
|
// surface index is part of the "hit-list"
|
||||||
if(segsel.Test(sei) && blp.surfid.Contains(mesh[sei].si))
|
if(segsel.Test(sei))
|
||||||
|
{
|
||||||
|
if(blp.surfid.Contains(mesh[sei].si))
|
||||||
{
|
{
|
||||||
// clear the bit to indicate that this segment has been processed
|
// clear the bit to indicate that this segment has been processed
|
||||||
segsel.Clear(sei);
|
segsel.Clear(sei);
|
||||||
|
|
||||||
// Find matching segment pair on other surface
|
// Find matching segment pair on other surface
|
||||||
for (SegmentIndex sej = 0; sej < nseg; sej++)
|
for(SegmentIndex sej = 0; sej < nseg; sej++)
|
||||||
{
|
{
|
||||||
PointIndex segpair_p1 = mesh[sej][1];
|
PointIndex segpair_p1 = mesh[sej][1];
|
||||||
PointIndex segpair_p2 = mesh[sej][0];
|
PointIndex segpair_p2 = mesh[sej][0];
|
||||||
@ -300,34 +324,20 @@ namespace netgen
|
|||||||
auto pnt2_elems = meshtopo.GetVertexSurfaceElements(segpair_p2);
|
auto pnt2_elems = meshtopo.GetVertexSurfaceElements(segpair_p2);
|
||||||
|
|
||||||
for(auto pnt1_sei : pnt1_elems)
|
for(auto pnt1_sei : pnt1_elems)
|
||||||
{
|
if(mesh[pnt1_sei].GetIndex() == mesh[sej].si)
|
||||||
const auto& pnt1_sel = mesh[pnt1_sei];
|
|
||||||
for(auto pnt2_sei : pnt2_elems)
|
for(auto pnt2_sei : pnt2_elems)
|
||||||
{
|
if(pnt1_sei == pnt2_sei)
|
||||||
const Element2d & pnt2_sel = mesh.SurfaceElement(pnt2_sei);
|
pnt_commelem = pnt1_sei;
|
||||||
if((pnt1_sel.GetIndex() == mesh[sej].si)
|
|
||||||
&& (pnt2_sel.GetIndex() == mesh[sej].si)
|
|
||||||
&& (pnt1_sei == pnt2_sei))
|
|
||||||
{
|
|
||||||
pnt_commelem = pnt1_sei;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
const Element2d & commsel = mesh.SurfaceElement(pnt_commelem);
|
if(IsInvalid(pnt_commelem))
|
||||||
|
throw Exception("Couldn't find element on other side for " + ToString(segpair_p1) + " to " + ToString(segpair_p2));
|
||||||
|
|
||||||
|
const auto& commsel = mesh[pnt_commelem];
|
||||||
auto surfelem_vect = GetSurfaceNormal(mesh, commsel);
|
auto surfelem_vect = GetSurfaceNormal(mesh, commsel);
|
||||||
if(blp.outside)
|
if(blp.outside)
|
||||||
surfelem_vect *= -1;
|
surfelem_vect *= -1;
|
||||||
|
|
||||||
double surfangle = Angle(growthvectors[segpair_p1],surfelem_vect);
|
double surfangle = Angle(growthvectors[segpair_p1],surfelem_vect);
|
||||||
// remap the segments to the new points
|
|
||||||
if(!blp.outside)
|
|
||||||
{
|
|
||||||
mesh[sei][0] = mapto[seg_p1];
|
|
||||||
mesh[sei][1] = mapto[seg_p2];
|
|
||||||
mesh[sej][1] = mapto[seg_p1];
|
|
||||||
mesh[sej][0] = mapto[seg_p2];
|
|
||||||
}
|
|
||||||
|
|
||||||
if((surfangle < (90 + angleThreshold) * 3.141592 / 180.0)
|
if((surfangle < (90 + angleThreshold) * 3.141592 / 180.0)
|
||||||
&& (surfangle > (90 - angleThreshold) * 3.141592 / 180.0))
|
&& (surfangle > (90 - angleThreshold) * 3.141592 / 180.0))
|
||||||
@ -342,109 +352,256 @@ namespace netgen
|
|||||||
Element2d sel(QUAD);
|
Element2d sel(QUAD);
|
||||||
if(blp.outside)
|
if(blp.outside)
|
||||||
Swap(seg_p1, seg_p2);
|
Swap(seg_p1, seg_p2);
|
||||||
sel.PNum(4) = mapto[seg_p1];
|
sel[0] = seg_p1;
|
||||||
sel.PNum(3) = mapto[seg_p2];
|
sel[1] = seg_p2;
|
||||||
sel.PNum(2) = seg_p2;
|
sel[2] = mapto[seg_p2];
|
||||||
sel.PNum(1) = seg_p1;
|
sel[3] = mapto[seg_p1];
|
||||||
sel.SetIndex(mesh[sej].si);
|
auto domains = make_tuple(commsel.GetIndex(), blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(commsel.GetIndex()).DomainOut());
|
||||||
|
|
||||||
|
if(domains_to_surf_index.find(domains) == domains_to_surf_index.end())
|
||||||
|
{
|
||||||
|
domains_to_surf_index[domains] = ++max_surface_index;
|
||||||
|
domains_to_surf_index[make_tuple(max_surface_index, get<1>(domains), get<2>(domains))] = max_surface_index;
|
||||||
|
FaceDescriptor fd(max_surface_index-1,
|
||||||
|
get<1>(domains),
|
||||||
|
get<2>(domains),
|
||||||
|
-1);
|
||||||
|
fd.SetBCProperty(max_surface_index);
|
||||||
|
mesh.AddFaceDescriptor(fd);
|
||||||
|
mesh.SetBCName(max_surface_index-1,
|
||||||
|
mesh.GetBCName(get<0>(domains)-1));
|
||||||
|
}
|
||||||
|
auto new_index = domains_to_surf_index[domains];
|
||||||
|
sel.SetIndex(new_index);
|
||||||
mesh.AddSurfaceElement(sel);
|
mesh.AddSurfaceElement(sel);
|
||||||
numquads++;
|
numquads++;
|
||||||
|
|
||||||
|
// Add segments
|
||||||
|
Segment seg_1, seg_2;
|
||||||
|
seg_1[0] = mapto[seg_p1];
|
||||||
|
seg_1[1] = seg_p1;
|
||||||
|
seg_2[0] = seg_p2;
|
||||||
|
seg_2[1] = mapto[seg_p2];
|
||||||
|
auto points = make_tuple(seg_p1, mapto[seg_p1]);
|
||||||
|
if(pi_to_edgenr.find(points) == pi_to_edgenr.end())
|
||||||
|
pi_to_edgenr[points] = ++max_edge_nr;
|
||||||
|
seg_1.edgenr = pi_to_edgenr[points];
|
||||||
|
seg_1[2] = PointIndex::INVALID;
|
||||||
|
seg_1.si = new_index;
|
||||||
|
mesh.AddSegment(seg_1);
|
||||||
|
|
||||||
|
points = make_tuple(seg_p2, mapto[seg_p2]);
|
||||||
|
if(pi_to_edgenr.find(points) == pi_to_edgenr.end())
|
||||||
|
pi_to_edgenr[points] = ++max_edge_nr;
|
||||||
|
|
||||||
|
seg_2[2] = PointIndex::INVALID;
|
||||||
|
seg_2.edgenr = pi_to_edgenr[points];
|
||||||
|
seg_2.si = new_index;
|
||||||
|
mesh.AddSegment(seg_2);
|
||||||
|
mesh[sej].si = new_index;
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
for (int k = 0; k < pnt1_elems.Size(); k++)
|
for(auto pnt1_ei : pnt1_elems)
|
||||||
{
|
{
|
||||||
Element2d & pnt_sel = mesh.SurfaceElement(pnt1_elems[k]);
|
auto& pnt_sel = mesh[pnt1_ei];
|
||||||
if(pnt_sel.GetIndex() == mesh[sej].si)
|
if(pnt_sel.GetIndex() == mesh[sej].si)
|
||||||
{
|
{
|
||||||
for(int l = 0; l < pnt_sel.GetNP(); l++)
|
for(auto& pi : pnt_sel.PNums())
|
||||||
{
|
{
|
||||||
if(pnt_sel[l] == segpair_p1)
|
if(pi == segpair_p1)
|
||||||
pnt_sel[l] = mapto[seg_p1];
|
pi = mapto[seg_p1];
|
||||||
else if (pnt_sel[l] == segpair_p2)
|
else if(pi == segpair_p2)
|
||||||
pnt_sel[l] = mapto[seg_p2];
|
pi = mapto[seg_p2];
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
for (int k = 0; k < pnt2_elems.Size(); k++)
|
for(auto sk : pnt2_elems)
|
||||||
{
|
{
|
||||||
Element2d & pnt_sel = mesh.SurfaceElement(pnt2_elems[k]);
|
auto& pnt_sel = mesh[sk];
|
||||||
if(pnt_sel.GetIndex() == mesh[sej].si)
|
if(pnt_sel.GetIndex() == mesh[sej].si)
|
||||||
{
|
{
|
||||||
for(int l = 0; l < pnt_sel.GetNP(); l++)
|
for(auto& p : pnt_sel.PNums())
|
||||||
{
|
{
|
||||||
if(pnt_sel[l] == segpair_p1)
|
if(p == segpair_p1)
|
||||||
pnt_sel[l] = mapto[seg_p1];
|
p = mapto[seg_p1];
|
||||||
else if (pnt_sel[l] == segpair_p2)
|
else if (p == segpair_p2)
|
||||||
pnt_sel[l] = mapto[seg_p2];
|
p = mapto[seg_p2];
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
// }
|
|
||||||
}
|
}
|
||||||
else
|
|
||||||
|
// in last layer insert new segments
|
||||||
|
if(layer == blp.heights.Size())
|
||||||
{
|
{
|
||||||
// If the code comes here, it indicates that we are at
|
Segment s1 = mesh[sei];
|
||||||
// a line segment pair which is at the intersection
|
Segment s2 = mesh[sej];
|
||||||
// of two surfaces, both of which have to grow boundary
|
s1.edgenr = ++max_edge_nr;
|
||||||
// layers.... here too, remapping the segments to the
|
s2.edgenr = max_edge_nr;
|
||||||
// new points is required
|
bool create_it = true;
|
||||||
mesh[sei][0] = mapto[seg_p1];
|
if(blp.surfid.Contains(mesh[sej].si))
|
||||||
mesh[sei][1] = mapto[seg_p2];
|
{
|
||||||
mesh[sej][1] = mapto[seg_p1];
|
if(last_layer_surface_index_map.find(s1.si) != last_layer_surface_index_map.end() &&
|
||||||
mesh[sej][0] = mapto[seg_p2];
|
last_layer_surface_index_map.find(s2.si) != last_layer_surface_index_map.end())
|
||||||
|
// edge already mapped
|
||||||
|
create_it = false;
|
||||||
|
s2.si = map_surface_index(s2.si);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
s2.si = domains_to_surf_index[make_tuple(s2.si,
|
||||||
|
blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(s2.si).DomainOut())];
|
||||||
|
}
|
||||||
|
s1.si = map_surface_index(s1.si);
|
||||||
|
if(create_it)
|
||||||
|
{
|
||||||
|
mesh.AddSegment(s1);
|
||||||
|
mesh.AddSegment(s2);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
// remap the segments to the new points
|
||||||
|
mesh[sei][0] = mapto[mesh[sei][0]];
|
||||||
|
mesh[sei][1] = mapto[mesh[sei][1]];
|
||||||
|
mesh[sej][1] = mapto[mesh[sej][1]];
|
||||||
|
mesh[sej][0] = mapto[mesh[sej][0]];
|
||||||
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
// check if it doesn't contain the other edge as well
|
||||||
|
// and if it doesn't contain both mark them as done and
|
||||||
|
// if necessary map them
|
||||||
|
for(SegmentIndex sej = 0; sej<nseg; sej++)
|
||||||
|
{
|
||||||
|
if(mesh[sej][0] == mesh[sei][1] &&
|
||||||
|
mesh[sej][1] == mesh[sei][0])
|
||||||
|
{
|
||||||
|
if(!blp.surfid.Contains(mesh[sej].si))
|
||||||
|
{
|
||||||
|
segsel.Clear(sei);
|
||||||
|
segsel.Clear(sej);
|
||||||
|
PointIndex mapped_point = PointIndex::INVALID;
|
||||||
|
auto p1 = mesh[sei][0];
|
||||||
|
auto p2 = mesh[sei][1];
|
||||||
|
if(mapto[p1].IsValid())
|
||||||
|
mapped_point = p1;
|
||||||
|
else if(mapto[p2].IsValid())
|
||||||
|
mapped_point = p2;
|
||||||
|
else
|
||||||
|
continue;
|
||||||
|
auto other_point = mapped_point == p1 ? p2 : p1;
|
||||||
|
if(growthvectors[mapped_point] * (mesh[other_point] - mesh[mapped_point]) < 0)
|
||||||
|
{
|
||||||
|
if(mapto[mesh[sei][0]].IsValid())
|
||||||
|
mesh[sei][0] = mapto[mesh[sei][0]];
|
||||||
|
if(mapto[mesh[sei][1]].IsValid())
|
||||||
|
mesh[sei][1] = mapto[mesh[sei][1]];
|
||||||
|
if(mapto[mesh[sej][0]].IsValid())
|
||||||
|
mesh[sej][0] = mapto[mesh[sej][0]];
|
||||||
|
if(mapto[mesh[sej][1]].IsValid())
|
||||||
|
mesh[sej][1] = mapto[mesh[sej][1]];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
else
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
// add surface elements between layer and old domain
|
||||||
|
if(layer == blp.heights.Size())
|
||||||
|
{
|
||||||
|
for(SurfaceElementIndex si = 0; si < nse; si++)
|
||||||
|
{
|
||||||
|
if(blp.surfid.Contains(mesh[si].GetIndex()))
|
||||||
|
{
|
||||||
|
const auto& sel = mesh[si];
|
||||||
|
Element2d newel = sel;
|
||||||
|
newel.SetIndex(map_surface_index(sel.GetIndex()));
|
||||||
|
mesh.AddSurfaceElement(newel);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
// Add prismatic cells at the boundaries
|
// Add prismatic cells at the boundaries
|
||||||
PrintMessage(3, "Generating prism boundary layer volume elements...");
|
PrintMessage(3, "Generating prism boundary layer volume elements...");
|
||||||
|
|
||||||
for (SurfaceElementIndex si = 0; si < nse; si++)
|
for (SurfaceElementIndex si = 0; si < nse; si++)
|
||||||
{
|
{
|
||||||
Element2d & sel = mesh.SurfaceElement(si);
|
const auto& sel = mesh[si];
|
||||||
if(blp.surfid.Contains(sel.GetIndex()))
|
if(blp.surfid.Contains(sel.GetIndex()))
|
||||||
{
|
{
|
||||||
int classify = 0;
|
int classify = 0;
|
||||||
for (int j = 0; j < 3; j++)
|
for(auto j : Range(sel.PNums()))
|
||||||
if (mapto[sel[j]].IsValid())
|
if (mapto[sel[j]].IsValid())
|
||||||
classify += (1 << j);
|
classify += (1 << j);
|
||||||
|
|
||||||
// cout << "classify = " << classify << endl;
|
if(classify == 0)
|
||||||
|
continue;
|
||||||
|
|
||||||
ELEMENT_TYPE types[] = { PRISM, TET, TET, PYRAMID,
|
Element el;
|
||||||
TET, PYRAMID, PYRAMID, PRISM };
|
|
||||||
int nums[] = { sel[0], sel[1], sel[2], mapto[sel[0]], mapto[sel[1]], mapto[sel[2]] };
|
|
||||||
int vertices[][6] =
|
|
||||||
{
|
|
||||||
{ 0, 1, 2, 0, 1, 2 }, // should not occur
|
|
||||||
{ 0, 2, 1, 3, 0, 0 },
|
|
||||||
{ 0, 2, 1, 4, 0, 0 },
|
|
||||||
{ 0, 1, 4, 3, 2, 0 },
|
|
||||||
|
|
||||||
{ 0, 2, 1, 5, 0, 0 },
|
if(sel.GetType() == TRIG)
|
||||||
{ 2, 0, 3, 5, 1, 0 },
|
|
||||||
{ 1, 2, 5, 4, 0, 0 },
|
|
||||||
{ 0, 2, 1, 3, 5, 4 }
|
|
||||||
};
|
|
||||||
if(blp.outside)
|
|
||||||
{
|
{
|
||||||
if(classify != 7)
|
ELEMENT_TYPE types[] = { PRISM, TET, TET, PYRAMID,
|
||||||
throw Exception("Outside with non prisms not yet implemented");
|
TET, PYRAMID, PYRAMID, PRISM };
|
||||||
for(auto i : Range(6))
|
int nums[] = { sel[0], sel[1], sel[2], mapto[sel[0]], mapto[sel[1]], mapto[sel[2]] };
|
||||||
vertices[7][i] = i;
|
int vertices[][6] =
|
||||||
|
{
|
||||||
|
{ 0, 1, 2, 0, 1, 2 }, // should not occur
|
||||||
|
{ 0, 2, 1, 3, 0, 0 },
|
||||||
|
{ 0, 2, 1, 4, 0, 0 },
|
||||||
|
{ 0, 1, 4, 3, 2, 0 },
|
||||||
|
|
||||||
|
{ 0, 2, 1, 5, 0, 0 },
|
||||||
|
{ 2, 0, 3, 5, 1, 0 },
|
||||||
|
{ 1, 2, 5, 4, 0, 0 },
|
||||||
|
{ 0, 2, 1, 3, 5, 4 }
|
||||||
|
};
|
||||||
|
if(blp.outside)
|
||||||
|
{
|
||||||
|
if(classify != 7)
|
||||||
|
throw Exception("Outside with non prisms not yet implemented");
|
||||||
|
for(auto i : Range(6))
|
||||||
|
vertices[7][i] = i;
|
||||||
|
}
|
||||||
|
|
||||||
|
el = Element(types[classify]);
|
||||||
|
for(auto i : Range(el.PNums()))
|
||||||
|
el.PNums()[i] = nums[vertices[classify][i]];
|
||||||
|
}
|
||||||
|
else // sel.GetType() == QUAD
|
||||||
|
{
|
||||||
|
int nums[] = { sel[0], sel[1], sel[2], sel[3],
|
||||||
|
mapto[sel[0]], mapto[sel[1]],
|
||||||
|
mapto[sel[2]], mapto[sel[3]] };
|
||||||
|
if(classify == 15)
|
||||||
|
{
|
||||||
|
int vertices[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
|
||||||
|
if(!blp.outside)
|
||||||
|
{
|
||||||
|
Swap(vertices[1], vertices[3]);
|
||||||
|
Swap(vertices[5], vertices[7]);
|
||||||
|
}
|
||||||
|
el = Element(HEX);
|
||||||
|
for(auto i : Range(el.PNums()))
|
||||||
|
el.PNums()[i] = nums[vertices[i]];
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
throw Exception("This type of quad layer not yet implemented!");
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
Element el(types[classify]);
|
|
||||||
for (int i = 0; i < 6; i++)
|
|
||||||
el[i] = nums[vertices[classify][i]];
|
|
||||||
el.SetIndex(blp.new_matnrs[layer-1]);
|
el.SetIndex(blp.new_matnrs[layer-1]);
|
||||||
if (classify != 0)
|
mesh.AddVolumeElement(el);
|
||||||
mesh.AddVolumeElement(el);
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -455,8 +612,16 @@ namespace netgen
|
|||||||
for(SurfaceElementIndex sei : Range(nse))
|
for(SurfaceElementIndex sei : Range(nse))
|
||||||
{
|
{
|
||||||
auto& sel = mesh[sei];
|
auto& sel = mesh[sei];
|
||||||
if((blp.outside && !blp.surfid.Contains(sel.GetIndex())) ||
|
bool to_move = blp.surfid.Contains(sel.GetIndex());
|
||||||
(!blp.outside && blp.surfid.Contains(sel.GetIndex())))
|
if(blp.domains.Size())
|
||||||
|
{
|
||||||
|
if(blp.outside)
|
||||||
|
to_move |= blp.domains[mesh.GetFaceDescriptor(sel.GetIndex()).DomainIn()];
|
||||||
|
else
|
||||||
|
to_move |= !blp.domains[mesh.GetFaceDescriptor(sel.GetIndex()).DomainIn()];
|
||||||
|
}
|
||||||
|
|
||||||
|
if(to_move)
|
||||||
{
|
{
|
||||||
for(auto& pnum : sel.PNums())
|
for(auto& pnum : sel.PNums())
|
||||||
// Check (Doublecheck) if the corresponding point has a
|
// Check (Doublecheck) if the corresponding point has a
|
||||||
@ -467,17 +632,18 @@ namespace netgen
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
if(blp.outside)
|
for(ElementIndex ei : Range(ne))
|
||||||
for(ElementIndex ei : Range(ne))
|
{
|
||||||
{
|
auto& el = mesh[ei];
|
||||||
auto& el = mesh[ei];
|
bool to_move = blp.outside ? blp.domains[el.GetIndex()] : !blp.domains[el.GetIndex()];
|
||||||
|
if(blp.domains.Size() == 0 || to_move)
|
||||||
for(auto& pnum : el.PNums())
|
for(auto& pnum : el.PNums())
|
||||||
// Check (Doublecheck) if the corresponding point has a
|
// Check (Doublecheck) if the corresponding point has a
|
||||||
// copy available for remapping
|
// copy available for remapping
|
||||||
if(mapto[pnum].IsValid())
|
if(mapto[pnum].IsValid())
|
||||||
// Map the surface elements to the new points
|
// Map the volume elements to the new points
|
||||||
pnum = mapto[pnum];
|
pnum = mapto[pnum];
|
||||||
}
|
}
|
||||||
|
|
||||||
// Lock all the prism points so that the rest of the mesh can be
|
// Lock all the prism points so that the rest of the mesh can be
|
||||||
// optimised without invalidating the entire mesh
|
// optimised without invalidating the entire mesh
|
||||||
@ -494,7 +660,6 @@ namespace netgen
|
|||||||
if(bndnodes.Test(i))
|
if(bndnodes.Test(i))
|
||||||
{
|
{
|
||||||
MeshPoint pointtomove;
|
MeshPoint pointtomove;
|
||||||
|
|
||||||
pointtomove = mesh.Point(i);
|
pointtomove = mesh.Point(i);
|
||||||
mesh.Point(i).SetPoint(pointtomove + layerht * growthvectors[i]);
|
mesh.Point(i).SetPoint(pointtomove + layerht * growthvectors[i]);
|
||||||
}
|
}
|
||||||
@ -505,12 +670,12 @@ namespace netgen
|
|||||||
for(int i=1; i <= mesh.GetNFD(); i++)
|
for(int i=1; i <= mesh.GetNFD(); i++)
|
||||||
{
|
{
|
||||||
auto& fd = mesh.GetFaceDescriptor(i);
|
auto& fd = mesh.GetFaceDescriptor(i);
|
||||||
if(blp.surfid.Contains(fd.SurfNr()))
|
if(blp.surfid.Contains(fd.BCProperty()))
|
||||||
{
|
{
|
||||||
if(blp.outside)
|
if(blp.outside)
|
||||||
fd.SetDomainOut(blp.new_matnrs[0]);
|
fd.SetDomainOut(blp.new_matnrs[blp.new_matnrs.Size()-1]);
|
||||||
else
|
else
|
||||||
fd.SetDomainIn(blp.new_matnrs[0]);
|
fd.SetDomainIn(blp.new_matnrs[blp.new_matnrs.Size()-1]);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -15,6 +15,7 @@ public:
|
|||||||
Array<int> surfid;
|
Array<int> surfid;
|
||||||
Array<double> heights;
|
Array<double> heights;
|
||||||
Array<size_t> new_matnrs;
|
Array<size_t> new_matnrs;
|
||||||
|
BitArray domains;
|
||||||
bool outside = false; // set the boundary layer on the outside
|
bool outside = false; // set the boundary layer on the outside
|
||||||
bool grow_edges = false;
|
bool grow_edges = false;
|
||||||
};
|
};
|
||||||
|
@ -966,12 +966,12 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
|
|||||||
|
|
||||||
auto prismlayers = blp.heights.Size();
|
auto prismlayers = blp.heights.Size();
|
||||||
auto first_new_mat = self.GetNDomains() + 1;
|
auto first_new_mat = self.GetNDomains() + 1;
|
||||||
for(auto i : Range(prismlayers))
|
auto max_dom_nr = first_new_mat;
|
||||||
blp.new_matnrs.Append(first_new_mat + i);
|
|
||||||
if(string* pmaterial = get_if<string>(&material); pmaterial)
|
if(string* pmaterial = get_if<string>(&material); pmaterial)
|
||||||
{
|
{
|
||||||
|
self.SetMaterial(first_new_mat, *pmaterial);
|
||||||
for(auto i : Range(prismlayers))
|
for(auto i : Range(prismlayers))
|
||||||
self.SetMaterial(first_new_mat + i, *pmaterial);
|
blp.new_matnrs.Append(first_new_mat);
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
@ -979,16 +979,41 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
|
|||||||
if(py::len(materials) != prismlayers)
|
if(py::len(materials) != prismlayers)
|
||||||
throw Exception("Length of thicknesses and materials must be same!");
|
throw Exception("Length of thicknesses and materials must be same!");
|
||||||
for(auto i : Range(prismlayers))
|
for(auto i : Range(prismlayers))
|
||||||
self.SetMaterial(first_new_mat + i, materials[i].cast<string>());
|
{
|
||||||
|
self.SetMaterial(first_new_mat+i, materials[i].cast<string>());
|
||||||
|
blp.new_matnrs.Append(first_new_mat + i);
|
||||||
|
}
|
||||||
|
max_dom_nr += prismlayers-1;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
blp.domains.SetSize(max_dom_nr + 1); // one based
|
||||||
|
blp.domains.Clear();
|
||||||
|
if(string* pdomain = get_if<string>(&domain); pdomain)
|
||||||
|
{
|
||||||
|
regex pattern(*pdomain);
|
||||||
|
for(auto i : Range(1, first_new_mat))
|
||||||
|
if(regex_match(self.GetMaterial(i), pattern))
|
||||||
|
blp.domains.SetBit(i);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
auto idomain = *get_if<int>(&domain);
|
||||||
|
blp.domains.SetBit(idomain);
|
||||||
|
}
|
||||||
|
// bits for new domains must be set
|
||||||
|
if(!outside)
|
||||||
|
for(auto i : Range(first_new_mat, max_dom_nr+1))
|
||||||
|
blp.domains.SetBit(i);
|
||||||
|
|
||||||
blp.outside = outside;
|
blp.outside = outside;
|
||||||
blp.grow_edges = grow_edges;
|
blp.grow_edges = grow_edges;
|
||||||
|
|
||||||
GenerateBoundaryLayer (self, blp);
|
GenerateBoundaryLayer (self, blp);
|
||||||
|
self.UpdateTopology();
|
||||||
}, py::arg("boundary"), py::arg("thickness"), py::arg("material"),
|
}, py::arg("boundary"), py::arg("thickness"), py::arg("material"),
|
||||||
py::arg("domain") = 1, py::arg("outside") = false,
|
py::arg("domains") = ".*", py::arg("outside") = false,
|
||||||
py::arg("grow_edges") = false, R"delimiter(
|
py::arg("grow_edges") = false,
|
||||||
|
R"delimiter(
|
||||||
Add boundary layer to mesh.
|
Add boundary layer to mesh.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
@ -1003,8 +1028,8 @@ thickness : float or List[float]
|
|||||||
material : str or List[str]
|
material : str or List[str]
|
||||||
Material name of boundary layer(s).
|
Material name of boundary layer(s).
|
||||||
|
|
||||||
domain : string or int
|
domain : str or int
|
||||||
Add layer into domain specified by name or number.
|
Regexp for domain boundarylayer is going into.
|
||||||
|
|
||||||
outside : bool = False
|
outside : bool = False
|
||||||
If true add the layer on the outside
|
If true add the layer on the outside
|
||||||
|
@ -2,16 +2,30 @@
|
|||||||
import pytest
|
import pytest
|
||||||
from netgen.csg import *
|
from netgen.csg import *
|
||||||
|
|
||||||
|
def GetNSurfaceElements(mesh, boundary):
|
||||||
|
nse_in_layer = 0
|
||||||
|
for el in mesh.Elements2D():
|
||||||
|
print(el.index)
|
||||||
|
if mesh.GetBCName(el.index-1) == boundary:
|
||||||
|
nse_in_layer += 1
|
||||||
|
return nse_in_layer
|
||||||
|
|
||||||
@pytest.mark.parametrize("outside", [True, False])
|
@pytest.mark.parametrize("outside", [True, False])
|
||||||
def test_boundarylayer(outside):
|
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
|
||||||
nse_in_layer = 0
|
|
||||||
layer_surfacenames = ["right", "top"]
|
layer_surfacenames = ["right", "top"]
|
||||||
for el in mesh.Elements2D():
|
|
||||||
if mesh.GetBCName(el.index-1) in layer_surfacenames:
|
|
||||||
nse_in_layer += 1
|
|
||||||
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)
|
||||||
assert mesh.ne == ne_before + 2 * nse_in_layer
|
|
||||||
|
|
||||||
|
should_ne = ne_before + 2 * sum([GetNSurfaceElements(mesh, surf) for surf in layer_surfacenames])
|
||||||
|
assert mesh.ne == should_ne
|
||||||
|
capture = capfd.readouterr()
|
||||||
|
assert not "elements are not matching" in capture.out
|
||||||
|
|
||||||
|
for side in ["front"]:
|
||||||
|
mesh.BoundaryLayer(side, [0.001], "layer", outside=outside, grow_edges=True)
|
||||||
|
should_ne += GetNSurfaceElements(mesh, side)
|
||||||
|
assert mesh.ne == should_ne
|
||||||
|
capture = capfd.readouterr()
|
||||||
|
assert not "elements are not matching" in capture.out
|
||||||
|
|
||||||
|
Loading…
Reference in New Issue
Block a user