ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-24 21:10:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

modernize and improve GenerateBoundaryLayer

Browse Source
This commit is contained in:
Christopher Lackner 2020-04-19 20:00:06 +02:00
parent 83a48af36a
commit 58e6e5dc18
6 changed files with 460 additions and 539 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

819
libsrc/meshing/boundarylayer.cpp
View File

@ -103,24 +103,16 @@ namespace netgen
function, in order to calculate the effective direction
in which the prismatic layer should grow
*/
void GetSurfaceNormal(Mesh & mesh, const Element2d & el, int Vertex, Vec3d & SurfaceNormal)
Vec<3> GetSurfaceNormal(Mesh & mesh, const Element2d & el)
{
int Vertex_A;
int Vertex_B;
Vertex_A = Vertex + 1;
if(Vertex_A > el.GetNP()) Vertex_A = 1;
Vertex_B = Vertex - 1;
if(Vertex_B <= 0) Vertex_B = el.GetNP();
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();
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;
}
@ -141,507 +133,390 @@ namespace netgen
Currently, the layer height is calculated using:
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");
// Angle between a surface element and a growth-vector below which
// Angle between a surface element and a growth-vector below which
// a prism is project onto that surface as a quad
// (in degrees)
double angleThreshold = 5.0;
NgArray<int> surfid (blp.surfid);
int prismlayers = blp.prismlayers;
double hfirst = blp.hfirst;
double growthfactor = blp.growthfactor;
NgArray<double> heights (blp.heights);
bool grow_edges = false; // grow layer at edges
// Monitor and print out the number of prism and quad elements
// 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--)
PrintMessage(1, "Generating boundary layer...");
PrintMessage(3, "Old NP: ", mesh.GetNP());
PrintMessage(3, "Old NSE: ",mesh.GetNSE());
for(int layer = blp.heights.Size(); layer >= 1; layer--)
{
cout << "Generating layer: " << layer << endl;
const MeshTopology& meshtopo = mesh.GetTopology();
const_cast<MeshTopology &> (meshtopo).SetBuildEdges(true);
const_cast<MeshTopology &> (meshtopo).SetBuildFaces(true);
const_cast<MeshTopology &> (meshtopo).Update();
PrintMessage(3, "Generating layer: ", layer);
double layerht = hfirst;
if(heights.Size()>0)
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 nseg = mesh.GetNSeg();
// Indicate which points need to be remapped
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<Array<Vec<3>>, PointIndex> all_growthvectors(np);
// Bit array to identify all the points belonging
// to the surface of interest
bndnodes.Clear();
// 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 normal = GetSurfaceNormal(mesh,sel);
if(!blp.outside)
normal *= -1;
for(int j : Range(sel.PNums()))
{
// Set the bitarray to indicate that the
// point is part of the required set
bndnodes.SetBit(sel[j]);
// Add the surface normal to the already existent one
// (This gives the effective normal direction at corners
// and curved areas)
all_growthvectors[sel[j]].Append(normal);
}
}
if (!blp.grow_edges)
for(const auto& sel : mesh.SurfaceElements())
{
bndnodes.Clear(sel[0]);
bndnodes.Clear(sel[1]);
}
// 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...");
for (PointIndex pi = 1; pi <= np; pi++)
{
layerht = heights[layer-1];
}
else
{
if(growthfactor == 1)
if (bndnodes.Test(pi))
{
layerht = layer * hfirst;
mapto[pi] = mesh.AddPoint(mesh[pi]);
growthvectors[pi] = all_growthvectors[pi][0];
for(int i = 1; i < all_growthvectors[pi].Size(); i++)
{
auto& veci = all_growthvectors[pi][i];
for(auto j : Range(i))
{
auto& vecj = all_growthvectors[pi][j];
veci -= 1./(vecj.Length()+1e-10) * (veci * vecj) * vecj;
}
growthvectors[pi] += veci;
}
// growthvectors[pi].Normalize();
// growthvectors[pi] *= -1.0;
}
else
{
layerht = hfirst*(pow(growthfactor,(layer+1)) - 1)/(growthfactor - 1);
mapto[pi].Invalidate();
growthvectors[pi] = {0,0,0};
}
}
cout << "Layer Height = " << layerht << endl;
// 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 nseg = mesh.GetNSeg();
// Indicate which points need to be remapped
NgBitArray bndnodes(np+1); // big enough for 1-based array
// Map of the old points to the new points
NgArray<PointIndex, PointIndex::BASE> mapto(np);
// Growth vectors for the prismatic layer based on
// the effective surface normal at a given point
NgArray<Vec3d, PointIndex::BASE> growthvectors(np);
// Bit array to identify all the points belonging
// to the surface of interest
bndnodes.Clear();
// 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
cout << "Marking points for remapping...." << endl;
for (SurfaceElementIndex si = 0; si < nse; si++)
if (surfid.Contains(mesh[si].GetIndex()))
{
const Element2d & sel = mesh[si];
for(int j = 0; j < sel.GetNP(); j++)
{
// Set the bitarray to indicate that the
// point is part of the required set
bndnodes.Set(sel[j]);
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
// (This gives the effective normal direction at corners
// and curved areas)
growthvectors[sel[j]] += surfacenormal;
}
}
if (!grow_edges)
for (SegmentIndex sei = 0; sei <= nseg; sei++)
{
bndnodes.Clear (mesh[sei][0]);
bndnodes.Clear (mesh[sei][1]);
}
// 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
cout << "Cloning points and calculating growth vectors...." << endl;
for (PointIndex pi = 1; pi <= np; pi++)
{
if (bndnodes.Test(pi))
{
mapto[pi] = mesh.AddPoint (mesh[pi]);
growthvectors[pi].Normalize();
growthvectors[pi] *= -1.0;
}
else
{
mapto[pi] = 0;
growthvectors[pi] = Vec3d(0,0,0);
}
}
// Add quad surface elements at edges for surfaces which
// don't have boundary layers
// Add quad surface elements at edges for surfaces which
// don't have boundary layers
// Bit array to keep track of segments already processed
NgBitArray segsel(nseg);
// Bit array to keep track of segments already processed
BitArray segsel(nseg);
// Set them all to "1" to initially activate all segments
segsel.Set();
// Set them all to "1" to initially activate all segments
segsel.Set();
cout << "Adding 2D Quad elements on required surfaces...." << endl;
PrintMessage(3, "Adding 2D Quad elements on required surfaces...");
if (grow_edges)
for (SegmentIndex sei = 0; sei <= nseg; sei++)
{
PointIndex seg_p1 = mesh[sei][0];
PointIndex seg_p2 = mesh[sei][1];
// 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) && surfid.Contains(mesh[sei].si))
{
// 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++)
{
PointIndex segpair_p1 = mesh[sej][1];
PointIndex segpair_p2 = mesh[sej][0];
// 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) && ((segpair_p1 == seg_p1) && (segpair_p2 == seg_p2)))
{
// clear bit to indicate that processing of this segment is done
segsel.Clear(sej);
// Only worry about those surfaces which are not in the
// boundary layer list
if(!surfid.Contains(mesh[sej].si))
{
SurfaceElementIndex pnt_commelem = 0;
NgArray<SurfaceElementIndex> pnt1_elems;
NgArray<SurfaceElementIndex> pnt2_elems;
meshtopo.GetVertexSurfaceElements(segpair_p1,pnt1_elems);
meshtopo.GetVertexSurfaceElements(segpair_p2,pnt2_elems);
if(blp.grow_edges)
for(SegmentIndex sei = 0; sei < nseg; sei++)
{
PointIndex seg_p1 = mesh[sei][0];
PointIndex seg_p2 = mesh[sei][1];
for(int k = 0; k < pnt1_elems.Size(); k++)
{
const Element2d & pnt1_sel = mesh.SurfaceElement(pnt1_elems[k]);
for(int l = 0; l < pnt2_elems.Size(); l++)
{
const Element2d & pnt2_sel = mesh.SurfaceElement(pnt2_elems[l]);
if((pnt1_sel.GetIndex() == mesh[sej].si)
&& (pnt2_sel.GetIndex() == mesh[sej].si)
&& (pnt1_elems[k] == pnt2_elems[l]))
{
pnt_commelem = pnt1_elems[k];
}
}
}
// 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) && blp.surfid.Contains(mesh[sei].si))
{
// clear the bit to indicate that this segment has been processed
segsel.Clear(sei);
/*
int pnum_commelem = 0;
for(int k = 1; k <= mesh.SurfaceElement(pnt_commelem).GetNP(); k++)
{
if((mesh.SurfaceElement(pnt_commelem).PNum(k) != segpair_p1)
&& (mesh.SurfaceElement(pnt_commelem).PNum(k) != segpair_p2))
{
pnum_commelem = mesh.SurfaceElement(pnt_commelem).PNum(k);
}
}
*/
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
mesh[sei][0] = mapto[seg_p1];
mesh[sei][1] = mapto[seg_p2];
mesh[sej][1] = mapto[seg_p1];
mesh[sej][0] = mapto[seg_p2];
if((surfangle < (90 + angleThreshold) * 3.141592 / 180.0)
&& (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
{
dbg << "\n";
for (int k = 0; k < pnt1_elems.Size(); k++)
{
Element2d & pnt_sel = mesh.SurfaceElement(pnt1_elems[k]);
if(pnt_sel.GetIndex() == mesh[sej].si)
{
for(int l = 0; l < pnt_sel.GetNP(); l++)
{
if(pnt_sel[l] == segpair_p1)
pnt_sel[l] = mapto[seg_p1];
else if (pnt_sel[l] == segpair_p2)
pnt_sel[l] = mapto[seg_p2];
}
}
}
for (int k = 0; k < pnt2_elems.Size(); k++)
{
Element2d & pnt_sel = mesh.SurfaceElement(pnt2_elems[k]);
if(pnt_sel.GetIndex() == mesh[sej].si)
{
for(int l = 0; l < pnt_sel.GetNP(); l++)
{
if(pnt_sel[l] == segpair_p1)
pnt_sel[l] = mapto.Get(seg_p1);
else if (pnt_sel[l] == segpair_p2)
pnt_sel[l] = mapto.Get(seg_p2);
}
}
}
}
// }
}
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
cout << "Generating prism boundary layer volume elements...." << endl;
// Find matching segment pair on other surface
for (SegmentIndex sej = 0; sej < nseg; sej++)
{
PointIndex segpair_p1 = mesh[sej][1];
PointIndex segpair_p2 = mesh[sej][0];
for (SurfaceElementIndex si = 0; si < nse; si++)
{
Element2d & sel = mesh.SurfaceElement(si);
if(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;
for (int j = 0; j < 3; j++)
if (mapto[sel[j]])
classify += (1 << j);
// 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) && ((segpair_p1 == seg_p1) && (segpair_p2 == seg_p2)))
{
// clear bit to indicate that processing of this segment is done
segsel.Clear(sej);
// cout << "classify = " << classify << endl;
// Only worry about those surfaces which are not in the
// boundary layer list
if(!blp.surfid.Contains(mesh[sej].si))
{
SurfaceElementIndex pnt_commelem;
SetInvalid(pnt_commelem);
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] =
{
auto pnt1_elems = meshtopo.GetVertexSurfaceElements(segpair_p1);
auto pnt2_elems = meshtopo.GetVertexSurfaceElements(segpair_p2);
for(auto pnt1_sei : pnt1_elems)
{
const auto& pnt1_sel = mesh[pnt1_sei];
for(auto pnt2_sei : pnt2_elems)
{
const Element2d & pnt2_sel = mesh.SurfaceElement(pnt2_sei);
if((pnt1_sel.GetIndex() == mesh[sej].si)
&& (pnt2_sel.GetIndex() == mesh[sej].si)
&& (pnt1_sei == pnt2_sei))
{
pnt_commelem = pnt1_sei;
}
}
}
const Element2d & commsel = mesh.SurfaceElement(pnt_commelem);
auto surfelem_vect = GetSurfaceNormal(mesh, commsel);
if(blp.outside)
surfelem_vect *= -1;
double surfangle = Angle(growthvectors[segpair_p1],surfelem_vect);
// remap the segments to the new points
if(!blp.outside)
{
mesh[sei][0] = mapto[seg_p1];
mesh[sei][1] = mapto[seg_p2];
mesh[sej][1] = mapto[seg_p1];
mesh[sej][0] = mapto[seg_p2];
}
if((surfangle < (90 + angleThreshold) * 3.141592 / 180.0)
&& (surfangle > (90 - angleThreshold) * 3.141592 / 180.0))
{
// 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);
if(blp.outside)
Swap(seg_p1, seg_p2);
sel.PNum(4) = mapto[seg_p1];
sel.PNum(3) = mapto[seg_p2];
sel.PNum(2) = seg_p2;
sel.PNum(1) = seg_p1;
sel.SetIndex(mesh[sej].si);
mesh.AddSurfaceElement(sel);
numquads++;
}
else
{
for (int k = 0; k < pnt1_elems.Size(); k++)
{
Element2d & pnt_sel = mesh.SurfaceElement(pnt1_elems[k]);
if(pnt_sel.GetIndex() == mesh[sej].si)
{
for(int l = 0; l < pnt_sel.GetNP(); l++)
{
if(pnt_sel[l] == segpair_p1)
pnt_sel[l] = mapto[seg_p1];
else if (pnt_sel[l] == segpair_p2)
pnt_sel[l] = mapto[seg_p2];
}
}
}
for (int k = 0; k < pnt2_elems.Size(); k++)
{
Element2d & pnt_sel = mesh.SurfaceElement(pnt2_elems[k]);
if(pnt_sel.GetIndex() == mesh[sej].si)
{
for(int l = 0; l < pnt_sel.GetNP(); l++)
{
if(pnt_sel[l] == segpair_p1)
pnt_sel[l] = mapto[seg_p1];
else if (pnt_sel[l] == segpair_p2)
pnt_sel[l] = mapto[seg_p2];
}
}
}
}
// }
}
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[seg_p1];
mesh[sei][1] = mapto[seg_p2];
mesh[sej][1] = mapto[seg_p1];
mesh[sej][0] = mapto[seg_p2];
}
}
}
}
}
// Add prismatic cells at the boundaries
PrintMessage(3, "Generating prism boundary layer volume elements...");
for (SurfaceElementIndex si = 0; si < nse; si++)
{
Element2d & sel = mesh.SurfaceElement(si);
if(blp.surfid.Contains(sel.GetIndex()))
{
int classify = 0;
for (int j = 0; j < 3; j++)
if (mapto[sel[j]].IsValid())
classify += (1 << j);
// cout << "classify = " << classify << endl;
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 },
{ 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;
}
Element el(types[classify]);
for (int i = 0; i < 6; i++)
el[i] = nums[vertices[classify][i]];
if(blp.new_matnrs.Size() > 0)
el.SetIndex(blp.new_matnrs[layer-1]);
else
el.SetIndex(blp.new_matnr);
// cout << "el = " << el << endl;
if (classify != 0)
mesh.AddVolumeElement(el);
}
}
// Finally switch the point indices of the surface elements
// to the newly added ones
cout << "Transferring boundary layer surface elements to new vertex references...." << endl;
for (int i = 1; i <= nse; i++)
{
Element2d & sel = mesh.SurfaceElement(i);
if(surfid.Contains(sel.GetIndex()))
{
for (int j = 1; j <= sel.GetNP(); j++)
{
// Check (Doublecheck) if the corresponding point has a
Element el(types[classify]);
for (int i = 0; i < 6; i++)
el[i] = nums[vertices[classify][i]];
el.SetIndex(blp.new_matnrs[layer-1]);
if (classify != 0)
mesh.AddVolumeElement(el);
}
}
// 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))
{
auto& sel = mesh[sei];
if((blp.outside && !blp.surfid.Contains(sel.GetIndex())) ||
(!blp.outside && blp.surfid.Contains(sel.GetIndex())))
{
for(auto& pnum : sel.PNums())
// Check (Doublecheck) if the corresponding point has a
// copy available for remapping
if (mapto.Get(sel.PNum(j)))
{
// Map the surface elements to the new points
sel.PNum(j) = mapto.Get(sel.PNum(j));
}
}
}
}
for (int i = 1; i <= ne; i++)
{
Element & el = mesh.VolumeElement(i);
if(el.GetIndex() != blp.bulk_matnr)
if(mapto[pnum].IsValid())
// Map the surface elements to the new points
pnum = mapto[pnum];
}
}
if(blp.outside)
for(ElementIndex ei : Range(ne))
{
for (int j = 1; j <= el.GetNP(); j++)
{
// Check (Doublecheck) if the corresponding point has a
// copy available for remapping
if (mapto.Get(el.PNum(j)))
{
// Map the surface elements to the new points
el.PNum(j) = mapto.Get(el.PNum(j));
}
}
auto& el = mesh[ei];
for(auto& pnum : el.PNums())
// Check (Doublecheck) if the corresponding point has a
// copy available for remapping
if(mapto[pnum].IsValid())
// Map the surface elements to the new points
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...");
// 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 : mesh.Points().Range())
{
if(bndnodes.Test(pi)) mesh.AddLockedPoint(pi);
}
for (PointIndex i = 1; i <= np; i++)
{
if(bndnodes.Test(i))
{
MeshPoint pointtomove;
// Now, actually pull back the old surface points to create
// the actual boundary layers
cout << "Moving and optimising boundary layer points...." << endl;
for (int i = 1; i <= np; i++)
{
NgArray<ElementIndex> vertelems;
pointtomove = mesh.Point(i);
mesh.Point(i).SetPoint(pointtomove + layerht * growthvectors[i]);
}
}
mesh.Compress();
}
if(bndnodes.Test(i))
{
MeshPoint pointtomove;
pointtomove = mesh.Point(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();
}
// Optimise the tet part of the volume mesh after all the modifications
// to the system are completed
//OptimizeVolume(mparam,mesh);
for(int i=1; i <= mesh.GetNFD(); i++)
{
auto& fd = mesh.GetFaceDescriptor(i);
if(blp.surfid.Contains(fd.SurfNr()))
{
if(blp.outside)
fd.SetDomainOut(blp.new_matnrs[0]);
else
fd.SetDomainIn(blp.new_matnrs[0]);
}
}
cout << "New NP: " << mesh.GetNP() << endl;
cout << "Num of Quads: " << numquads << endl;
cout << "Num of Prisms: " << numprisms << endl;
cout << "Boundary Layer Generation....Done!" << endl;
dbg.close();
PrintMessage(3, "New NP: ", mesh.GetNP());
PrintMessage(3, "Num of Quads: ", numquads);
PrintMessage(3, "Num of Prisms: ", numprisms);
PrintMessage(1, "Boundary Layer Generation....Done!");
}
}

17
libsrc/meshing/boundarylayer.hpp
View File

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

127
libsrc/meshing/python_mesh.cpp
View File

@ -1,5 +1,7 @@
#ifdef NG_PYTHON
#include <regex>
#include <../general/ngpython.hpp>
#include <core/python_ngcore.hpp>
#include "python_mesh.hpp"
@ -885,12 +887,13 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
meshingparameter_description.c_str(),
py::call_guard<py::gil_scoped_release>())
.def ("OptimizeVolumeMesh", [](Mesh & self)
.def ("OptimizeVolumeMesh", [](Mesh & self, MeshingParameters* pars)
{
MeshingParameters mp;
mp.optsteps3d = 5;
if(pars) mp = *pars;
else mp.optsteps3d = 5;
OptimizeVolume (mp, self);
},py::call_guard<py::gil_scoped_release>())
}, py::arg("mp"), py::call_guard<py::gil_scoped_release>())
.def ("OptimizeMesh2d", [](Mesh & self)
{
@ -929,63 +932,87 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
.def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard<py::gil_scoped_release>())
.def ("BoundaryLayer", FunctionPointer
([](Mesh & self, int bc, py::list thicknesses, int volnr, py::list materials)
.def ("BoundaryLayer", [](Mesh & self, variant<string, int> boundary,
variant<double, py::list> thickness,
variant<string, py::list> material,
variant<string, int> domain, bool outside,
bool grow_edges)
{
int n = py::len(thicknesses);
BoundaryLayerParameters blp;
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++ )
if(int* bc = std::get_if<int>(&boundary); bc)
{
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());
for (int i = 1; i <= self.GetNFD(); i++)
if(self.GetFaceDescriptor(i).BCProperty() == *bc)
blp.surfid.Append (i);
}
blp.bulk_matnr = volnr;
else
{
regex pattern(std::get<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<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;
for(auto i : Range(prismlayers))
blp.new_matnrs.Append(first_new_mat + i);
if(string* pmaterial = get_if<string>(&material); pmaterial)
{
for(auto i : Range(prismlayers))
self.SetMaterial(first_new_mat + i, *pmaterial);
}
else
{
auto materials = get<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.outside = outside;
blp.grow_edges = grow_edges;
GenerateBoundaryLayer (self, blp);
}
))
}, py::arg("boundary"), py::arg("thickness"), py::arg("material"),
py::arg("domain") = 1, py::arg("outside") = false,
py::arg("grow_edges") = false, R"delimiter(
Add boundary layer to mesh.
.def ("BoundaryLayer", FunctionPointer
([](Mesh & self, int bc, double thickness, int volnr, string material)
{
BoundaryLayerParameters blp;
Parameters
----------
for (int i = 1; i <= self.GetNFD(); i++)
if (self.GetFaceDescriptor(i).BCProperty() == bc)
blp.surfid.Append (i);
boundary : string or int
Boundary name or number.
cout << "add layer at surfaces: " << blp.surfid << endl;
thickness : float or List[float]
Thickness of boundary layer(s).
blp.prismlayers = 1;
blp.hfirst = thickness;
blp.growthfactor = 1.0;
material : str or List[str]
Material name of boundary layer(s).
// 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);
}
))
domain : string or int
Add layer into domain specified by name or number.
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)
{

2
libsrc/meshing/topology.hpp
View File

@ -166,7 +166,7 @@ public:
{ return vert2element[vnr]; }
void GetVertexSurfaceElements( int vnr, NgArray<SurfaceElementIndex>& elements ) const;
NgFlatArray<SurfaceElementIndex> GetVertexSurfaceElements (int vnr) const
NgFlatArray<SurfaceElementIndex> GetVertexSurfaceElements(PointIndex vnr) const
{ return vert2surfelement[vnr]; }
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
// layers for multiple surfaces in the future...
NgArray<int> surfid;
Array<int> surfid;
int surfinp = 0;
int prismlayers = 1;
double hfirst = 0.01;
@ -1172,8 +1172,8 @@ namespace netgen
cout << "Number of surfaces entered = " << surfid.Size() << endl;
cout << "Selected surfaces are:" << endl;
for(int i = 1; i <= surfid.Size(); i++)
cout << "Surface " << i << ": " << surfid.Elem(i) << endl;
for(auto i : Range(surfid))
cout << "Surface " << i << ": " << surfid[i] << endl;
cout << endl << "Enter number of prism layers: ";
cin >> prismlayers;
@ -1189,9 +1189,14 @@ namespace netgen
BoundaryLayerParameters blp;
blp.surfid = surfid;
blp.prismlayers = prismlayers;
blp.hfirst = blp.hfirst;
blp.growthfactor = growthfactor;
for(auto i : Range(prismlayers))
{
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);
return TCL_OK;
}

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

@ -0,0 +1,17 @@
import pytest
from netgen.csg import *
@pytest.mark.parametrize("outside", [True, False])
def test_boundarylayer(outside):
mesh = unit_cube.GenerateMesh(maxh=0.3)
ne_before = mesh.ne
nse_in_layer = 0
layer_surfacenames = ["right", "top"]
for el in mesh.Elements2D():
if mesh.GetBCName(el.index-1) in layer_surfacenames:
nse_in_layer += 1
mesh.BoundaryLayer("|".join(layer_surfacenames), [0.01, 0.02], "layer", outside=outside, grow_edges=True)
assert mesh.ne == ne_before + 2 * nse_in_layer