mirror of
https://github.com/NGSolve/netgen.git
synced 2024-12-25 05:20:34 +05:00
rewrite create boundarylayer function (now more efficient and stable
and easier)
This commit is contained in:
Christopher Lackner
parent
0d48924392
commit
609cbbcadf
@ -88,167 +88,175 @@ namespace netgen
|
|||||||
cout << "Quads: " << nq << endl;
|
cout << "Quads: " << nq << endl;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
void GenerateBoundaryLayer(Mesh& mesh, const BoundaryLayerParameters& blp)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/*
|
|
||||||
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)
|
|
||||||
{
|
{
|
||||||
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;
|
int max_edge_nr = -1;
|
||||||
for(const auto& seg : mesh.LineSegments())
|
for(const auto& seg : mesh.LineSegments())
|
||||||
if(seg.edgenr > max_edge_nr)
|
if(seg.edgenr > max_edge_nr)
|
||||||
max_edge_nr = seg.edgenr;
|
max_edge_nr = seg.edgenr;
|
||||||
|
|
||||||
for(int layer = blp.heights.Size(); layer >= 1; layer--)
|
int new_mat_nr = mesh.GetNDomains() +1;
|
||||||
{
|
mesh.SetMaterial(new_mat_nr, blp.new_mat);
|
||||||
PrintMessage(3, "Generating layer: ", layer);
|
|
||||||
|
|
||||||
auto map_surface_index = [&](auto si)
|
auto domains = blp.domains;
|
||||||
{
|
if(!blp.outside)
|
||||||
if(last_layer_surface_index_map.find(si) == last_layer_surface_index_map.end())
|
domains.Invert();
|
||||||
{
|
|
||||||
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];
|
|
||||||
};
|
|
||||||
|
|
||||||
mesh.UpdateTopology();
|
mesh.UpdateTopology();
|
||||||
auto& meshtopo = mesh.GetTopology();
|
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 np = mesh.GetNP();
|
||||||
int nse = mesh.GetNSE();
|
|
||||||
int ne = mesh.GetNE();
|
int ne = mesh.GetNE();
|
||||||
|
int nse = mesh.GetNSE();
|
||||||
// Safety measure to ensure no issues with mesh
|
|
||||||
// consistency
|
|
||||||
int nseg = mesh.GetNSeg();
|
int nseg = mesh.GetNSeg();
|
||||||
|
|
||||||
// Indicate which points need to be remapped
|
Array<Array<PointIndex>, PointIndex> mapto(np);
|
||||||
BitArray bndnodes(np+1); // big enough for 1-based array
|
|
||||||
|
|
||||||
// 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);
|
Array<Vec<3>, PointIndex> growthvectors(np);
|
||||||
growthvectors = 0.;
|
growthvectors = 0.;
|
||||||
|
|
||||||
// Bit array to identify all the points belonging
|
Array<double> surfacefacs(mesh.GetNFD()+1);
|
||||||
// to the surface of interest
|
surfacefacs = 0.;
|
||||||
bndnodes.Clear();
|
|
||||||
|
|
||||||
// Run through all the surface elements and mark the points
|
auto getSurfaceNormal = [&mesh] (const Element2d& el)
|
||||||
// belonging to those where a boundary layer has to be created.
|
{
|
||||||
// In addition, also calculate the effective surface normal
|
auto v0 = mesh[el[0]];
|
||||||
// vectors at each of those points to determine the mesh motion
|
return Cross(mesh[el[1]]-v0, mesh[el[2]]-v0).Normalize();
|
||||||
// direction
|
};
|
||||||
PrintMessage(3, "Marking points for remapping...");
|
|
||||||
|
// surface index map
|
||||||
for(const auto& sel : mesh.SurfaceElements())
|
Array<int> si_map(mesh.GetNFD()+1);
|
||||||
if (blp.surfid.Contains(sel.GetIndex()))
|
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);
|
|
||||||
auto& fd = mesh.GetFaceDescriptor(sel.GetIndex());
|
|
||||||
auto domin = fd.DomainIn();
|
|
||||||
auto domout = fd.DomainOut();
|
|
||||||
if(blp.domains.Test(domout) && !blp.domains.Test(domin))
|
|
||||||
n2 *= -1;
|
|
||||||
if(!blp.outside)
|
|
||||||
n2 *= -1;
|
|
||||||
for(auto pi : sel.PNums())
|
for(auto pi : sel.PNums())
|
||||||
{
|
{
|
||||||
// Set the bitarray to indicate that the
|
auto & np = growthvectors[pi];
|
||||||
// point is part of the required set
|
if(np.Length() == 0) { np = n; continue; }
|
||||||
bndnodes.SetBit(pi);
|
auto npn = np * n;
|
||||||
|
auto npnp = np * np;
|
||||||
// Add the surface normal to the already existent one
|
auto nn = n * n;
|
||||||
// (This gives the effective normal direction at corners
|
if(nn-npn*npn/npnp == 0) { np = n; continue; }
|
||||||
// and curved areas)
|
np += (nn - npn)/(nn - npn*npn/npnp) * (n - npn/npnp * np);
|
||||||
|
|
||||||
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);
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// 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
|
// project growthvector on surface for inner angles
|
||||||
for(const auto& sel : mesh.SurfaceElements())
|
if(blp.grow_edges)
|
||||||
// if(!blp.surfid.Contains(sel.GetIndex()))
|
|
||||||
// {
|
|
||||||
if(blp.project_boundaries.Contains(sel.GetIndex()))
|
|
||||||
{
|
{
|
||||||
auto n = GetSurfaceNormal(mesh, sel);
|
for(const auto& sel : mesh.SurfaceElements())
|
||||||
|
if(project_boundaries.Test(sel.GetIndex()))
|
||||||
|
{
|
||||||
|
auto n = getSurfaceNormal(sel);
|
||||||
for(auto i : Range(sel.PNums()))
|
for(auto i : Range(sel.PNums()))
|
||||||
{
|
{
|
||||||
auto pi = sel.PNums()[i];
|
auto pi = sel.PNums()[i];
|
||||||
@ -260,491 +268,205 @@ namespace netgen
|
|||||||
auto v2 = (mesh[prev] - mesh[pi]).Normalize();
|
auto v2 = (mesh[prev] - mesh[pi]).Normalize();
|
||||||
auto v3 = growthvectors[pi];
|
auto v3 = growthvectors[pi];
|
||||||
v3.Normalize();
|
v3.Normalize();
|
||||||
// only project if sel goes in boundarylayer direction
|
if((v1 * v3 > 1e-12) || (v2 * v3 > 1e-12))
|
||||||
if((v1 * v3 < 1e-12) || (v2 * v3 < 1e-12))
|
in_surface_direction.SetBit(sel.GetIndex());
|
||||||
continue;
|
|
||||||
|
|
||||||
auto& g = growthvectors[pi];
|
auto& g = growthvectors[pi];
|
||||||
auto ng = n * g;
|
auto ng = n * g;
|
||||||
auto gg = g * g;
|
auto gg = g * g;
|
||||||
auto nn = n * n;
|
auto nn = n * n;
|
||||||
if(fabs(ng*ng-nn*gg) < 1e-12 || fabs(ng) < 1e-12) continue;
|
// if(fabs(ng*ng-nn*gg) < 1e-12 || fabs(ng) < 1e-12) continue;
|
||||||
auto a = -ng*ng/(ng*ng-nn * gg);
|
auto a = -ng*ng/(ng*ng-nn * gg);
|
||||||
auto b = ng*gg/(ng*ng-nn*gg);
|
auto b = ng*gg/(ng*ng-nn*gg);
|
||||||
g += a*g + b*n;
|
g += a*g + b*n;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
}
|
||||||
if (!blp.grow_edges)
|
else
|
||||||
{
|
{
|
||||||
for(const auto& sel : mesh.LineSegments())
|
for(const auto& seg : mesh.LineSegments())
|
||||||
{
|
{
|
||||||
int count = 0;
|
int count = 0;
|
||||||
for(const auto& sel2 : mesh.LineSegments())
|
for(const auto& seg2 : mesh.LineSegments())
|
||||||
if(((sel[0] == sel2[0] && sel[1] == sel2[1]) || (sel[0] == sel2[1] && sel[1] == sel2[0])) && blp.surfid.Contains(sel2.si))
|
if(((seg[0] == seg2[0] && seg[1] == seg2[1]) || (seg[0] == seg2[1] && seg[1] == seg2[0])) && blp.surfid.Contains(seg2.si))
|
||||||
count++;
|
count++;
|
||||||
if(count == 1)
|
if(count == 1)
|
||||||
{
|
{
|
||||||
bndnodes.Clear(sel[0]);
|
growthvectors[seg[0]] = {0., 0., 0.};
|
||||||
bndnodes.Clear(sel[1]);
|
growthvectors[seg[1]] = {0., 0., 0.};
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
// Add additional points into the mesh structure in order to
|
// insert new points
|
||||||
// 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...");
|
|
||||||
|
|
||||||
for (PointIndex pi = 1; pi <= np; pi++)
|
for (PointIndex pi = 1; pi <= np; pi++)
|
||||||
|
if (growthvectors[pi].Length2() != 0)
|
||||||
{
|
{
|
||||||
if (bndnodes.Test(pi))
|
Point<3> p = mesh[pi];
|
||||||
mapto[pi] = mesh.AddPoint(mesh[pi]);
|
for(auto i : Range(blp.heights))
|
||||||
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++)
|
|
||||||
{
|
{
|
||||||
if(!segsel.Test(sei)) continue;
|
p += blp.heights[i] * growthvectors[pi];
|
||||||
const auto& segi = mesh[sei];
|
mapto[pi].Append(mesh.AddPoint(p));
|
||||||
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;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
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)
|
if(blp.grow_edges)
|
||||||
for(SegmentIndex sei = 0; sei < nseg; sei++)
|
|
||||||
{
|
{
|
||||||
// Only go in if the segment is still active, and if both its
|
for(auto sei : moved_segs)
|
||||||
// surface index is part of the "hit-list"
|
|
||||||
if(segsel.Test(sei))
|
|
||||||
{
|
{
|
||||||
// copy here since we will add segments and this would
|
// copy here since we will add segments and this would
|
||||||
// invalidate a reference!
|
// invalidate a reference!
|
||||||
auto segi = mesh[sei];
|
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];
|
auto segj = mesh[sej];
|
||||||
// Find the segment pair on the neighbouring surface element
|
if(type == 0)
|
||||||
// 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])))
|
|
||||||
{
|
{
|
||||||
// clear bit to indicate that processing of this segment is done
|
Segment s;
|
||||||
segsel.Clear(sej);
|
s[0] = mapto[segj[0]].Last();
|
||||||
|
s[1] = mapto[segj[1]].Last();
|
||||||
// if segj is not in surfel list we nned to add quads
|
s[2] = PointIndex::INVALID;
|
||||||
if(!blp.surfid.Contains(segj.si))
|
auto pair = s[0] < s[1] ? make_pair(s[0], s[1]) : make_pair(s[1], s[0]);
|
||||||
{
|
if(seg2edge.find(pair) == seg2edge.end())
|
||||||
SurfaceElementIndex pnt_commelem;
|
seg2edge[pair] = ++max_edge_nr;
|
||||||
SetInvalid(pnt_commelem);
|
s.edgenr = seg2edge[pair];
|
||||||
|
s.si = si_map[segj.si];
|
||||||
auto pnt1_elems = meshtopo.GetVertexSurfaceElements(segj[0]);
|
mesh.AddSegment(s);
|
||||||
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));
|
|
||||||
}
|
}
|
||||||
auto new_index = domains_to_surf_index[domains];
|
// here we need to grow the quad elements
|
||||||
sel.SetIndex(new_index);
|
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);
|
mesh.AddSurfaceElement(sel);
|
||||||
|
|
||||||
// Add segments
|
// TODO: Too many, would be enough to only add outermost ones
|
||||||
Segment seg_1, seg_2;
|
Segment s1;
|
||||||
seg_1[0] = mapto[seg_p1];
|
s1[0] = p2;
|
||||||
seg_1[1] = seg_p1;
|
s1[1] = p3;
|
||||||
seg_2[0] = seg_p2;
|
s1[2] = PointIndex::INVALID;
|
||||||
seg_2[1] = mapto[seg_p2];
|
auto pair = make_pair(p2, p3);
|
||||||
auto points = make_tuple(seg_p1, mapto[seg_p1]);
|
if(seg2edge.find(pair) == seg2edge.end())
|
||||||
if(pi_to_edgenr.find(points) == pi_to_edgenr.end())
|
seg2edge[pair] = ++max_edge_nr;
|
||||||
pi_to_edgenr[points] = ++max_edge_nr;
|
s1.edgenr = seg2edge[pair];
|
||||||
seg_1.edgenr = pi_to_edgenr[points];
|
s1.si = segj.si;
|
||||||
seg_1[2] = PointIndex::INVALID;
|
mesh.AddSegment(s1);
|
||||||
seg_1.si = new_index;
|
Segment s2;
|
||||||
mesh.AddSegment(seg_1);
|
s2[0] = p4;
|
||||||
|
s2[1] = p1;
|
||||||
points = make_tuple(seg_p2, mapto[seg_p2]);
|
s2[2] = PointIndex::INVALID;
|
||||||
if(pi_to_edgenr.find(points) == pi_to_edgenr.end())
|
pair = make_pair(p1, p4);
|
||||||
pi_to_edgenr[points] = ++max_edge_nr;
|
if(seg2edge.find(pair) == seg2edge.end())
|
||||||
|
seg2edge[pair] = ++max_edge_nr;
|
||||||
seg_2[2] = PointIndex::INVALID;
|
s2.edgenr = seg2edge[pair];
|
||||||
seg_2.edgenr = pi_to_edgenr[points];
|
s2.si = segj.si;
|
||||||
seg_2.si = new_index;
|
mesh.AddSegment(s2);
|
||||||
mesh.AddSegment(seg_2);
|
p1 = p4;
|
||||||
|
p2 = p3;
|
||||||
}
|
}
|
||||||
|
Segment s3;
|
||||||
// in last layer insert new segments
|
s3[0] = p3;
|
||||||
if(layer == blp.heights.Size())
|
s3[1] = p4;
|
||||||
{
|
s3[2] = PointIndex::INVALID;
|
||||||
max_edge_nr++;
|
auto pair = p3 < p4 ? make_pair(p3, p4) : make_pair(p4, p3);
|
||||||
if(!blp.surfid.Contains(segj.si))
|
if(seg2edge.find(pair) == seg2edge.end())
|
||||||
{
|
seg2edge[pair] = ++max_edge_nr;
|
||||||
Segment s3 = segj;
|
s3.edgenr = seg2edge[pair];
|
||||||
s3.si = map_surface_index(segj.si)-1;
|
s3.si = segj.si;
|
||||||
Swap(s3[0], s3[1]);
|
|
||||||
if(blp.outside)
|
|
||||||
{
|
|
||||||
s3[0] = mapto[s3[0]];
|
|
||||||
s3[1] = mapto[s3[1]];
|
|
||||||
}
|
|
||||||
else
|
|
||||||
s3.edgenr = max_edge_nr;
|
|
||||||
mesh.AddSegment(s3);
|
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++)
|
for(SurfaceElementIndex si = 0; si < nse; si++)
|
||||||
{
|
{
|
||||||
const auto& sel = mesh[si];
|
auto& sel = mesh[si];
|
||||||
if(blp.surfid.Contains(sel.GetIndex()))
|
if(si_map[sel.GetIndex()] != -1)
|
||||||
{
|
{
|
||||||
Element2d newel = sel;
|
Array<PointIndex> points(sel.PNums());
|
||||||
newel.SetIndex(map_surface_index(sel.GetIndex()));
|
if(surfacefacs[sel.GetIndex()] > 0) Swap(points[0], points[2]);
|
||||||
mesh.AddSurfaceElement(newel);
|
for(auto j : Range(blp.heights))
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// Add prismatic cells at the boundaries
|
|
||||||
PrintMessage(3, "Generating prism boundary layer volume elements...");
|
|
||||||
|
|
||||||
for (SurfaceElementIndex si = 0; si < nse; si++)
|
|
||||||
{
|
{
|
||||||
const auto& sel = mesh[si];
|
auto eltype = points.Size() == 3 ? PRISM : HEX;
|
||||||
if(blp.surfid.Contains(sel.GetIndex()))
|
Element el(eltype);
|
||||||
{
|
for(auto i : Range(points))
|
||||||
int classify = 0;
|
el[i] = points[i];
|
||||||
for(auto j : Range(sel.PNums()))
|
for(auto i : Range(points))
|
||||||
if (mapto[sel[j]].IsValid())
|
points[i] = mapto[sel.PNums()[i]][j];
|
||||||
classify += (1 << j);
|
if(surfacefacs[sel.GetIndex()] > 0) Swap(points[0], points[2]);
|
||||||
|
for(auto i : Range(points))
|
||||||
if(classify == 0)
|
el[sel.PNums().Size() + i] = points[i];
|
||||||
continue;
|
el.SetIndex(new_mat_nr);
|
||||||
|
|
||||||
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]);
|
|
||||||
mesh.AddVolumeElement(el);
|
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
|
for(SegmentIndex sei = 0; sei < nseg; sei++)
|
||||||
// to the newly added ones
|
|
||||||
PrintMessage(3, "Transferring boundary layer surface elements to new vertex references...");
|
|
||||||
|
|
||||||
for(SurfaceElementIndex sei : Range(nse))
|
|
||||||
{
|
{
|
||||||
auto& sel = mesh[sei];
|
auto& seg = mesh[sei];
|
||||||
if(!blp.surfid.Contains(sel.GetIndex()))
|
if(move_boundaries.Test(seg.si))
|
||||||
{
|
for(auto& p : seg.PNums())
|
||||||
const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex());
|
if(mapto[p].Size())
|
||||||
if(blp.outside &&
|
p = mapto[p].Last();
|
||||||
(!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];
|
|
||||||
}
|
}
|
||||||
|
|
||||||
for(ElementIndex ei : Range(ne))
|
for(ElementIndex ei = 0; ei < ne; ei++)
|
||||||
{
|
{
|
||||||
auto& el = mesh[ei];
|
auto& el = mesh[ei];
|
||||||
// only move the elements on the correct side
|
if(!domains[el.GetIndex()])
|
||||||
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(bndnodes.Test(i))
|
for(auto& p : el.PNums())
|
||||||
{
|
if(mapto[p].Size())
|
||||||
MeshPoint pointtomove;
|
p = mapto[p].Last();
|
||||||
pointtomove = mesh.Point(i);
|
|
||||||
mesh.Point(i).SetPoint(pointtomove + layerht * growthvectors[i]);
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
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(mesh.GetFaceDescriptor(mesh.GetNFD()).DomainIn() == new_mat_nr)
|
||||||
if(blp.surfid.Contains(fd.BCProperty()))
|
mesh.GetFaceDescriptor(i).SetDomainOut(new_mat_nr);
|
||||||
{
|
|
||||||
if(blp.outside)
|
|
||||||
fd.SetDomainOut(blp.new_matnrs[blp.new_matnrs.Size()-1]);
|
|
||||||
else
|
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!");
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|||||||
@ -14,7 +14,7 @@ public:
|
|||||||
// parameters by Philippose ..
|
// parameters by Philippose ..
|
||||||
Array<int> surfid;
|
Array<int> surfid;
|
||||||
Array<double> heights;
|
Array<double> heights;
|
||||||
Array<size_t> new_matnrs;
|
string new_mat;
|
||||||
BitArray domains;
|
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;
|
||||||
|
|||||||
@ -1017,9 +1017,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
|
|||||||
|
|
||||||
.def ("BoundaryLayer", [](Mesh & self, variant<string, int> boundary,
|
.def ("BoundaryLayer", [](Mesh & self, variant<string, int> boundary,
|
||||||
variant<double, py::list> thickness,
|
variant<double, py::list> thickness,
|
||||||
variant<string, py::list> material,
|
string material,
|
||||||
variant<string, int> domain, bool outside,
|
variant<string, int> domain, bool outside,
|
||||||
bool grow_edges,
|
|
||||||
optional<string> project_boundaries)
|
optional<string> project_boundaries)
|
||||||
{
|
{
|
||||||
BoundaryLayerParameters blp;
|
BoundaryLayerParameters blp;
|
||||||
@ -1050,6 +1049,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
blp.new_mat = material;
|
||||||
|
|
||||||
if(project_boundaries.has_value())
|
if(project_boundaries.has_value())
|
||||||
{
|
{
|
||||||
@ -1070,34 +1070,13 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
|
|||||||
blp.heights.Append(val.cast<double>());
|
blp.heights.Append(val.cast<double>());
|
||||||
}
|
}
|
||||||
|
|
||||||
auto prismlayers = blp.heights.Size();
|
int nr_domains = self.GetNDomains();
|
||||||
auto first_new_mat = self.GetNDomains() + 1;
|
blp.domains.SetSize(nr_domains + 1); // one based
|
||||||
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
|
|
||||||
blp.domains.Clear();
|
blp.domains.Clear();
|
||||||
if(string* pdomain = get_if<string>(&domain); pdomain)
|
if(string* pdomain = get_if<string>(&domain); pdomain)
|
||||||
{
|
{
|
||||||
regex pattern(*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))
|
if(regex_match(self.GetMaterial(i), pattern))
|
||||||
blp.domains.SetBit(i);
|
blp.domains.SetBit(i);
|
||||||
}
|
}
|
||||||
@ -1106,19 +1085,15 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
|
|||||||
auto idomain = *get_if<int>(&domain);
|
auto idomain = *get_if<int>(&domain);
|
||||||
blp.domains.SetBit(idomain);
|
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 = true;
|
||||||
|
|
||||||
GenerateBoundaryLayer (self, blp);
|
GenerateBoundaryLayer (self, blp);
|
||||||
self.UpdateTopology();
|
self.UpdateTopology();
|
||||||
}, py::arg("boundary"), py::arg("thickness"), py::arg("material"),
|
}, py::arg("boundary"), py::arg("thickness"), py::arg("material"),
|
||||||
py::arg("domains") = ".*", py::arg("outside") = false,
|
py::arg("domains") = ".*", py::arg("outside") = false,
|
||||||
py::arg("grow_edges") = false, py::arg("project_boundaries")=nullopt,
|
py::arg("project_boundaries")=nullopt,
|
||||||
R"delimiter(
|
R"delimiter(
|
||||||
Add boundary layer to mesh.
|
Add boundary layer to mesh.
|
||||||
|
|
||||||
|
|||||||
@ -18,7 +18,7 @@ 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
|
||||||
layer_surfacenames = ["right", "top", "left", "back", "bottom"]
|
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)
|
should_ne = ne_before + 2 * GetNSurfaceElements(mesh, layer_surfacenames)
|
||||||
assert mesh.ne == should_ne
|
assert mesh.ne == should_ne
|
||||||
@ -26,7 +26,7 @@ def test_boundarylayer(outside, capfd):
|
|||||||
assert not "elements are not matching" in capture.out
|
assert not "elements are not matching" in capture.out
|
||||||
|
|
||||||
for side in ["front"]:
|
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])
|
should_ne += 2 * GetNSurfaceElements(mesh, [side])
|
||||||
assert mesh.ne == should_ne
|
assert mesh.ne == should_ne
|
||||||
capture = capfd.readouterr()
|
capture = capfd.readouterr()
|
||||||
@ -53,8 +53,8 @@ def test_boundarylayer2(outside, version, capfd):
|
|||||||
geo.CloseSurfaces(top, bot, [])
|
geo.CloseSurfaces(top, bot, [])
|
||||||
mesh = geo.GenerateMesh()
|
mesh = geo.GenerateMesh()
|
||||||
should_ne = mesh.ne + 2 * GetNSurfaceElements(mesh, ["default"], "part")
|
should_ne = mesh.ne + 2 * GetNSurfaceElements(mesh, ["default"], "part")
|
||||||
layersize = 0.05
|
layersize = 0.025
|
||||||
mesh.BoundaryLayer("default", [0.5 * layersize, layersize], "part", domains="part", outside=outside, grow_edges=True)
|
mesh.BoundaryLayer("default", [layersize, layersize], "part", domains="part", outside=outside)
|
||||||
assert mesh.ne == should_ne
|
assert mesh.ne == should_ne
|
||||||
assert not "elements are not matching" in capfd.readouterr().out
|
assert not "elements are not matching" in capfd.readouterr().out
|
||||||
import netgen.gui
|
import netgen.gui
|
||||||
@ -73,8 +73,7 @@ def test_wrong_orientation(outside):
|
|||||||
|
|
||||||
mesh = geo.GenerateMesh()
|
mesh = geo.GenerateMesh()
|
||||||
|
|
||||||
mesh.BoundaryLayer(".*", 0.1, "air", domains="air", outside=outside,
|
mesh.BoundaryLayer(".*", 0.1, "air", domains="air", outside=outside)
|
||||||
grow_edges=True)
|
|
||||||
ngs = pytest.importorskip("ngsolve")
|
ngs = pytest.importorskip("ngsolve")
|
||||||
mesh = ngs.Mesh(mesh)
|
mesh = ngs.Mesh(mesh)
|
||||||
assert ngs.Integrate(1, mesh) == pytest.approx(1.2**3 if outside else 1)
|
assert ngs.Integrate(1, mesh) == pytest.approx(1.2**3 if outside else 1)
|
||||||
@ -88,8 +87,7 @@ def test_splitted_surface():
|
|||||||
geo.Add((brick*slots).mat("slot"))
|
geo.Add((brick*slots).mat("slot"))
|
||||||
|
|
||||||
mesh = geo.GenerateMesh()
|
mesh = geo.GenerateMesh()
|
||||||
mesh.BoundaryLayer(".*", [0.001, 0.002], "block", "block", outside=False,
|
mesh.BoundaryLayer(".*", [0.001, 0.001], "block", "block", outside=False)
|
||||||
grow_edges=True)
|
|
||||||
ngs = pytest.importorskip("ngsolve")
|
ngs = pytest.importorskip("ngsolve")
|
||||||
mesh = ngs.Mesh(mesh)
|
mesh = ngs.Mesh(mesh)
|
||||||
assert ngs.Integrate(1, mesh) == pytest.approx(1)
|
assert ngs.Integrate(1, mesh) == pytest.approx(1)
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user