From 58e6e5dc182971efba7183cfed5ee49f224c382f Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Sun, 19 Apr 2020 20:00:06 +0200 Subject: [PATCH] modernize and improve GenerateBoundaryLayer --- libsrc/meshing/boundarylayer.cpp | 819 ++++++++++++----------------- libsrc/meshing/boundarylayer.hpp | 17 +- libsrc/meshing/python_mesh.cpp | 127 +++-- libsrc/meshing/topology.hpp | 2 +- ng/ngpkg.cpp | 17 +- tests/pytest/test_boundarylayer.py | 17 + 6 files changed, 460 insertions(+), 539 deletions(-) create mode 100644 tests/pytest/test_boundarylayer.py diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index 4d132615..0d61d4cc 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -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 surfid (blp.surfid); - int prismlayers = blp.prismlayers; - double hfirst = blp.hfirst; - double growthfactor = blp.growthfactor; - NgArray 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 (meshtopo).SetBuildEdges(true); - const_cast (meshtopo).SetBuildFaces(true); - const_cast (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 mapto(np); + + // Growth vectors for the prismatic layer based on + // the effective surface normal at a given point + Array, PointIndex> growthvectors(np); + Array>, 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 mapto(np); - - // Growth vectors for the prismatic layer based on - // the effective surface normal at a given point - NgArray 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 pnt1_elems; - NgArray 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 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!"); } - } - diff --git a/libsrc/meshing/boundarylayer.hpp b/libsrc/meshing/boundarylayer.hpp index a174a55a..4c7d62a8 100644 --- a/libsrc/meshing/boundarylayer.hpp +++ b/libsrc/meshing/boundarylayer.hpp @@ -12,18 +12,15 @@ class BoundaryLayerParameters { public: // parameters by Philippose .. - NgArray surfid; - NgArray heights; - NgArray new_matnrs; - int prismlayers = 1; - int bulk_matnr = 1; - int new_matnr = 1; - double hfirst = 0.01; - double growthfactor = 1; - bool optimize = true; + Array surfid; + Array heights; + Array 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 diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 7e1f4821..05e0edbb 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -1,5 +1,7 @@ #ifdef NG_PYTHON +#include + #include <../general/ngpython.hpp> #include #include "python_mesh.hpp" @@ -885,12 +887,13 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) meshingparameter_description.c_str(), py::call_guard()) - .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::arg("mp"), py::call_guard()) .def ("OptimizeMesh2d", [](Mesh & self) { @@ -929,63 +932,87 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard()) - .def ("BoundaryLayer", FunctionPointer - ([](Mesh & self, int bc, py::list thicknesses, int volnr, py::list materials) + .def ("BoundaryLayer", [](Mesh & self, variant boundary, + variant thickness, + variant material, + variant 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(&boundary); bc) { - blp.heights.Append( py::extract(thicknesses[i])()) ; - blp.new_matnrs.Append( maxind+1+i ); - self.SetMaterial (maxind+1+i, py::extract(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(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(&thickness); pthickness) + { + blp.heights.Append(*pthickness); + } + else + { + auto thicknesses = get(thickness); + for(auto val : thicknesses) + blp.heights.Append(val.cast()); + } + + 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(&material); pmaterial) + { + for(auto i : Range(prismlayers)) + self.SetMaterial(first_new_mat + i, *pmaterial); + } + else + { + auto materials = get(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()); + } + + 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) { diff --git a/libsrc/meshing/topology.hpp b/libsrc/meshing/topology.hpp index 558646a2..8335f101 100644 --- a/libsrc/meshing/topology.hpp +++ b/libsrc/meshing/topology.hpp @@ -166,7 +166,7 @@ public: { return vert2element[vnr]; } void GetVertexSurfaceElements( int vnr, NgArray& elements ) const; - NgFlatArray GetVertexSurfaceElements (int vnr) const + NgFlatArray GetVertexSurfaceElements(PointIndex vnr) const { return vert2surfelement[vnr]; } NgFlatArray GetVertexSegments (int vnr) const diff --git a/ng/ngpkg.cpp b/ng/ngpkg.cpp index afd4d0e1..52279b19 100644 --- a/ng/ngpkg.cpp +++ b/ng/ngpkg.cpp @@ -1155,7 +1155,7 @@ namespace netgen // Use an array to support creation of boundary // layers for multiple surfaces in the future... - NgArray surfid; + Array 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; } diff --git a/tests/pytest/test_boundarylayer.py b/tests/pytest/test_boundarylayer.py new file mode 100644 index 00000000..9be268c0 --- /dev/null +++ b/tests/pytest/test_boundarylayer.py @@ -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 + +