Merge branch 'boundarylayers' into 'master'

Boundarylayers

See merge request jschoeberl/netgen!349
This commit is contained in:
Joachim Schöberl 2020-11-19 20:26:13 +00:00
commit 79542999f5
4 changed files with 417 additions and 651 deletions

View File

@ -88,645 +88,385 @@ namespace netgen
cout << "Quads: " << nq << endl;
}
/*
Philippose Rajan - 11 June 2009
Function to calculate the surface normal at a given
vertex of a surface element, with respect to that
surface element.
This function is used by the boundary layer generation
function, in order to calculate the effective direction
in which the prismatic layer should grow
*/
inline Vec<3> GetSurfaceNormal(Mesh & mesh, const Element2d & el)
void GenerateBoundaryLayer(Mesh& mesh, const BoundaryLayerParameters& blp)
{
auto v0 = mesh[el[0]];
auto v1 = mesh[el[1]];
auto v2 = mesh[el[2]];
Vec<3> vec1 = v1-v0;
Vec<3> vec2 = v2-v0;
Vec<3> normal = Cross(vec1, vec2);
normal.Normalize();
return normal;
}
/*
Philippose Rajan - 11 June 2009
modified by Christopher Lackner Apr 2020
Added an initial experimental function for
generating prismatic boundary layers on
a given set of surfaces.
The number of layers, height of the first layer
and the growth / shrink factor can be specified
by the user
Currently, the layer height is calculated using:
height = h_first_layer * (growth_factor^(num_layers - 1))
*/
void GenerateBoundaryLayer (Mesh & mesh, const BoundaryLayerParameters & blp)
{
PrintMessage(1, "Generating boundary layer...");
PrintMessage(3, "Old NP: ", mesh.GetNP());
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--)
{
PrintMessage(3, "Generating layer: ", layer);
int new_mat_nr = mesh.GetNDomains() +1;
mesh.SetMaterial(new_mat_nr, blp.new_mat);
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();
// -1 surf nr is so that curving does not do anything
FaceDescriptor fd(-1,
domin, domout, -1);
fd.SetBCProperty(max_surface_index);
mesh.AddFaceDescriptor(fd);
mesh.SetBCName(max_surface_index-1,
"mapped_" + old_fd.GetBCName());
return max_surface_index;
}
return last_layer_surface_index_map[si];
};
auto domains = blp.domains;
if(!blp.outside)
domains.Invert();
mesh.UpdateTopology();
auto& meshtopo = mesh.GetTopology();
auto layerht = blp.heights[layer-1];
PrintMessage(5, "Layer Height = ", layerht);
// Need to store the old number of points and
// surface elements because there are new points and
// surface elements being added during the process
int np = mesh.GetNP();
int nse = mesh.GetNSE();
int ne = mesh.GetNE();
// Safety measure to ensure no issues with mesh
// consistency
int nse = mesh.GetNSE();
int nseg = mesh.GetNSeg();
// Indicate which points need to be remapped
BitArray bndnodes(np+1); // big enough for 1-based array
Array<Array<PointIndex>, PointIndex> mapto(np);
// Map of the old points to the new points
Array<PointIndex, PointIndex> mapto(np);
// Growth vectors for the prismatic layer based on
// the effective surface normal at a given point
Array<Vec<3>, PointIndex> growthvectors(np);
growthvectors = 0.;
// Bit array to identify all the points belonging
// to the surface of interest
bndnodes.Clear();
Array<double> surfacefacs(mesh.GetNFD()+1);
surfacefacs = 0.;
// Run through all the surface elements and mark the points
// belonging to those where a boundary layer has to be created.
// In addition, also calculate the effective surface normal
// vectors at each of those points to determine the mesh motion
// direction
PrintMessage(3, "Marking points for remapping...");
for(const auto& sel : mesh.SurfaceElements())
if (blp.surfid.Contains(sel.GetIndex()))
auto getSurfaceNormal = [&mesh] (const Element2d& el)
{
auto v0 = mesh[el[0]];
return Cross(mesh[el[1]]-v0, mesh[el[2]]-v0).Normalize();
};
// surface index map
Array<int> si_map(mesh.GetNFD()+1);
si_map = -1;
int fd_old = mesh.GetNFD();
// create new FaceDescriptors
for(auto i : Range(1, fd_old+1))
{
auto& fd = mesh.GetFaceDescriptor(i);
if(blp.surfid.Contains(i))
{
if(auto isIn = domains.Test(fd.DomainIn()); isIn != domains.Test(fd.DomainOut()))
{
int new_si = mesh.GetNFD()+1;
surfacefacs[i] = isIn ? 1. : -1.;
// -1 surf nr is so that curving does not do anything
FaceDescriptor new_fd(-1, isIn ? new_mat_nr : fd.DomainIn(),
isIn ? fd.DomainOut() : new_mat_nr, -1);
new_fd.SetBCProperty(new_si);
mesh.AddFaceDescriptor(new_fd);
si_map[i] = new_si;
mesh.SetBCName(new_si-1, "mapped_" + fd.GetBCName());
}
}
}
// mark points for remapping
for(const auto& sel : mesh.SurfaceElements())
{
auto n = surfacefacs[sel.GetIndex()] * getSurfaceNormal(sel);
if(n.Length2() != 0.)
{
auto n2 = GetSurfaceNormal(mesh,sel);
if(!blp.outside)
n2 *= -1;
for(auto pi : sel.PNums())
{
// Set the bitarray to indicate that the
// point is part of the required set
bndnodes.SetBit(pi);
// Add the surface normal to the already existent one
// (This gives the effective normal direction at corners
// and curved areas)
auto& n1 = growthvectors[pi];
if(n1.Length() == 0) { n1 = n2; continue; }
auto n1n2 = n1 * n2;
auto n1n1 = n1 * n1;
auto n2n2 = n2 * n2;
if(n2n2 - n1n2*n1n2/n1n1 == 0) { n1 = n2; continue; }
n1 += (n2n2 - n1n2)/(n2n2 - n1n2*n1n2/n1n1) * (n2 - n1n2/n1n1 * n1);
auto & np = growthvectors[pi];
if(np.Length() == 0) { np = n; continue; }
auto npn = np * n;
auto npnp = np * np;
auto nn = n * n;
if(nn-npn*npn/npnp == 0) { np = n; continue; }
np += (nn - npn)/(nn - npn*npn/npnp) * (n - npn/npnp * np);
}
}
}
// Bit array to keep track of segments already processed
BitArray segs_done(nseg);
segs_done.Clear();
// map for all segments with same points
// points to pair of SegmentIndex, int
// int is type of other segment, either:
// 0 == adjacent surface grows layer
// 1 == adjacent surface doesn't grow layer, but layer ends on it
// 2 == adjacent surface is interior surface that ends on layer
// 3 == adjacent surface is exterior surface that ends on layer (not allowed yet)
Array<Array<pair<SegmentIndex, int>>, SegmentIndex> segmap(mesh.GetNSeg());
// moved segments
Array<SegmentIndex> moved_segs;
// boundaries to project endings to
BitArray project_boundaries(fd_old+1);
BitArray move_boundaries(fd_old+1);
project_boundaries.Clear();
move_boundaries.Clear();
Array<SurfaceElementIndex, SegmentIndex> seg2surfel(mesh.GetNSeg());
for(auto si : Range(mesh.SurfaceElements()))
{
NgArray<int> surfeledges;
meshtopo.GetSurfaceElementEdges(si+1, surfeledges);
for(auto edgenr : surfeledges)
for(auto sei : Range(mesh.LineSegments()))
if(meshtopo.GetEdge(sei)+1 == edgenr &&
mesh[sei].si == mesh[si].GetIndex())
seg2surfel[sei] = si;
}
for(auto si : Range(mesh.LineSegments()))
{
if(segs_done[si]) continue;
const auto& segi = mesh[si];
if(si_map[segi.si] == -1) continue;
segs_done.SetBit(si);
segmap[si].Append(make_pair(si, 0));
moved_segs.Append(si);
for(auto sj : Range(mesh.LineSegments()))
{
if(segs_done.Test(sj)) continue;
const auto& segj = mesh[sj];
if((segi[0] == segj[0] && segi[1] == segj[1]) ||
(segi[0] == segj[1] && segi[1] == segj[0]))
{
segs_done.SetBit(sj);
int type;
if(si_map[segj.si] != -1)
type = 0;
else if(const auto& fd = mesh.GetFaceDescriptor(segj.si); domains.Test(fd.DomainIn()) && domains.Test(fd.DomainOut()))
{
type = 2;
if(fd.DomainIn() == 0 || fd.DomainOut() == 0)
project_boundaries.SetBit(segj.si);
}
else if(const auto& fd = mesh.GetFaceDescriptor(segj.si); !domains.Test(fd.DomainIn()) && !domains.Test(fd.DomainOut()))
{
type = 3;
if(fd.DomainIn() == 0 || fd.DomainOut() == 0)
project_boundaries.SetBit(segj.si);
move_boundaries.SetBit(segj.si);
}
else
{
type = 1;
// in case 1 we project the growthvector onto the surface
project_boundaries.SetBit(segj.si);
}
segmap[si].Append(make_pair(sj, type));
}
}
}
BitArray in_surface_direction(fd_old+1);
in_surface_direction.Clear();
// project growthvector on surface for inner angles
// for(const auto& sel : mesh.SurfaceElements())
// if(!blp.surfid.Contains(sel.GetIndex()))
// {
// auto n = GetSurfaceNormal(mesh, sel);
// for(auto pi : sel.PNums())
// {
// if(growthvectors[pi].Length2() == 0.)
// continue;
// auto& g = growthvectors[pi];
// auto ng = n * g;
// auto gg = g * g;
// auto nn = n * n;
// if(fabs(ng*ng-nn*gg) < 1e-12 || fabs(ng) < 1e-12) continue;
// auto a = -ng*ng/(ng*ng-nn * gg);
// auto b = ng*gg/(ng*ng-nn*gg);
// g += a*g + b*n;
// }
// }
if (!blp.grow_edges)
if(blp.grow_edges)
{
for(const auto& sel : mesh.LineSegments())
for(const auto& sel : mesh.SurfaceElements())
if(project_boundaries.Test(sel.GetIndex()))
{
auto n = getSurfaceNormal(sel);
for(auto i : Range(sel.PNums()))
{
auto pi = sel.PNums()[i];
if(growthvectors[pi].Length2() == 0.)
continue;
auto next = sel.PNums()[(i+1)%sel.GetNV()];
auto prev = sel.PNums()[i == 0 ? sel.GetNV()-1 : i-1];
auto v1 = (mesh[next] - mesh[pi]).Normalize();
auto v2 = (mesh[prev] - mesh[pi]).Normalize();
auto v3 = growthvectors[pi];
v3.Normalize();
if((v1 * v3 > 1e-12) || (v2 * v3 > 1e-12))
in_surface_direction.SetBit(sel.GetIndex());
auto& g = growthvectors[pi];
auto ng = n * g;
auto gg = g * g;
auto nn = n * n;
// if(fabs(ng*ng-nn*gg) < 1e-12 || fabs(ng) < 1e-12) continue;
auto a = -ng*ng/(ng*ng-nn * gg);
auto b = ng*gg/(ng*ng-nn*gg);
g += a*g + b*n;
}
}
}
else
{
for(const auto& seg : mesh.LineSegments())
{
int count = 0;
for(const auto& sel2 : mesh.LineSegments())
if(((sel[0] == sel2[0] && sel[1] == sel2[1]) || (sel[0] == sel2[1] && sel[1] == sel2[0])) && blp.surfid.Contains(sel2.si))
for(const auto& seg2 : mesh.LineSegments())
if(((seg[0] == seg2[0] && seg[1] == seg2[1]) || (seg[0] == seg2[1] && seg[1] == seg2[0])) && blp.surfid.Contains(seg2.si))
count++;
if(count == 1)
{
bndnodes.Clear(sel[0]);
bndnodes.Clear(sel[1]);
growthvectors[seg[0]] = {0., 0., 0.};
growthvectors[seg[1]] = {0., 0., 0.};
}
}
}
// Add additional points into the mesh structure in order to
// clone the surface elements.
// Also invert the growth vectors so that they point inwards,
// and normalize them
PrintMessage(3, "Cloning points and calculating growth vectors...");
// insert new points
for (PointIndex pi = 1; pi <= np; pi++)
if (growthvectors[pi].Length2() != 0)
{
if (bndnodes.Test(pi))
mapto[pi] = mesh.AddPoint(mesh[pi]);
else
mapto[pi].Invalidate();
}
// Add quad surface elements at edges for surfaces which
// don't have boundary layers
// Bit array to keep track of segments already processed
BitArray segsel(nseg);
// Set them all to "1" to initially activate all segments
segsel.Set();
// remove double segments (if more than 2 surfaces come together
// in one edge. If one of them is mapped, keep that one and
// map the others to it.
Array<Array<SegmentIndex>> segmap(nseg);
for(SegmentIndex sei = 0; sei < nseg; sei++)
Point<3> p = mesh[pi];
for(auto i : Range(blp.heights))
{
if(!segsel.Test(sei)) continue;
const auto& segi = mesh[sei];
for(SegmentIndex sej = 0; sej < nseg; sej++)
{
if(sej == sei || !segsel.Test(sej)) continue;
const auto& segj = mesh[sej];
if(segi[0] == segj[0] && segi[1] == segj[1])
{
SegmentIndex main, other;
if(blp.surfid.Contains(segi.si))
{ main = sei; other = sej; }
else { main = sej; other = sei; }
segsel.Clear(other);
for(auto& s : segmap[other])
segmap[main].Append(s);
segmap[other].SetSize(0);
segmap[main].Append(other);
if(other == sei) sej = nseg;
}
p += blp.heights[i] * growthvectors[pi];
mapto[pi].Append(mesh.AddPoint(p));
}
}
PrintMessage(3, "Adding 2D Quad elements on required surfaces...");
// add 2d quads on required surfaces
map<pair<PointIndex, PointIndex>, int> seg2edge;
if(blp.grow_edges)
for(SegmentIndex sei = 0; sei < nseg; sei++)
{
// Only go in if the segment is still active, and if both its
// surface index is part of the "hit-list"
if(segsel.Test(sei))
for(auto sei : moved_segs)
{
// copy here since we will add segments and this would
// invalidate a reference!
auto segi = mesh[sei];
if(blp.surfid.Contains(segi.si))
for(auto [sej, type] : segmap[sei])
{
// clear the bit to indicate that this segment has been processed
segsel.Clear(sei);
// Find matching segment pair on other surface
for(SegmentIndex sej = 0; sej < nseg; sej++)
{
// copy here since we will add segments and this would
// invalidate a reference!
auto segj = mesh[sej];
// Find the segment pair on the neighbouring surface element
// Identified by: seg1[0] = seg_pair[1] and seg1[1] = seg_pair[0]
if(segsel.Test(sej) && ((segi[0] == segj[1]) && (segi[1] == segj[0])))
if(type == 0)
{
// clear bit to indicate that processing of this segment is done
segsel.Clear(sej);
// if segj is not in surfel list we nned to add quads
if(!blp.surfid.Contains(segj.si))
{
SurfaceElementIndex pnt_commelem;
SetInvalid(pnt_commelem);
auto pnt1_elems = meshtopo.GetVertexSurfaceElements(segj[0]);
auto pnt2_elems = meshtopo.GetVertexSurfaceElements(segj[1]);
for(auto pnt1_sei : pnt1_elems)
if(mesh[pnt1_sei].GetIndex() == segj.si)
for(auto pnt2_sei : pnt2_elems)
if(pnt1_sei == pnt2_sei)
pnt_commelem = pnt1_sei;
if(IsInvalid(pnt_commelem))
throw Exception("Couldn't find element on other side for " + ToString(segj[0]) + " to " + ToString(segj[1]));
const auto& commsel = mesh[pnt_commelem];
Element2d sel(QUAD);
auto seg_p1 = segi[0];
auto seg_p2 = segi[1];
if(blp.outside)
Swap(seg_p1, seg_p2);
sel[0] = seg_p1;
sel[1] = seg_p2;
sel[2] = mapto[seg_p2];
sel[3] = mapto[seg_p1];
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(-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));
Segment s;
s[0] = mapto[segj[0]].Last();
s[1] = mapto[segj[1]].Last();
s[2] = PointIndex::INVALID;
auto pair = s[0] < s[1] ? make_pair(s[0], s[1]) : make_pair(s[1], s[0]);
if(seg2edge.find(pair) == seg2edge.end())
seg2edge[pair] = ++max_edge_nr;
s.edgenr = seg2edge[pair];
s.si = si_map[segj.si];
mesh.AddSegment(s);
}
auto new_index = domains_to_surf_index[domains];
sel.SetIndex(new_index);
// here we need to grow the quad elements
else if(type == 1)
{
PointIndex pp1 = segj[1];
PointIndex pp2 = segj[0];
if(in_surface_direction.Test(segj.si))
{
Swap(pp1, pp2);
move_boundaries.SetBit(segj.si);
}
PointIndex p1 = pp1;
PointIndex p2 = pp2;
PointIndex p3, p4;
Segment s0;
s0[0] = p1;
s0[1] = p2;
s0[2] = PointIndex::INVALID;
s0.edgenr = segj.edgenr;
s0.si = segj.si;
mesh.AddSegment(s0);
for(auto i : Range(blp.heights))
{
Element2d sel(QUAD);
p3 = mapto[pp2][i];
p4 = mapto[pp1][i];
sel[0] = p1;
sel[1] = p2;
sel[2] = p3;
sel[3] = p4;
sel.SetIndex(segj.si);
mesh.AddSurfaceElement(sel);
// 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);
// TODO: Too many, would be enough to only add outermost ones
Segment s1;
s1[0] = p2;
s1[1] = p3;
s1[2] = PointIndex::INVALID;
auto pair = make_pair(p2, p3);
if(seg2edge.find(pair) == seg2edge.end())
seg2edge[pair] = ++max_edge_nr;
s1.edgenr = seg2edge[pair];
s1.si = segj.si;
mesh.AddSegment(s1);
Segment s2;
s2[0] = p4;
s2[1] = p1;
s2[2] = PointIndex::INVALID;
pair = make_pair(p1, p4);
if(seg2edge.find(pair) == seg2edge.end())
seg2edge[pair] = ++max_edge_nr;
s2.edgenr = seg2edge[pair];
s2.si = segj.si;
mesh.AddSegment(s2);
p1 = p4;
p2 = p3;
}
// in last layer insert new segments
if(layer == blp.heights.Size())
{
max_edge_nr++;
if(!blp.surfid.Contains(segj.si))
{
Segment s3 = segj;
s3.si = map_surface_index(segj.si)-1;
Swap(s3[0], s3[1]);
if(blp.outside)
{
s3[0] = mapto[s3[0]];
s3[1] = mapto[s3[1]];
}
else
s3.edgenr = max_edge_nr;
Segment s3;
s3[0] = p3;
s3[1] = p4;
s3[2] = PointIndex::INVALID;
auto pair = p3 < p4 ? make_pair(p3, p4) : make_pair(p4, p3);
if(seg2edge.find(pair) == seg2edge.end())
seg2edge[pair] = ++max_edge_nr;
s3.edgenr = seg2edge[pair];
s3.si = segj.si;
mesh.AddSegment(s3);
}
Segment s1 = segi;
Segment s2 = segj;
s1.edgenr = max_edge_nr;
s2.edgenr = max_edge_nr;
auto side_surf = domains_to_surf_index[make_tuple(s2.si, blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(s2.si).DomainOut())];
if(blp.surfid.Contains(segj.si))
s2.si = map_surface_index(segj.si);
else
{
if(blp.outside)
{
s2.si = side_surf;
}
else
mesh[sej].si = side_surf;
}
s1.si = map_surface_index(s1.si);
s1.surfnr1 = s1.surfnr2 = s2.surfnr1 = s2.surfnr2 = -1;
mesh.AddSegment(s1);
mesh.AddSegment(s2);
}
segmap.SetSize(mesh.LineSegments().Size());
for(auto sei2 : segmap[sei])
{
auto& s = mesh[sei2];
if(blp.outside && layer == blp.heights.Size())
{
if(blp.surfid.Contains(s.si))
s.si = map_surface_index(s.si);
s.edgenr = max_edge_nr;
}
else
{
s[0] = mapto[s[0]];
s[1] = mapto[s[1]];
}
}
for(auto sej2 : segmap[sej])
{
auto& s = mesh[sej2];
if(blp.outside && layer == blp.heights.Size())
{
if(blp.surfid.Contains(s.si))
s.si = map_surface_index(s.si);
s.edgenr = max_edge_nr;
}
else
{
s[0] = mapto[s[0]];
s[1] = mapto[s[1]];
}
}
// do not use segi (not even with reference, since
// mesh.AddSegment will resize segment array and
// invalidate reference), this is why we copy it!!!
mesh[sei][0] = mapto[segi[0]];
mesh[sei][1] = mapto[segi[1]];
mesh[sej][0] = mapto[segj[0]];
mesh[sej][1] = mapto[segj[1]];
}
}
}
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(segsel.Test(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]];
}
}
}
}
}
}
}
}
// add surface elements between layer and old domain
if(layer == blp.heights.Size())
{
for(SurfaceElementIndex si = 0; si < nse; si++)
{
const auto& sel = mesh[si];
if(blp.surfid.Contains(sel.GetIndex()))
auto& sel = mesh[si];
if(si_map[sel.GetIndex()] != -1)
{
Element2d newel = sel;
newel.SetIndex(map_surface_index(sel.GetIndex()));
mesh.AddSurfaceElement(newel);
}
}
}
// Add prismatic cells at the boundaries
PrintMessage(3, "Generating prism boundary layer volume elements...");
for (SurfaceElementIndex si = 0; si < nse; si++)
Array<PointIndex> points(sel.PNums());
if(surfacefacs[sel.GetIndex()] > 0) Swap(points[0], points[2]);
for(auto j : Range(blp.heights))
{
const auto& sel = mesh[si];
if(blp.surfid.Contains(sel.GetIndex()))
{
int classify = 0;
for(auto j : Range(sel.PNums()))
if (mapto[sel[j]].IsValid())
classify += (1 << j);
if(classify == 0)
continue;
Element el;
if(sel.GetType() == TRIG)
{
ELEMENT_TYPE types[] = { PRISM, TET, TET, PYRAMID,
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 },
{ 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]] };
ArrayMem<int, 8> vertices;
switch(classify)
{
case 6:
{
if(blp.outside)
throw Exception("Type 6 quad outside layer is not yet implemented!");
el = Element(PRISM);
vertices = {0, 1, 5, 3, 2, 6};
break;
}
case 9:
{
if(blp.outside)
throw Exception("Type 9 quad outside layer is not yet implemented!");
el = Element(PRISM);
vertices = { 1, 4, 0, 2, 7, 3 };
break;
}
case 15:
{
vertices = { 0, 1, 2, 3, 4, 5, 6, 7 };
if(!blp.outside)
{
Swap(vertices[1], vertices[3]);
Swap(vertices[5], vertices[7]);
}
el = Element(HEX);
break;
}
default:
throw Exception("Type " + ToString(classify) + " for quad layer not yet implemented!");
}
for(auto i : Range(el.PNums()))
el.PNums()[i] = nums[vertices[i]];
}
el.SetIndex(blp.new_matnrs[layer-1]);
auto eltype = points.Size() == 3 ? PRISM : HEX;
Element el(eltype);
for(auto i : Range(points))
el[i] = points[i];
for(auto i : Range(points))
points[i] = mapto[sel.PNums()[i]][j];
if(surfacefacs[sel.GetIndex()] > 0) Swap(points[0], points[2]);
for(auto i : Range(points))
el[sel.PNums().Size() + i] = points[i];
el.SetIndex(new_mat_nr);
mesh.AddVolumeElement(el);
}
Element2d newel = sel;
for(auto& p : newel.PNums())
p = mapto[p].Last();
newel.SetIndex(si_map[sel.GetIndex()]);
mesh.AddSurfaceElement(newel);
}
if(move_boundaries.Test(sel.GetIndex()))
for(auto& p : sel.PNums())
if(mapto[p].Size())
p = mapto[p].Last();
}
// Finally switch the point indices of the surface elements
// to the newly added ones
PrintMessage(3, "Transferring boundary layer surface elements to new vertex references...");
for(SurfaceElementIndex sei : Range(nse))
for(SegmentIndex sei = 0; sei < nseg; sei++)
{
auto& sel = mesh[sei];
if(!blp.surfid.Contains(sel.GetIndex()))
{
const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex());
if(blp.outside &&
(!blp.domains[fd.DomainIn()] && !blp.domains[fd.DomainOut()]))
continue;
if(!blp.outside &&
(blp.domains[fd.DomainIn()] || blp.domains[fd.DomainOut()]))
continue;
}
for(auto& pnum : sel.PNums())
if(mapto[pnum].IsValid())
pnum = mapto[pnum];
auto& seg = mesh[sei];
if(move_boundaries.Test(seg.si))
for(auto& p : seg.PNums())
if(mapto[p].Size())
p = mapto[p].Last();
}
for(ElementIndex ei : Range(ne))
for(ElementIndex ei = 0; ei < ne; ei++)
{
auto& el = mesh[ei];
// only move the elements on the correct side
if(blp.outside ? blp.domains[el.GetIndex()] : !blp.domains[el.GetIndex()])
for(auto& pnum : el.PNums())
if(mapto[pnum].IsValid())
pnum = mapto[pnum];
}
// Lock all the prism points so that the rest of the mesh can be
// optimised without invalidating the entire mesh
// for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++)
for (PointIndex pi = 1; pi <= np; pi++)
if(bndnodes.Test(pi)) mesh.AddLockedPoint(pi);
// Now, actually pull back the old surface points to create
// the actual boundary layers
PrintMessage(3, "Moving and optimising boundary layer points...");
for (PointIndex i = 1; i <= np; i++)
if(!domains[el.GetIndex()])
{
if(bndnodes.Test(i))
{
MeshPoint pointtomove;
pointtomove = mesh.Point(i);
mesh.Point(i).SetPoint(pointtomove + layerht * growthvectors[i]);
for(auto& p : el.PNums())
if(mapto[p].Size())
p = mapto[p].Last();
}
}
mesh.Compress();
}
for(int i=1; i <= mesh.GetNFD(); i++)
for(auto i : Range(1, fd_old+1))
if(si_map[i] != -1)
{
auto& fd = mesh.GetFaceDescriptor(i);
if(blp.surfid.Contains(fd.BCProperty()))
{
if(blp.outside)
fd.SetDomainOut(blp.new_matnrs[blp.new_matnrs.Size()-1]);
if(mesh.GetFaceDescriptor(mesh.GetNFD()).DomainIn() == new_mat_nr)
mesh.GetFaceDescriptor(i).SetDomainOut(new_mat_nr);
else
fd.SetDomainIn(blp.new_matnrs[blp.new_matnrs.Size()-1]);
mesh.GetFaceDescriptor(i).SetDomainIn(new_mat_nr);
}
}
PrintMessage(3, "New NP: ", mesh.GetNP());
PrintMessage(1, "Boundary Layer Generation....Done!");
}
}

