diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index 4d132615..ae7c7a30 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -103,32 +103,21 @@ 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) + inline 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; } - - - - /* Philippose Rajan - 11 June 2009 + modified by Christopher Lackner Apr 2020 Added an initial experimental function for generating prismatic boundary layers on @@ -141,507 +130,490 @@ 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"); + 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 - // a prism is project onto that surface as a quad - // (in degrees) - double angleThreshold = 5.0; + map, int> domains_to_surf_index; + map, int> pi_to_edgenr; + map last_layer_surface_index_map; + int max_surface_index = mesh.GetNFD(); - NgArray surfid (blp.surfid); - int prismlayers = blp.prismlayers; - double hfirst = blp.hfirst; - double growthfactor = blp.growthfactor; - NgArray heights (blp.heights); + int max_edge_nr = -1; + for(const auto& seg : mesh.LineSegments()) + if(seg.edgenr > max_edge_nr) + max_edge_nr = seg.edgenr; - bool grow_edges = false; // grow layer at edges - - - // 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--) + 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) + auto map_surface_index = [&](auto si) { - layerht = heights[layer-1]; - } - else - { - if(growthfactor == 1) + if(last_layer_surface_index_map.find(si) == last_layer_surface_index_map.end()) { - layerht = layer * hfirst; + last_layer_surface_index_map[si] = ++max_surface_index; + auto& old_fd = mesh.GetFaceDescriptor(si); + int domout = blp.outside ? old_fd.DomainOut() : blp.new_matnrs[layer-1]; + int domin = blp.outside ? blp.new_matnrs[layer-1] : old_fd.DomainIn(); + FaceDescriptor fd(max_surface_index-1, + domin, domout, -1); + fd.SetBCProperty(max_surface_index); + mesh.AddFaceDescriptor(fd); + return max_surface_index; } + return last_layer_surface_index_map[si]; + }; + + 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); + growthvectors = 0.; + + // 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 n2 = GetSurfaceNormal(mesh,sel); + if(!blp.outside) + n2 *= -1; + for(auto pi : sel.PNums()) + { + // Set the bitarray to indicate that the + // point is part of the required set + bndnodes.SetBit(pi); + + // Add the surface normal to the already existent one + // (This gives the effective normal direction at corners + // and curved areas) + + auto& n1 = growthvectors[pi]; + if(n1.Length() == 0) { n1 = n2; continue; } + auto n1n2 = n1 * n2; + auto n1n1 = n1 * n1; + auto n2n2 = n2 * n2; + if(n2n2 - n1n2*n1n2/n1n1 == 0) { n1 = n2; continue; } + n1 += (n2n2 - n1n2)/(n2n2 - n1n2*n1n2/n1n1) * (n2 - n1n2/n1n1 * n1); + } + } + + if (!blp.grow_edges) + for(const auto& sel : mesh.LineSegments()) + { + 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++) + { + if (bndnodes.Test(pi)) + mapto[pi] = mesh.AddPoint(mesh[pi]); else - { - layerht = hfirst*(pow(growthfactor,(layer+1)) - 1)/(growthfactor - 1); - } + mapto[pi].Invalidate(); } - - 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(); + // Add quad surface elements at edges for surfaces which + // don't have boundary layers - // Indicate which points need to be remapped - NgBitArray bndnodes(np+1); // big enough for 1-based array + // Bit array to keep track of segments already processed + BitArray segsel(nseg); - // Map of the old points to the new points - NgArray mapto(np); + // Set them all to "1" to initially activate all segments + segsel.Set(); - // Growth vectors for the prismatic layer based on - // the effective surface normal at a given point - NgArray growthvectors(np); + PrintMessage(3, "Adding 2D Quad elements on required surfaces..."); - // 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 - - // Bit array to keep track of segments already processed - NgBitArray segsel(nseg); - - // Set them all to "1" to initially activate all segments - segsel.Set(); - - cout << "Adding 2D Quad elements on required surfaces...." << endl; - - 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); - - 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]; - } - } - } - - /* - 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; - - 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); - - // 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 }, - { 1, 2, 5, 4, 0, 0 }, - { 0, 2, 1, 3, 5, 4 } - }; - - 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())) + if(blp.grow_edges) + for(SegmentIndex sei = 0; sei < nseg; sei++) { - for (int j = 1; j <= sel.GetNP(); j++) + 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)) { - // 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) - { - for (int j = 1; j <= el.GetNP(); j++) + if(blp.surfid.Contains(mesh[sei].si)) { - // Check (Doublecheck) if the corresponding point has a - // copy available for remapping - if (mapto.Get(el.PNum(j))) + // 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++) { - // Map the surface elements to the new points - el.PNum(j) = mapto.Get(el.PNum(j)); - } - } - } - } + PointIndex segpair_p1 = mesh[sej][1]; + PointIndex segpair_p2 = mesh[sej][0]; - - - - // 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); - } - - // 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; - - 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())) + // 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))) { - //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(); + // 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(!blp.surfid.Contains(mesh[sej].si)) + { + SurfaceElementIndex pnt_commelem; + SetInvalid(pnt_commelem); + + auto pnt1_elems = meshtopo.GetVertexSurfaceElements(segpair_p1); + auto pnt2_elems = meshtopo.GetVertexSurfaceElements(segpair_p2); + + for(auto pnt1_sei : pnt1_elems) + if(mesh[pnt1_sei].GetIndex() == mesh[sej].si) + for(auto pnt2_sei : pnt2_elems) + if(pnt1_sei == pnt2_sei) + pnt_commelem = pnt1_sei; + + if(IsInvalid(pnt_commelem)) + throw Exception("Couldn't find element on other side for " + ToString(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()) + { + domains_to_surf_index[domains] = ++max_surface_index; + domains_to_surf_index[make_tuple(max_surface_index, get<1>(domains), get<2>(domains))] = max_surface_index; + FaceDescriptor fd(max_surface_index-1, + get<1>(domains), + get<2>(domains), + -1); + fd.SetBCProperty(max_surface_index); + mesh.AddFaceDescriptor(fd); + mesh.SetBCName(max_surface_index-1, + mesh.GetBCName(get<0>(domains)-1)); + } + auto new_index = domains_to_surf_index[domains]; + 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 + if(layer == blp.heights.Size()) + { + Segment s1 = mesh[sei]; + 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)) + { + 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); + } + } + + // remap the segments to the new points + mesh[sei][0] = mapto[mesh[sei][0]]; + mesh[sei][1] = mapto[mesh[sei][1]]; + mesh[sej][1] = mapto[mesh[sej][1]]; + mesh[sej][0] = mapto[mesh[sej][0]]; + } } - } + } else { - mesh.Point(i).SetPoint(pointtomove + layerht * growthvectors.Elem(i)); + // check if it doesn't contain the other edge as well + // and if it doesn't contain both mark them as done and + // if necessary map them + for(SegmentIndex sej = 0; sej 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; + BitArray domains; + 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..fdfb2c40 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,112 @@ 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(*get_if(&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_if(&thickness); + for(auto val : thicknesses) + blp.heights.Append(val.cast()); + } + + auto prismlayers = blp.heights.Size(); + auto first_new_mat = self.GetNDomains() + 1; + auto max_dom_nr = first_new_mat; + if(string* pmaterial = get_if(&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(&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.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(&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(&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. - .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 : 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) { 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..1733323d --- /dev/null +++ b/tests/pytest/test_boundarylayer.py @@ -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 +