ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-01-12 14:10:34 +05:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'boundarylayer' into 'master'

Browse Source 
modernize and improve GenerateBoundaryLayer

See merge request jschoeberl/netgen!318
This commit is contained in:
Joachim Schöberl 2020-05-13 18:01:34 +00:00
commit 319b6dc600
6 changed files with 617 additions and 559 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

630
libsrc/meshing/boundarylayer.cpp
View File

@ -103,32 +103,21 @@ 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
*/ */
void GetSurfaceNormal(Mesh & mesh, const Element2d & el, int Vertex, Vec3d & SurfaceNormal) inline Vec<3> GetSurfaceNormal(Mesh & mesh, const Element2d & el)
{ {
int Vertex_A; auto v0 = mesh[el[0]];
int Vertex_B; auto v1 = mesh[el[1]];
auto v2 = mesh[el[2]];
Vertex_A = Vertex + 1; Vec<3> vec1 = v1-v0;
if(Vertex_A > el.GetNP()) Vertex_A = 1; Vec<3> vec2 = v2-v0;
Vec<3> normal = Cross(vec1, vec2);
Vertex_B = Vertex - 1; normal.Normalize();
if(Vertex_B <= 0) Vertex_B = el.GetNP(); return normal;
Vec3d Vect_A,Vect_B;
Vect_A = mesh[el.PNum(Vertex_A)] - mesh[el.PNum(Vertex)];
Vect_B = mesh[el.PNum(Vertex_B)] - mesh[el.PNum(Vertex)];
SurfaceNormal = Cross(Vect_A,Vect_B);
SurfaceNormal.Normalize();
} }
/* /*
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
@ -141,62 +130,50 @@ namespace netgen
Currently, the layer height is calculated using: Currently, the layer height is calculated using:
height = h_first_layer * (growth_factor^(num_layers - 1)) height = h_first_layer * (growth_factor^(num_layers - 1))
*/ */
void GenerateBoundaryLayer (Mesh & mesh, BoundaryLayerParameters & blp) void GenerateBoundaryLayer (Mesh & mesh, const BoundaryLayerParameters & blp)
{ {
ofstream dbg("BndLayerDebug.log"); PrintMessage(1, "Generating boundary layer...");
PrintMessage(3, "Old NP: ", mesh.GetNP());
PrintMessage(3, "Old NSE: ",mesh.GetNSE());
// Angle between a surface element and a growth-vector below which map<tuple<int, int, int>, int> domains_to_surf_index;
// a prism is project onto that surface as a quad map<tuple<PointIndex, PointIndex>, int> pi_to_edgenr;
// (in degrees)
double angleThreshold = 5.0;
map<int, int> last_layer_surface_index_map;
int max_surface_index = mesh.GetNFD();
NgArray<int> surfid (blp.surfid); int max_edge_nr = -1;
int prismlayers = blp.prismlayers; for(const auto& seg : mesh.LineSegments())
double hfirst = blp.hfirst; if(seg.edgenr > max_edge_nr)
double growthfactor = blp.growthfactor; max_edge_nr = seg.edgenr;
NgArray<double> heights (blp.heights);
bool grow_edges = false; // grow layer at edges for(int layer = blp.heights.Size(); layer >= 1; layer--)
// Monitor and print out the number of prism and quad elements
// added to the mesh
int numprisms = 0;
int numquads = 0;
cout << "Old NP: " << mesh.GetNP() << endl;
cout << "Old NSE: " << mesh.GetNSE() << endl;
for(int layer = prismlayers; layer >= 1; layer--)
{ {
cout << "Generating layer: " << layer << endl; PrintMessage(3, "Generating layer: ", layer);
const MeshTopology& meshtopo = mesh.GetTopology(); auto map_surface_index = [&](auto si)
const_cast<MeshTopology &> (meshtopo).SetBuildEdges(true);
const_cast<MeshTopology &> (meshtopo).SetBuildFaces(true);
const_cast<MeshTopology &> (meshtopo).Update();
double layerht = hfirst;
if(heights.Size()>0)
{ {
layerht = heights[layer-1]; if(last_layer_surface_index_map.find(si) == last_layer_surface_index_map.end())
} {
else last_layer_surface_index_map[si] = ++max_surface_index;
{ auto& old_fd = mesh.GetFaceDescriptor(si);
if(growthfactor == 1) int domout = blp.outside ? old_fd.DomainOut() : blp.new_matnrs[layer-1];
{ int domin = blp.outside ? blp.new_matnrs[layer-1] : old_fd.DomainIn();
layerht = layer * hfirst; FaceDescriptor fd(max_surface_index-1,
} domin, domout, -1);
else fd.SetBCProperty(max_surface_index);
{ mesh.AddFaceDescriptor(fd);
layerht = hfirst*(pow(growthfactor,(layer+1)) - 1)/(growthfactor - 1); return max_surface_index;
}
} }
return last_layer_surface_index_map[si];
};
cout << "Layer Height = " << layerht << endl; 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 // Need to store the old number of points and
// surface elements because there are new points and // surface elements because there are new points and
@ -210,14 +187,15 @@ namespace netgen
int nseg = mesh.GetNSeg(); int nseg = mesh.GetNSeg();
// Indicate which points need to be remapped // Indicate which points need to be remapped
NgBitArray bndnodes(np+1); // big enough for 1-based array BitArray bndnodes(np+1); // big enough for 1-based array
// Map of the old points to the new points // Map of the old points to the new points
NgArray<PointIndex, PointIndex::BASE> mapto(np); Array<PointIndex, PointIndex> mapto(np);
// Growth vectors for the prismatic layer based on // Growth vectors for the prismatic layer based on
// the effective surface normal at a given point // the effective surface normal at a given point
NgArray<Vec3d, PointIndex::BASE> growthvectors(np); Array<Vec<3>, PointIndex> growthvectors(np);
growthvectors = 0.;
// Bit array to identify all the points belonging // Bit array to identify all the points belonging
// to the surface of interest // to the surface of interest
@ -228,86 +206,83 @@ namespace netgen
// In addition, also calculate the effective surface normal // In addition, also calculate the effective surface normal
// vectors at each of those points to determine the mesh motion // vectors at each of those points to determine the mesh motion
// direction // direction
cout << "Marking points for remapping...." << endl; PrintMessage(3, "Marking points for remapping...");
for (SurfaceElementIndex si = 0; si < nse; si++) for(const auto& sel : mesh.SurfaceElements())
if (surfid.Contains(mesh[si].GetIndex())) if (blp.surfid.Contains(sel.GetIndex()))
{ {
const Element2d & sel = mesh[si]; auto n2 = GetSurfaceNormal(mesh,sel);
for(int j = 0; j < sel.GetNP(); j++) if(!blp.outside)
n2 *= -1;
for(auto pi : sel.PNums())
{ {
// Set the bitarray to indicate that the // Set the bitarray to indicate that the
// point is part of the required set // point is part of the required set
bndnodes.Set(sel[j]); bndnodes.SetBit(pi);
Vec3d surfacenormal;
// Calculate the surface normal at the current point
// with respect to the current surface element
GetSurfaceNormal(mesh,sel,j+1,surfacenormal);
// Add the surface normal to the already existent one // Add the surface normal to the already existent one
// (This gives the effective normal direction at corners // (This gives the effective normal direction at corners
// and curved areas) // and curved areas)
growthvectors[sel[j]] += surfacenormal;
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);
} }
} }
if (!grow_edges) if (!blp.grow_edges)
for (SegmentIndex sei = 0; sei <= nseg; sei++) for(const auto& sel : mesh.LineSegments())
{ {
bndnodes.Clear (mesh[sei][0]); bndnodes.Clear(sel[0]);
bndnodes.Clear (mesh[sei][1]); bndnodes.Clear(sel[1]);
} }
// Add additional points into the mesh structure in order to // Add additional points into the mesh structure in order to
// clone the surface elements. // clone the surface elements.
// Also invert the growth vectors so that they point inwards, // Also invert the growth vectors so that they point inwards,
// and normalize them // and normalize them
cout << "Cloning points and calculating growth vectors...." << endl; PrintMessage(3, "Cloning points and calculating growth vectors...");
for (PointIndex pi = 1; pi <= np; pi++) for (PointIndex pi = 1; pi <= np; pi++)
{ {
if (bndnodes.Test(pi)) if (bndnodes.Test(pi))
{ mapto[pi] = mesh.AddPoint(mesh[pi]);
mapto[pi] = mesh.AddPoint (mesh[pi]);
growthvectors[pi].Normalize();
growthvectors[pi] *= -1.0;
}
else else
{ mapto[pi].Invalidate();
mapto[pi] = 0;
growthvectors[pi] = Vec3d(0,0,0);
} }
}
// 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
// Bit array to keep track of segments already processed // Bit array to keep track of segments already processed
NgBitArray segsel(nseg); BitArray segsel(nseg);
// Set them all to "1" to initially activate all segments // Set them all to "1" to initially activate all segments
segsel.Set(); segsel.Set();
cout << "Adding 2D Quad elements on required surfaces...." << endl; PrintMessage(3, "Adding 2D Quad elements on required surfaces...");
if (grow_edges) if(blp.grow_edges)
for (SegmentIndex sei = 0; sei <= nseg; sei++) for(SegmentIndex sei = 0; sei < nseg; sei++)
{ {
PointIndex seg_p1 = mesh[sei][0]; PointIndex seg_p1 = mesh[sei][0];
PointIndex seg_p2 = mesh[sei][1]; PointIndex seg_p2 = mesh[sei][1];
// 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) && 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];
@ -321,184 +296,195 @@ namespace netgen
// Only worry about those surfaces which are not in the // Only worry about those surfaces which are not in the
// boundary layer list // boundary layer list
if(!surfid.Contains(mesh[sej].si)) if(!blp.surfid.Contains(mesh[sej].si))
{ {
SurfaceElementIndex pnt_commelem = 0; SurfaceElementIndex pnt_commelem;
NgArray<SurfaceElementIndex> pnt1_elems; SetInvalid(pnt_commelem);
NgArray<SurfaceElementIndex> pnt2_elems;
auto pnt1_elems = meshtopo.GetVertexSurfaceElements(segpair_p1);
auto pnt2_elems = meshtopo.GetVertexSurfaceElements(segpair_p2);
meshtopo.GetVertexSurfaceElements(segpair_p1,pnt1_elems); for(auto pnt1_sei : pnt1_elems)
meshtopo.GetVertexSurfaceElements(segpair_p2,pnt2_elems); if(mesh[pnt1_sei].GetIndex() == mesh[sej].si)
for(auto pnt2_sei : pnt2_elems)
if(pnt1_sei == pnt2_sei)
pnt_commelem = pnt1_sei;
for(int k = 0; k < pnt1_elems.Size(); k++) 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);
if(blp.outside)
surfelem_vect *= -1;
Element2d sel(QUAD);
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())
{ {
const Element2d & pnt1_sel = mesh.SurfaceElement(pnt1_elems[k]); domains_to_surf_index[domains] = ++max_surface_index;
for(int l = 0; l < pnt2_elems.Size(); l++) domains_to_surf_index[make_tuple(max_surface_index, get<1>(domains), get<2>(domains))] = max_surface_index;
{ FaceDescriptor fd(max_surface_index-1,
const Element2d & pnt2_sel = mesh.SurfaceElement(pnt2_elems[l]); get<1>(domains),
if((pnt1_sel.GetIndex() == mesh[sej].si) get<2>(domains),
&& (pnt2_sel.GetIndex() == mesh[sej].si) -1);
&& (pnt1_elems[k] == pnt2_elems[l])) fd.SetBCProperty(max_surface_index);
{ mesh.AddFaceDescriptor(fd);
pnt_commelem = pnt1_elems[k]; 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);
// 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;
} }
/* // in last layer insert new segments
int pnum_commelem = 0; if(layer == blp.heights.Size())
for(int k = 1; k <= mesh.SurfaceElement(pnt_commelem).GetNP(); k++)
{ {
if((mesh.SurfaceElement(pnt_commelem).PNum(k) != segpair_p1) Segment s1 = mesh[sei];
&& (mesh.SurfaceElement(pnt_commelem).PNum(k) != segpair_p2)) Segment s2 = mesh[sej];
s1.edgenr = ++max_edge_nr;
s2.edgenr = max_edge_nr;
bool create_it = true;
if(blp.surfid.Contains(mesh[sej].si))
{ {
pnum_commelem = mesh.SurfaceElement(pnt_commelem).PNum(k); if(last_layer_surface_index_map.find(s1.si) != last_layer_surface_index_map.end() &&
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);
} }
} }
*/
Vec3d surfelem_vect, surfelem_vect1;
const Element2d & commsel = mesh.SurfaceElement(pnt_commelem);
dbg << "NP= " << commsel.GetNP() << " : ";
for(int k = 1; k <= commsel.GetNP(); k++)
{
GetSurfaceNormal(mesh,commsel,k,surfelem_vect1);
surfelem_vect += surfelem_vect1;
}
surfelem_vect.Normalize();
double surfangle = Angle(growthvectors.Elem(segpair_p1),surfelem_vect);
dbg << "V1= " << surfelem_vect1
<< " : V2= " << surfelem_vect1
<< " : V= " << surfelem_vect
<< " : GV= " << growthvectors.Elem(segpair_p1)
<< " : Angle= " << surfangle * 180 / 3.141592;
// remap the segments to the new points // remap the segments to the new points
mesh[sei][0] = mapto[seg_p1]; mesh[sei][0] = mapto[mesh[sei][0]];
mesh[sei][1] = mapto[seg_p2]; mesh[sei][1] = mapto[mesh[sei][1]];
mesh[sej][1] = mapto[seg_p1]; mesh[sej][1] = mapto[mesh[sej][1]];
mesh[sej][0] = mapto[seg_p2]; mesh[sej][0] = mapto[mesh[sej][0]];
if((surfangle < (90 + angleThreshold) * 3.141592 / 180.0) }
&& (surfangle > (90 - angleThreshold) * 3.141592 / 180.0)) }
{
dbg << " : quad\n";
// Since the surface is lower than the threshold, change the effective
// prism growth vector to match with the surface vector, so that
// the Quad which is created lies on the original surface
//growthvectors.Elem(segpair_p1) = surfelem_vect;
// Add a quad element to account for the prism volume
// element which is going to be added
Element2d sel(QUAD);
sel.PNum(4) = mapto[seg_p1];
sel.PNum(3) = mapto[seg_p2];
sel.PNum(2) = segpair_p2;
sel.PNum(1) = segpair_p1;
sel.SetIndex(mesh[sej].si);
mesh.AddSurfaceElement(sel);
numquads++;
} }
else else
{ {
dbg << "\n"; // check if it doesn't contain the other edge as well
for (int k = 0; k < pnt1_elems.Size(); k++) // and if it doesn't contain both mark them as done and
// if necessary map them
for(SegmentIndex sej = 0; sej<nseg; sej++)
{ {
Element2d & pnt_sel = mesh.SurfaceElement(pnt1_elems[k]); if(mesh[sej][0] == mesh[sei][1] &&
if(pnt_sel.GetIndex() == mesh[sej].si) mesh[sej][1] == mesh[sei][0])
{ {
for(int l = 0; l < pnt_sel.GetNP(); l++) if(!blp.surfid.Contains(mesh[sej].si))
{ {
if(pnt_sel[l] == segpair_p1) segsel.Clear(sei);
pnt_sel[l] = mapto[seg_p1]; segsel.Clear(sej);
else if (pnt_sel[l] == segpair_p2) PointIndex mapped_point = PointIndex::INVALID;
pnt_sel[l] = mapto[seg_p2]; 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;
}
}
} }
} }
} }
for (int k = 0; k < pnt2_elems.Size(); k++) // add surface elements between layer and old domain
if(layer == blp.heights.Size())
{ {
Element2d & pnt_sel = mesh.SurfaceElement(pnt2_elems[k]); for(SurfaceElementIndex si = 0; si < nse; si++)
if(pnt_sel.GetIndex() == mesh[sej].si)
{ {
for(int l = 0; l < pnt_sel.GetNP(); l++) if(blp.surfid.Contains(mesh[si].GetIndex()))
{ {
if(pnt_sel[l] == segpair_p1) const auto& sel = mesh[si];
pnt_sel[l] = mapto.Get(seg_p1); Element2d newel = sel;
else if (pnt_sel[l] == segpair_p2) newel.SetIndex(map_surface_index(sel.GetIndex()));
pnt_sel[l] = mapto.Get(seg_p2); mesh.AddSurfaceElement(newel);
}
}
}
}
// }
}
else
{
// If the code comes here, it indicates that we are at
// a line segment pair which is at the intersection
// of two surfaces, both of which have to grow boundary
// layers.... here too, remapping the segments to the
// new points is required
mesh[sei][0] = mapto.Get(seg_p1);
mesh[sei][1] = mapto.Get(seg_p2);
mesh[sej][1] = mapto.Get(seg_p1);
mesh[sej][0] = mapto.Get(seg_p2);
}
}
} }
} }
} }
// Add prismatic cells at the boundaries // Add prismatic cells at the boundaries
cout << "Generating prism boundary layer volume elements...." << endl; 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(surfid.Contains(sel.GetIndex())) if(blp.surfid.Contains(sel.GetIndex()))
{ {
/*
Element el(PRISM);
for (int j = 0; j < sel.GetNP(); j++)
{
// Check (Doublecheck) if the corresponding point has a
// copy available for remapping
if (mapto.Get(sel[j]))
{
// Define the points of the newly added Prism cell
el[j+3] = mapto[sel[j]];
el[j] = sel[j];
}
else
{
el[j+3] = sel[j];
el[j] = sel[j];
}
}
el.SetIndex(1);
el.Invert();
mesh.AddVolumeElement(el);
numprisms++;
*/
// cout << "add element: " << endl;
int classify = 0; int classify = 0;
for (int j = 0; j < 3; j++) for(auto j : Range(sel.PNums()))
if (mapto[sel[j]]) if (mapto[sel[j]].IsValid())
classify += (1 << j); classify += (1 << j);
// cout << "classify = " << classify << endl; if(classify == 0)
continue;
Element el;
if(sel.GetType() == TRIG)
{
ELEMENT_TYPE types[] = { PRISM, TET, TET, PYRAMID, ELEMENT_TYPE types[] = { PRISM, TET, TET, PYRAMID,
TET, PYRAMID, PYRAMID, PRISM }; TET, PYRAMID, PYRAMID, PRISM };
int nums[] = { sel[0], sel[1], sel[2], mapto[sel[0]], mapto[sel[1]], mapto[sel[2]] }; int nums[] = { sel[0], sel[1], sel[2], mapto[sel[0]], mapto[sel[1]], mapto[sel[2]] };
@ -514,134 +500,120 @@ namespace netgen
{ 1, 2, 5, 4, 0, 0 }, { 1, 2, 5, 4, 0, 0 },
{ 0, 2, 1, 3, 5, 4 } { 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;
}
Element el(types[classify]); el = Element(types[classify]);
for (int i = 0; i < 6; i++) for(auto i : Range(el.PNums()))
el[i] = nums[vertices[classify][i]]; el.PNums()[i] = nums[vertices[classify][i]];
if(blp.new_matnrs.Size() > 0) }
el.SetIndex(blp.new_matnrs[layer-1]); 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 else
el.SetIndex(blp.new_matnr); {
// cout << "el = " << el << endl; throw Exception("This type of quad layer not yet implemented!");
if (classify != 0) }
}
el.SetIndex(blp.new_matnrs[layer-1]);
mesh.AddVolumeElement(el); mesh.AddVolumeElement(el);
} }
} }
// Finally switch the point indices of the surface elements // Finally switch the point indices of the surface elements
// to the newly added ones // to the newly added ones
cout << "Transferring boundary layer surface elements to new vertex references...." << endl; PrintMessage(3, "Transferring boundary layer surface elements to new vertex references...");
for (int i = 1; i <= nse; i++) for(SurfaceElementIndex sei : Range(nse))
{ {
Element2d & sel = mesh.SurfaceElement(i); auto& sel = mesh[sei];
if(surfid.Contains(sel.GetIndex())) bool to_move = blp.surfid.Contains(sel.GetIndex());
if(blp.domains.Size())
{ {
for (int j = 1; j <= sel.GetNP(); j++) 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())
// 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.Get(sel.PNum(j))) if(mapto[pnum].IsValid())
{
// Map the surface elements to the new points // Map the surface elements to the new points
sel.PNum(j) = mapto.Get(sel.PNum(j)); pnum = mapto[pnum];
} }
} }
}
} for(ElementIndex ei : Range(ne))
for (int i = 1; i <= ne; i++)
{
Element & el = mesh.VolumeElement(i);
if(el.GetIndex() != blp.bulk_matnr)
{
for (int j = 1; j <= el.GetNP(); j++)
{ {
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())
// 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.Get(el.PNum(j))) if(mapto[pnum].IsValid())
{ // Map the volume elements to the new points
// Map the surface elements to the new points pnum = mapto[pnum];
el.PNum(j) = mapto.Get(el.PNum(j));
} }
}
}
}
// 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
// for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++) // for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++)
for (PointIndex pi : mesh.Points().Range()) for (PointIndex pi = 1; pi <= np; pi++)
{
if(bndnodes.Test(pi)) mesh.AddLockedPoint(pi); if(bndnodes.Test(pi)) mesh.AddLockedPoint(pi);
}
// Now, actually pull back the old surface points to create // Now, actually pull back the old surface points to create
// the actual boundary layers // the actual boundary layers
cout << "Moving and optimising boundary layer points...." << endl; PrintMessage(3, "Moving and optimising boundary layer points...");
for (int i = 1; i <= np; i++) for (PointIndex i = 1; i <= np; i++)
{ {
NgArray<ElementIndex> vertelems;
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]);
if(layer == prismlayers)
{
mesh.Point(i).SetPoint(pointtomove + layerht * growthvectors.Elem(i));
meshtopo.GetVertexElements(i,vertelems);
for(int j = 1; j <= vertelems.Size(); j++)
{
// double sfact = 0.9;
Element volel = mesh.VolumeElement(vertelems.Elem(j));
if(((volel.GetType() == TET) || (volel.GetType() == TET10)) && (!volel.IsDeleted()))
{
//while((volel.Volume(mesh.Points()) <= 0.0) && (sfact >= 0.0))
//{
// mesh.Point(i).SetPoint(pointtomove + (sfact * layerht * growthvectors.Elem(i)));
// mesh.ImproveMesh();
// // Try to move the point back by one step but
// // if the volume drops to below zero, double back
// mesh.Point(i).SetPoint(pointtomove + ((sfact + 0.1) * layerht * growthvectors.Elem(i)));
// if(volel.Volume(mesh.Points()) <= 0.0)
// {
// mesh.Point(i).SetPoint(pointtomove + (sfact * layerht * growthvectors.Elem(i)));
// }
// sfact -= 0.1;
//}
volel.Delete();
}
}
}
else
{
mesh.Point(i).SetPoint(pointtomove + layerht * growthvectors.Elem(i));
}
} }
} }
mesh.Compress(); mesh.Compress();
} }
// Optimise the tet part of the volume mesh after all the modifications for(int i=1; i <= mesh.GetNFD(); i++)
// to the system are completed {
//OptimizeVolume(mparam,mesh); auto& fd = mesh.GetFaceDescriptor(i);
if(blp.surfid.Contains(fd.BCProperty()))
cout << "New NP: " << mesh.GetNP() << endl; {
cout << "Num of Quads: " << numquads << endl; if(blp.outside)
cout << "Num of Prisms: " << numprisms << endl; fd.SetDomainOut(blp.new_matnrs[blp.new_matnrs.Size()-1]);
cout << "Boundary Layer Generation....Done!" << endl; else
fd.SetDomainIn(blp.new_matnrs[blp.new_matnrs.Size()-1]);
dbg.close(); }
} }
PrintMessage(3, "New NP: ", mesh.GetNP());
PrintMessage(1, "Boundary Layer Generation....Done!");
}
} }