View File

@ -14,10 +14,11 @@ public:
// parameters by Philippose ..
Array<int> surfid;
Array<double> heights;
Array<size_t> new_matnrs;
string new_mat;
BitArray domains;
bool outside = false; // set the boundary layer on the outside
bool grow_edges = false;
Array<size_t> project_boundaries;
};
DLL_HEADER void GenerateBoundaryLayer (Mesh & mesh,

View File

@ -839,6 +839,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
.def("FaceDescriptor", static_cast<FaceDescriptor&(Mesh::*)(int)> (&Mesh::GetFaceDescriptor),
py::return_value_policy::reference)
.def("GetNFaceDescriptors", &Mesh::GetNFD)
.def("GetNDomains", &Mesh::GetNDomains)
.def("GetVolumeNeighboursOfSurfaceElement", [](Mesh & self, size_t sel)
{
@ -1016,9 +1017,9 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
.def ("BoundaryLayer", [](Mesh & self, variant<string, int> boundary,
variant<double, py::list> thickness,
variant<string, py::list> material,
string material,
variant<string, int> domain, bool outside,
bool grow_edges)
optional<string> project_boundaries)
{
BoundaryLayerParameters blp;
if(int* bc = get_if<int>(&boundary); bc)
@ -1040,7 +1041,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
if(dom_pattern)
{
regex pattern(*dom_pattern);
if(regex_match(self.GetMaterial(fd.DomainIn()), pattern) || (fd.DomainOut() > 0 ? regex_match(self.GetMaterial(fd.DomainOut()), pattern) : false))
if((fd.DomainIn() > 0 && regex_match(self.GetMaterial(fd.DomainIn()), pattern)) || (fd.DomainOut() > 0 && regex_match(self.GetMaterial(fd.DomainOut()), pattern)))
blp.surfid.Append(i);
}
else
@ -1048,6 +1049,15 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
}
}
}
blp.new_mat = material;
if(project_boundaries.has_value())
{
regex pattern(*project_boundaries);
for(int i = 1; i<=self.GetNFD(); i++)
if(regex_match(self.GetFaceDescriptor(i).GetBCName(), pattern))
blp.project_boundaries.Append(i);
}
if(double* pthickness = get_if<double>(&thickness); pthickness)
{
@ -1060,34 +1070,13 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
blp.heights.Append(val.cast<double>());
}
auto prismlayers = blp.heights.Size();
auto first_new_mat = self.GetNDomains() + 1;
auto max_dom_nr = first_new_mat;
if(string* pmaterial = get_if<string>(&material); pmaterial)
{
self.SetMaterial(first_new_mat, *pmaterial);
for(auto i : Range(prismlayers))
blp.new_matnrs.Append(first_new_mat);
}
else
{
auto materials = *get_if<py::list>(&material);
if(py::len(materials) != prismlayers)
throw Exception("Length of thicknesses and materials must be same!");
for(auto i : Range(prismlayers))
{
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
int nr_domains = self.GetNDomains();
blp.domains.SetSize(nr_domains + 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))
for(auto i : Range(1, nr_domains+1))
if(regex_match(self.GetMaterial(i), pattern))
blp.domains.SetBit(i);
}
@ -1096,19 +1085,15 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
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.grow_edges = grow_edges;
blp.grow_edges = true;
GenerateBoundaryLayer (self, blp);
self.UpdateTopology();
}, py::arg("boundary"), py::arg("thickness"), py::arg("material"),
py::arg("domains") = ".*", py::arg("outside") = false,
py::arg("grow_edges") = false,
py::arg("project_boundaries")=nullopt,
R"delimiter(
Add boundary layer to mesh.
@ -1133,6 +1118,11 @@ outside : bool = False
grow_edges : bool = False
Grow boundary layer over edges.
project_boundaries : Optional[str] = None
Project boundarylayer to these boundaries if they meet them. Set
to boundaries that meet boundarylayer at a non-orthogonal edge and
layer-ending should be projected to that boundary.
)delimiter")
.def ("EnableTable", [] (Mesh & self, string name, bool set)

View File

@ -18,7 +18,7 @@ def test_boundarylayer(outside, capfd):
mesh = unit_cube.GenerateMesh(maxh=0.3)
ne_before = mesh.ne
layer_surfacenames = ["right", "top", "left", "back", "bottom"]
mesh.BoundaryLayer("|".join(layer_surfacenames), [0.01, 0.02], "layer", outside=outside, grow_edges=True)
mesh.BoundaryLayer("|".join(layer_surfacenames), [0.01, 0.01], "layer", outside=outside)
should_ne = ne_before + 2 * GetNSurfaceElements(mesh, layer_surfacenames)
assert mesh.ne == should_ne
@ -26,7 +26,7 @@ def test_boundarylayer(outside, capfd):
assert not "elements are not matching" in capture.out
for side in ["front"]:
mesh.BoundaryLayer(side, [0.001, 0.002], "layer", outside=outside, grow_edges=True)
mesh.BoundaryLayer(side, [0.001, 0.001], "layer", outside=outside)
should_ne += 2 * GetNSurfaceElements(mesh, [side])
assert mesh.ne == should_ne
capture = capfd.readouterr()
@ -53,7 +53,42 @@ def test_boundarylayer2(outside, version, capfd):
geo.CloseSurfaces(top, bot, [])
mesh = geo.GenerateMesh()
should_ne = mesh.ne + 2 * GetNSurfaceElements(mesh, ["default"], "part")
layersize = 0.05
mesh.BoundaryLayer("default", [0.5 * layersize, layersize], "layer", domains="part", outside=outside, grow_edges=True)
layersize = 0.025
mesh.BoundaryLayer("default", [layersize, layersize], "part", domains="part", outside=outside)
assert mesh.ne == should_ne
assert not "elements are not matching" in capfd.readouterr().out
import netgen.gui
ngs = pytest.importorskip("ngsolve")
ngs.Draw(ngs.Mesh(mesh))
mesh = ngs.Mesh(mesh)
assert ngs.Integrate(1, mesh.Materials("part")) == pytest.approx(0.5*2.05*2.05 if outside else 0.4*2*2)
assert ngs.Integrate(1, mesh) == pytest.approx(3**3)
@pytest.mark.parametrize("outside", [True, False])
def test_wrong_orientation(outside):
geo = CSGeometry()
brick = OrthoBrick((-1,0,0),(1,1,1)) - Plane((0,0,0), (1,0,0))
geo.Add(brick.mat("air"))
mesh = geo.GenerateMesh()
mesh.BoundaryLayer(".*", 0.1, "air", domains="air", outside=outside)
ngs = pytest.importorskip("ngsolve")
mesh = ngs.Mesh(mesh)
assert ngs.Integrate(1, mesh) == pytest.approx(1.2**3 if outside else 1)
def test_splitted_surface():
geo = CSGeometry()
brick = OrthoBrick((0,0,0), (1,1,1))
slots = OrthoBrick((0.2,0,-1), (0.4, 1, 2)) + OrthoBrick((0.6, 0,-1), (0.8, 1,2))
geo.Add((brick-slots).mat("block"))
geo.Add((brick*slots).mat("slot"))
mesh = geo.GenerateMesh()
mesh.BoundaryLayer(".*", [0.001, 0.001], "block", "block", outside=False)
ngs = pytest.importorskip("ngsolve")
mesh = ngs.Mesh(mesh)
assert ngs.Integrate(1, mesh) == pytest.approx(1)
assert ngs.Integrate(1, mesh.Materials("slot")) == pytest.approx(0.4)