18
libsrc/meshing/boundarylayer.hpp
View File

@ -12,18 +12,16 @@ class BoundaryLayerParameters
{ {
public: public:
// parameters by Philippose .. // parameters by Philippose ..
NgArray<int> surfid; Array<int> surfid;
NgArray<double> heights; Array<double> heights;
NgArray<double> new_matnrs; Array<size_t> new_matnrs;
int prismlayers = 1; BitArray domains;
int bulk_matnr = 1; bool outside = false; // set the boundary layer on the outside
int new_matnr = 1; bool grow_edges = false;
double hfirst = 0.01;
double growthfactor = 1;
bool optimize = true;
}; };
DLL_HEADER extern void GenerateBoundaryLayer (Mesh & mesh, BoundaryLayerParameters & blp); DLL_HEADER void GenerateBoundaryLayer (Mesh & mesh,
const BoundaryLayerParameters & blp);
#endif #endif

162
libsrc/meshing/python_mesh.cpp
View File

@ -1,5 +1,7 @@
#ifdef NG_PYTHON #ifdef NG_PYTHON
#include <regex>
#include <../general/ngpython.hpp> #include <../general/ngpython.hpp>
#include <core/python_ngcore.hpp> #include <core/python_ngcore.hpp>
#include "python_mesh.hpp" #include "python_mesh.hpp"
@ -885,12 +887,13 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
meshingparameter_description.c_str(), meshingparameter_description.c_str(),
py::call_guard<py::gil_scoped_release>()) py::call_guard<py::gil_scoped_release>())
.def ("OptimizeVolumeMesh", [](Mesh & self) .def ("OptimizeVolumeMesh", [](Mesh & self, MeshingParameters* pars)
{ {
MeshingParameters mp; MeshingParameters mp;
mp.optsteps3d = 5; if(pars) mp = *pars;
else mp.optsteps3d = 5;
OptimizeVolume (mp, self); OptimizeVolume (mp, self);
},py::call_guard<py::gil_scoped_release>()) }, py::arg("mp"), py::call_guard<py::gil_scoped_release>())
.def ("OptimizeMesh2d", [](Mesh & self) .def ("OptimizeMesh2d", [](Mesh & self)
{ {
@ -929,63 +932,112 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
.def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard<py::gil_scoped_release>()) .def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard<py::gil_scoped_release>())
.def ("BoundaryLayer", FunctionPointer .def ("BoundaryLayer", [](Mesh & self, variant<string, int> boundary,
([](Mesh & self, int bc, py::list thicknesses, int volnr, py::list materials) variant<double, py::list> thickness,
{ variant<string, py::list> material,
int n = py::len(thicknesses); variant<string, int> domain, bool outside,
BoundaryLayerParameters blp; bool grow_edges)
for (int i = 1; i <= self.GetNFD(); i++)
if (self.GetFaceDescriptor(i).BCProperty() == bc)
blp.surfid.Append (i);
cout << "add layer at surfaces: " << blp.surfid << endl;
blp.prismlayers = n;
blp.growthfactor = 1.0;
// find max domain nr
int maxind = 0;
for (ElementIndex ei = 0; ei < self.GetNE(); ei++)
maxind = max (maxind, self[ei].GetIndex());
cout << "maxind = " << maxind << endl;
for ( int i=0; i<n; i++ )
{
blp.heights.Append( py::extract<double>(thicknesses[i])()) ;
blp.new_matnrs.Append( maxind+1+i );
self.SetMaterial (maxind+1+i, py::extract<string>(materials[i])().c_str());
}
blp.bulk_matnr = volnr;
GenerateBoundaryLayer (self, blp);
}
))
.def ("BoundaryLayer", FunctionPointer
([](Mesh & self, int bc, double thickness, int volnr, string material)
{ {
BoundaryLayerParameters blp; BoundaryLayerParameters blp;
if(int* bc = get_if<int>(&boundary); bc)
{
for (int i = 1; i <= self.GetNFD(); i++) for (int i = 1; i <= self.GetNFD(); i++)
if (self.GetFaceDescriptor(i).BCProperty() == bc) if(self.GetFaceDescriptor(i).BCProperty() == *bc)
blp.surfid.Append (i); blp.surfid.Append (i);
cout << "add layer at surfaces: " << blp.surfid << endl;
blp.prismlayers = 1;
blp.hfirst = thickness;
blp.growthfactor = 1.0;
// find max domain nr
int maxind = 0;
for (ElementIndex ei = 0; ei < self.GetNE(); ei++)
maxind = max (maxind, self[ei].GetIndex());
cout << "maxind = " << maxind << endl;
self.SetMaterial (maxind+1, material.c_str());
blp.new_matnr = maxind+1;
blp.bulk_matnr = volnr;
GenerateBoundaryLayer (self, blp);
} }
)) else
{
regex pattern(*get_if<string>(&boundary));
for(int i = 1; i<=self.GetNFD(); i++)
if(regex_match(self.GetFaceDescriptor(i).GetBCName(), pattern))
blp.surfid.Append(i);
}
if(double* pthickness = get_if<double>(&thickness); pthickness)
{
blp.heights.Append(*pthickness);
}
else
{
auto thicknesses = *get_if<py::list>(&thickness);
for(auto val : thicknesses)
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
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.grow_edges = grow_edges;
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,
R"delimiter(
Add boundary layer to mesh.
Parameters
----------
boundary : string or int
Boundary name or number.
thickness : float or List[float]
Thickness of boundary layer(s).
material : str or List[str]
Material name of boundary layer(s).
domain : str or int
Regexp for domain boundarylayer is going into.
outside : bool = False
If true add the layer on the outside
grow_edges : bool = False
Grow boundary layer over edges.
)delimiter")
.def ("EnableTable", [] (Mesh & self, string name, bool set) .def ("EnableTable", [] (Mesh & self, string name, bool set)
{ {

2
libsrc/meshing/topology.hpp
View File

@ -166,7 +166,7 @@ public:
{ return vert2element[vnr]; } { return vert2element[vnr]; }
void GetVertexSurfaceElements( int vnr, NgArray<SurfaceElementIndex>& elements ) const; void GetVertexSurfaceElements( int vnr, NgArray<SurfaceElementIndex>& elements ) const;
NgFlatArray<SurfaceElementIndex> GetVertexSurfaceElements (int vnr) const NgFlatArray<SurfaceElementIndex> GetVertexSurfaceElements(PointIndex vnr) const
{ return vert2surfelement[vnr]; } { return vert2surfelement[vnr]; }
NgFlatArray<SegmentIndex> GetVertexSegments (int vnr) const NgFlatArray<SegmentIndex> GetVertexSegments (int vnr) const

17
ng/ngpkg.cpp
View File

@ -1155,7 +1155,7 @@ namespace netgen
// Use an array to support creation of boundary // Use an array to support creation of boundary
// layers for multiple surfaces in the future... // layers for multiple surfaces in the future...
NgArray<int> surfid; Array<int> surfid;
int surfinp = 0; int surfinp = 0;
int prismlayers = 1; int prismlayers = 1;
double hfirst = 0.01; double hfirst = 0.01;
@ -1172,8 +1172,8 @@ namespace netgen
cout << "Number of surfaces entered = " << surfid.Size() << endl; cout << "Number of surfaces entered = " << surfid.Size() << endl;
cout << "Selected surfaces are:" << endl; cout << "Selected surfaces are:" << endl;
for(int i = 1; i <= surfid.Size(); i++) for(auto i : Range(surfid))
cout << "Surface " << i << ": " << surfid.Elem(i) << endl; cout << "Surface " << i << ": " << surfid[i] << endl;
cout << endl << "Enter number of prism layers: "; cout << endl << "Enter number of prism layers: ";
cin >> prismlayers; cin >> prismlayers;
@ -1189,9 +1189,14 @@ namespace netgen
BoundaryLayerParameters blp; BoundaryLayerParameters blp;
blp.surfid = surfid; blp.surfid = surfid;
blp.prismlayers = prismlayers; for(auto i : Range(prismlayers))
blp.hfirst = blp.hfirst; {
blp.growthfactor = growthfactor; auto layer = i+1;
if(growthfactor == 1)
blp.heights.Append(layer * hfirst);
else
blp.heights.Append(hfirst * (pow(growthfactor, (layer+1))-1)/(growthfactor-1));
}
GenerateBoundaryLayer (*mesh, blp); GenerateBoundaryLayer (*mesh, blp);
return TCL_OK; return TCL_OK;
} }

31
tests/pytest/test_boundarylayer.py Normal file
View File

@ -0,0 +1,31 @@
import pytest
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])
def test_boundarylayer(outside, capfd):
mesh = unit_cube.GenerateMesh(maxh=0.3)
ne_before = mesh.ne
layer_surfacenames = ["right", "top"]
mesh.BoundaryLayer("|".join(layer_surfacenames), [0.01, 0.02], "layer", outside=outside, grow_edges=True)
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