From 609cbbcadfb23e90ff9a619ec22daa4b8dc4c0ad Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Tue, 17 Nov 2020 18:43:39 +0100 Subject: [PATCH] rewrite create boundarylayer function (now more efficient and stable and easier) --- libsrc/meshing/boundarylayer.cpp | 968 ++++++++++------------------- libsrc/meshing/boundarylayer.hpp | 2 +- libsrc/meshing/python_mesh.cpp | 39 +- tests/pytest/test_boundarylayer.py | 14 +- 4 files changed, 359 insertions(+), 664 deletions(-) diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index 81c466a8..63c90ba9 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -88,663 +88,385 @@ namespace netgen cout << "Quads: " << nq << endl; } + void GenerateBoundaryLayer(Mesh& mesh, const BoundaryLayerParameters& blp) + { + int max_edge_nr = -1; + for(const auto& seg : mesh.LineSegments()) + if(seg.edgenr > max_edge_nr) + max_edge_nr = seg.edgenr; + int new_mat_nr = mesh.GetNDomains() +1; + mesh.SetMaterial(new_mat_nr, blp.new_mat); + auto domains = blp.domains; + if(!blp.outside) + domains.Invert(); + mesh.UpdateTopology(); + auto& meshtopo = mesh.GetTopology(); -/* - Philippose Rajan - 11 June 2009 + int np = mesh.GetNP(); + int ne = mesh.GetNE(); + int nse = mesh.GetNSE(); + int nseg = mesh.GetNSeg(); - Function to calculate the surface normal at a given - vertex of a surface element, with respect to that - surface element. + Array, PointIndex> mapto(np); - This function is used by the boundary layer generation - function, in order to calculate the effective direction - in which the prismatic layer should grow -*/ - inline Vec<3> GetSurfaceNormal(Mesh & mesh, const Element2d & el) - { - auto v0 = mesh[el[0]]; - auto v1 = mesh[el[1]]; - auto v2 = mesh[el[2]]; - Vec<3> vec1 = v1-v0; - Vec<3> vec2 = v2-v0; - Vec<3> normal = Cross(vec1, vec2); - normal.Normalize(); - return normal; - } + Array, PointIndex> growthvectors(np); + growthvectors = 0.; -/* - Philippose Rajan - 11 June 2009 - modified by Christopher Lackner Apr 2020 - - Added an initial experimental function for - generating prismatic boundary layers on - a given set of surfaces. - - The number of layers, height of the first layer - and the growth / shrink factor can be specified - by the user + Array surfacefacs(mesh.GetNFD()+1); + surfacefacs = 0.; - Currently, the layer height is calculated using: - height = h_first_layer * (growth_factor^(num_layers - 1)) -*/ - void GenerateBoundaryLayer (Mesh & mesh, const BoundaryLayerParameters & blp) - { - PrintMessage(1, "Generating boundary layer..."); - PrintMessage(3, "Old NP: ", mesh.GetNP()); - PrintMessage(3, "Old NSE: ",mesh.GetNSE()); + auto getSurfaceNormal = [&mesh] (const Element2d& el) + { + auto v0 = mesh[el[0]]; + return Cross(mesh[el[1]]-v0, mesh[el[2]]-v0).Normalize(); + }; - map, int> domains_to_surf_index; - map, int> pi_to_edgenr; + // surface index map + Array si_map(mesh.GetNFD()+1); + si_map = -1; - map last_layer_surface_index_map; - int max_surface_index = mesh.GetNFD(); + int fd_old = mesh.GetNFD(); - int max_edge_nr = -1; - for(const auto& seg : mesh.LineSegments()) - if(seg.edgenr > max_edge_nr) - max_edge_nr = seg.edgenr; - - for(int layer = blp.heights.Size(); layer >= 1; layer--) - { - PrintMessage(3, "Generating layer: ", layer); - - auto map_surface_index = [&](auto si) - { - if(last_layer_surface_index_map.find(si) == last_layer_surface_index_map.end()) - { - last_layer_surface_index_map[si] = ++max_surface_index; - auto& old_fd = mesh.GetFaceDescriptor(si); - int domout = blp.outside ? old_fd.DomainOut() : blp.new_matnrs[layer-1]; - int domin = blp.outside ? blp.new_matnrs[layer-1] : old_fd.DomainIn(); - // -1 surf nr is so that curving does not do anything - FaceDescriptor fd(-1, - domin, domout, -1); - fd.SetBCProperty(max_surface_index); - mesh.AddFaceDescriptor(fd); - mesh.SetBCName(max_surface_index-1, - "mapped_" + old_fd.GetBCName()); - return max_surface_index; - } - return last_layer_surface_index_map[si]; - }; - - mesh.UpdateTopology(); - 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())) + // create new FaceDescriptors + for(auto i : Range(1, fd_old+1)) + { + auto& fd = mesh.GetFaceDescriptor(i); + if(blp.surfid.Contains(i)) + { + if(auto isIn = domains.Test(fd.DomainIn()); isIn != domains.Test(fd.DomainOut())) { - auto n2 = GetSurfaceNormal(mesh,sel); - auto& fd = mesh.GetFaceDescriptor(sel.GetIndex()); - auto domin = fd.DomainIn(); - auto domout = fd.DomainOut(); - if(blp.domains.Test(domout) && !blp.domains.Test(domin)) - n2 *= -1; - if(!blp.outside) - n2 *= -1; - for(auto pi : sel.PNums()) - { - // 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); - } + int new_si = mesh.GetNFD()+1; + surfacefacs[i] = isIn ? 1. : -1.; + // -1 surf nr is so that curving does not do anything + FaceDescriptor new_fd(-1, isIn ? new_mat_nr : fd.DomainIn(), + isIn ? fd.DomainOut() : new_mat_nr, -1); + new_fd.SetBCProperty(new_si); + mesh.AddFaceDescriptor(new_fd); + si_map[i] = new_si; + mesh.SetBCName(new_si-1, "mapped_" + fd.GetBCName()); } + } + } - // project growthvector on surface for inner angles - for(const auto& sel : mesh.SurfaceElements()) - // if(!blp.surfid.Contains(sel.GetIndex())) - // { - if(blp.project_boundaries.Contains(sel.GetIndex())) + // mark points for remapping + for(const auto& sel : mesh.SurfaceElements()) + { + auto n = surfacefacs[sel.GetIndex()] * getSurfaceNormal(sel); + if(n.Length2() != 0.) + { + for(auto pi : sel.PNums()) { - auto n = GetSurfaceNormal(mesh, sel); - for(auto i : Range(sel.PNums())) - { - auto pi = sel.PNums()[i]; - if(growthvectors[pi].Length2() == 0.) - continue; - auto next = sel.PNums()[(i+1)%sel.GetNV()]; - auto prev = sel.PNums()[i == 0 ? sel.GetNV()-1 : i-1]; - auto v1 = (mesh[next] - mesh[pi]).Normalize(); - auto v2 = (mesh[prev] - mesh[pi]).Normalize(); - auto v3 = growthvectors[pi]; - v3.Normalize(); - // only project if sel goes in boundarylayer direction - if((v1 * v3 < 1e-12) || (v2 * v3 < 1e-12)) - continue; - - auto& g = growthvectors[pi]; - auto ng = n * g; - auto gg = g * g; - auto nn = n * n; - if(fabs(ng*ng-nn*gg) < 1e-12 || fabs(ng) < 1e-12) continue; - auto a = -ng*ng/(ng*ng-nn * gg); - auto b = ng*gg/(ng*ng-nn*gg); - g += a*g + b*n; - } + auto & np = growthvectors[pi]; + if(np.Length() == 0) { np = n; continue; } + auto npn = np * n; + auto npnp = np * np; + auto nn = n * n; + if(nn-npn*npn/npnp == 0) { np = n; continue; } + np += (nn - npn)/(nn - npn*npn/npnp) * (n - npn/npnp * np); } + } + } - if (!blp.grow_edges) - { - for(const auto& sel : mesh.LineSegments()) - { - int count = 0; - for(const auto& sel2 : mesh.LineSegments()) - if(((sel[0] == sel2[0] && sel[1] == sel2[1]) || (sel[0] == sel2[1] && sel[1] == sel2[0])) && blp.surfid.Contains(sel2.si)) - count++; - if(count == 1) - { - bndnodes.Clear(sel[0]); - bndnodes.Clear(sel[1]); - } - } - } + // Bit array to keep track of segments already processed + BitArray segs_done(nseg); + segs_done.Clear(); - // 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..."); + // map for all segments with same points + // points to pair of SegmentIndex, int + // int is type of other segment, either: + // 0 == adjacent surface grows layer + // 1 == adjacent surface doesn't grow layer, but layer ends on it + // 2 == adjacent surface is interior surface that ends on layer + // 3 == adjacent surface is exterior surface that ends on layer (not allowed yet) + Array>, SegmentIndex> segmap(mesh.GetNSeg()); - for (PointIndex pi = 1; pi <= np; pi++) - { - if (bndnodes.Test(pi)) - mapto[pi] = mesh.AddPoint(mesh[pi]); - else - mapto[pi].Invalidate(); - } + // moved segments + Array moved_segs; - // Add quad surface elements at edges for surfaces which - // don't have boundary layers + // boundaries to project endings to + BitArray project_boundaries(fd_old+1); + BitArray move_boundaries(fd_old+1); + project_boundaries.Clear(); + move_boundaries.Clear(); - // Bit array to keep track of segments already processed - BitArray segsel(nseg); + Array seg2surfel(mesh.GetNSeg()); + for(auto si : Range(mesh.SurfaceElements())) + { + NgArray surfeledges; + meshtopo.GetSurfaceElementEdges(si+1, surfeledges); + for(auto edgenr : surfeledges) + for(auto sei : Range(mesh.LineSegments())) + if(meshtopo.GetEdge(sei)+1 == edgenr && + mesh[sei].si == mesh[si].GetIndex()) + seg2surfel[sei] = si; + } - // Set them all to "1" to initially activate all segments - segsel.Set(); - - // remove double segments (if more than 2 surfaces come together - // in one edge. If one of them is mapped, keep that one and - // map the others to it. - Array> segmap(nseg); - for(SegmentIndex sei = 0; sei < nseg; sei++) - { - if(!segsel.Test(sei)) continue; - const auto& segi = mesh[sei]; - for(SegmentIndex sej = 0; sej < nseg; sej++) - { - if(sej == sei || !segsel.Test(sej)) continue; - const auto& segj = mesh[sej]; - if(segi[0] == segj[0] && segi[1] == segj[1]) - { - SegmentIndex main, other; - if(blp.surfid.Contains(segi.si)) - { main = sei; other = sej; } - else { main = sej; other = sei; } - segsel.Clear(other); - for(auto& s : segmap[other]) - segmap[main].Append(s); - segmap[other].SetSize(0); - segmap[main].Append(other); - if(other == sei) sej = nseg; - } - } - } - - PrintMessage(3, "Adding 2D Quad elements on required surfaces..."); - - if(blp.grow_edges) - for(SegmentIndex sei = 0; sei < nseg; sei++) + for(auto si : Range(mesh.LineSegments())) + { + if(segs_done[si]) continue; + const auto& segi = mesh[si]; + if(si_map[segi.si] == -1) continue; + segs_done.SetBit(si); + segmap[si].Append(make_pair(si, 0)); + moved_segs.Append(si); + for(auto sj : Range(mesh.LineSegments())) + { + if(segs_done.Test(sj)) continue; + const auto& segj = mesh[sj]; + if((segi[0] == segj[0] && segi[1] == segj[1]) || + (segi[0] == segj[1] && segi[1] == segj[0])) { - // 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)) + segs_done.SetBit(sj); + int type; + if(si_map[segj.si] != -1) + type = 0; + else if(const auto& fd = mesh.GetFaceDescriptor(segj.si); domains.Test(fd.DomainIn()) && domains.Test(fd.DomainOut())) { - // copy here since we will add segments and this would - // invalidate a reference! - auto segi = mesh[sei]; - if(blp.surfid.Contains(segi.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++) - { - // copy here since we will add segments and this would - // invalidate a reference! - auto segj = mesh[sej]; - // 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) && ((segi[0] == segj[1]) && (segi[1] == segj[0]))) - { - // clear bit to indicate that processing of this segment is done - segsel.Clear(sej); - - // if segj is not in surfel list we nned to add quads - if(!blp.surfid.Contains(segj.si)) - { - SurfaceElementIndex pnt_commelem; - SetInvalid(pnt_commelem); - - auto pnt1_elems = meshtopo.GetVertexSurfaceElements(segj[0]); - auto pnt2_elems = meshtopo.GetVertexSurfaceElements(segj[1]); - - for(auto pnt1_sei : pnt1_elems) - if(mesh[pnt1_sei].GetIndex() == segj.si) - for(auto pnt2_sei : pnt2_elems) - if(pnt1_sei == pnt2_sei) - pnt_commelem = pnt1_sei; - - if(IsInvalid(pnt_commelem)) - throw Exception("Couldn't find element on other side for " + ToString(segj[0]) + " to " + ToString(segj[1])); - - const auto& commsel = mesh[pnt_commelem]; - Element2d sel(QUAD); - auto seg_p1 = segi[0]; - auto seg_p2 = segi[1]; - if(blp.outside) - Swap(seg_p1, seg_p2); - sel[0] = seg_p1; - sel[1] = seg_p2; - sel[2] = mapto[seg_p2]; - sel[3] = mapto[seg_p1]; - auto domains = make_tuple(commsel.GetIndex(), blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(commsel.GetIndex()).DomainOut()); - - if(domains_to_surf_index.find(domains) == domains_to_surf_index.end()) - { - domains_to_surf_index[domains] = ++max_surface_index; - domains_to_surf_index[make_tuple(max_surface_index, get<1>(domains), get<2>(domains))] = max_surface_index; - FaceDescriptor fd(-1, - get<1>(domains), - get<2>(domains), - -1); - fd.SetBCProperty(max_surface_index); - mesh.AddFaceDescriptor(fd); - mesh.SetBCName(max_surface_index-1, - mesh.GetBCName(get<0>(domains)-1)); - } - auto new_index = domains_to_surf_index[domains]; - 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); - } - - // in last layer insert new segments - if(layer == blp.heights.Size()) - { - max_edge_nr++; - if(!blp.surfid.Contains(segj.si)) - { - Segment s3 = segj; - s3.si = map_surface_index(segj.si)-1; - Swap(s3[0], s3[1]); - if(blp.outside) - { - s3[0] = mapto[s3[0]]; - s3[1] = mapto[s3[1]]; - } - else - s3.edgenr = max_edge_nr; - mesh.AddSegment(s3); - } - Segment s1 = segi; - Segment s2 = segj; - s1.edgenr = max_edge_nr; - s2.edgenr = max_edge_nr; - auto side_surf = domains_to_surf_index[make_tuple(s2.si, blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(s2.si).DomainOut())]; - if(blp.surfid.Contains(segj.si)) - s2.si = map_surface_index(segj.si); - else - { - if(blp.outside) - { - s2.si = side_surf; - } - else - mesh[sej].si = side_surf; - } - s1.si = map_surface_index(s1.si); - s1.surfnr1 = s1.surfnr2 = s2.surfnr1 = s2.surfnr2 = -1; - mesh.AddSegment(s1); - mesh.AddSegment(s2); - } - - segmap.SetSize(mesh.LineSegments().Size()); - for(auto sei2 : segmap[sei]) - { - auto& s = mesh[sei2]; - if(blp.outside && layer == blp.heights.Size()) - { - if(blp.surfid.Contains(s.si)) - s.si = map_surface_index(s.si); - s.edgenr = max_edge_nr; - } - else - { - s[0] = mapto[s[0]]; - s[1] = mapto[s[1]]; - } - } - for(auto sej2 : segmap[sej]) - { - auto& s = mesh[sej2]; - if(blp.outside && layer == blp.heights.Size()) - { - if(blp.surfid.Contains(s.si)) - s.si = map_surface_index(s.si); - s.edgenr = max_edge_nr; - } - else - { - s[0] = mapto[s[0]]; - s[1] = mapto[s[1]]; - } - } - - // do not use segi (not even with reference, since - // mesh.AddSegment will resize segment array and - // invalidate reference), this is why we copy it!!! - mesh[sei][0] = mapto[segi[0]]; - mesh[sei][1] = mapto[segi[1]]; - mesh[sej][0] = mapto[segj[0]]; - mesh[sej][1] = mapto[segj[1]]; - } - } + type = 2; + if(fd.DomainIn() == 0 || fd.DomainOut() == 0) + project_boundaries.SetBit(segj.si); + } + else if(const auto& fd = mesh.GetFaceDescriptor(segj.si); !domains.Test(fd.DomainIn()) && !domains.Test(fd.DomainOut())) + { + type = 3; + if(fd.DomainIn() == 0 || fd.DomainOut() == 0) + project_boundaries.SetBit(segj.si); + move_boundaries.SetBit(segj.si); } else { - // 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 1e-12) || (v2 * v3 > 1e-12)) + in_surface_direction.SetBit(sel.GetIndex()); + + auto& g = growthvectors[pi]; + auto ng = n * g; + auto gg = g * g; + auto nn = n * n; + // if(fabs(ng*ng-nn*gg) < 1e-12 || fabs(ng) < 1e-12) continue; + auto a = -ng*ng/(ng*ng-nn * gg); + auto b = ng*gg/(ng*ng-nn*gg); + g += a*g + b*n; + } + } + } + else + { + for(const auto& seg : mesh.LineSegments()) + { + int count = 0; + for(const auto& seg2 : mesh.LineSegments()) + if(((seg[0] == seg2[0] && seg[1] == seg2[1]) || (seg[0] == seg2[1] && seg[1] == seg2[0])) && blp.surfid.Contains(seg2.si)) + count++; + if(count == 1) + { + growthvectors[seg[0]] = {0., 0., 0.}; + growthvectors[seg[1]] = {0., 0., 0.}; + } + } + } + + // insert new points + for (PointIndex pi = 1; pi <= np; pi++) + if (growthvectors[pi].Length2() != 0) + { + Point<3> p = mesh[pi]; + for(auto i : Range(blp.heights)) + { + p += blp.heights[i] * growthvectors[pi]; + mapto[pi].Append(mesh.AddPoint(p)); + } + } + + // add 2d quads on required surfaces + map, int> seg2edge; + if(blp.grow_edges) + { + for(auto sei : moved_segs) + { + // copy here since we will add segments and this would + // invalidate a reference! + auto segi = mesh[sei]; + for(auto [sej, type] : segmap[sei]) + { + auto segj = mesh[sej]; + if(type == 0) + { + Segment s; + s[0] = mapto[segj[0]].Last(); + s[1] = mapto[segj[1]].Last(); + s[2] = PointIndex::INVALID; + auto pair = s[0] < s[1] ? make_pair(s[0], s[1]) : make_pair(s[1], s[0]); + if(seg2edge.find(pair) == seg2edge.end()) + seg2edge[pair] = ++max_edge_nr; + s.edgenr = seg2edge[pair]; + s.si = si_map[segj.si]; + mesh.AddSegment(s); + } + // here we need to grow the quad elements + else if(type == 1) + { + PointIndex pp1 = segj[1]; + PointIndex pp2 = segj[0]; + if(in_surface_direction.Test(segj.si)) + { + Swap(pp1, pp2); + move_boundaries.SetBit(segj.si); + } + PointIndex p1 = pp1; + PointIndex p2 = pp2; + PointIndex p3, p4; + Segment s0; + s0[0] = p1; + s0[1] = p2; + s0[2] = PointIndex::INVALID; + s0.edgenr = segj.edgenr; + s0.si = segj.si; + mesh.AddSegment(s0); + + for(auto i : Range(blp.heights)) + { + Element2d sel(QUAD); + p3 = mapto[pp2][i]; + p4 = mapto[pp1][i]; + sel[0] = p1; + sel[1] = p2; + sel[2] = p3; + sel[3] = p4; + sel.SetIndex(segj.si); + mesh.AddSurfaceElement(sel); + + // TODO: Too many, would be enough to only add outermost ones + Segment s1; + s1[0] = p2; + s1[1] = p3; + s1[2] = PointIndex::INVALID; + auto pair = make_pair(p2, p3); + if(seg2edge.find(pair) == seg2edge.end()) + seg2edge[pair] = ++max_edge_nr; + s1.edgenr = seg2edge[pair]; + s1.si = segj.si; + mesh.AddSegment(s1); + Segment s2; + s2[0] = p4; + s2[1] = p1; + s2[2] = PointIndex::INVALID; + pair = make_pair(p1, p4); + if(seg2edge.find(pair) == seg2edge.end()) + seg2edge[pair] = ++max_edge_nr; + s2.edgenr = seg2edge[pair]; + s2.si = segj.si; + mesh.AddSegment(s2); + p1 = p4; + p2 = p3; + } + Segment s3; + s3[0] = p3; + s3[1] = p4; + s3[2] = PointIndex::INVALID; + auto pair = p3 < p4 ? make_pair(p3, p4) : make_pair(p4, p3); + if(seg2edge.find(pair) == seg2edge.end()) + seg2edge[pair] = ++max_edge_nr; + s3.edgenr = seg2edge[pair]; + s3.si = segj.si; + mesh.AddSegment(s3); } } + } + } - // add surface elements between layer and old domain - if(layer == blp.heights.Size()) - { - for(SurfaceElementIndex si = 0; si < nse; si++) - { - const auto& sel = mesh[si]; - if(blp.surfid.Contains(sel.GetIndex())) - { - Element2d newel = sel; - newel.SetIndex(map_surface_index(sel.GetIndex())); - mesh.AddSurfaceElement(newel); - } - } - } + for(SurfaceElementIndex si = 0; si < nse; si++) + { + auto& sel = mesh[si]; + if(si_map[sel.GetIndex()] != -1) + { + Array points(sel.PNums()); + if(surfacefacs[sel.GetIndex()] > 0) Swap(points[0], points[2]); + for(auto j : Range(blp.heights)) + { + auto eltype = points.Size() == 3 ? PRISM : HEX; + Element el(eltype); + for(auto i : Range(points)) + el[i] = points[i]; + for(auto i : Range(points)) + points[i] = mapto[sel.PNums()[i]][j]; + if(surfacefacs[sel.GetIndex()] > 0) Swap(points[0], points[2]); + for(auto i : Range(points)) + el[sel.PNums().Size() + i] = points[i]; + el.SetIndex(new_mat_nr); + mesh.AddVolumeElement(el); + } + Element2d newel = sel; + for(auto& p : newel.PNums()) + p = mapto[p].Last(); + newel.SetIndex(si_map[sel.GetIndex()]); + mesh.AddSurfaceElement(newel); + } + if(move_boundaries.Test(sel.GetIndex())) + for(auto& p : sel.PNums()) + if(mapto[p].Size()) + p = mapto[p].Last(); + } - // Add prismatic cells at the boundaries - PrintMessage(3, "Generating prism boundary layer volume elements..."); + for(SegmentIndex sei = 0; sei < nseg; sei++) + { + auto& seg = mesh[sei]; + if(move_boundaries.Test(seg.si)) + for(auto& p : seg.PNums()) + if(mapto[p].Size()) + p = mapto[p].Last(); + } - for (SurfaceElementIndex si = 0; si < nse; si++) - { - const auto& sel = mesh[si]; - if(blp.surfid.Contains(sel.GetIndex())) - { - int classify = 0; - for(auto j : Range(sel.PNums())) - if (mapto[sel[j]].IsValid()) - classify += (1 << j); + for(ElementIndex ei = 0; ei < ne; ei++) + { + auto& el = mesh[ei]; + if(!domains[el.GetIndex()]) + { + for(auto& p : el.PNums()) + if(mapto[p].Size()) + p = mapto[p].Last(); + } + } - if(classify == 0) - continue; - - Element el; - - if(sel.GetType() == TRIG) - { - ELEMENT_TYPE types[] = { PRISM, TET, TET, PYRAMID, - TET, PYRAMID, PYRAMID, PRISM }; - int nums[] = { sel[0], sel[1], sel[2], mapto[sel[0]], mapto[sel[1]], mapto[sel[2]] }; - int vertices[][6] = - { - { 0, 1, 2, 0, 1, 2 }, // should not occur - { 0, 2, 1, 3, 0, 0 }, - { 0, 2, 1, 4, 0, 0 }, - { 0, 1, 4, 3, 2, 0 }, - - { 0, 2, 1, 5, 0, 0 }, - { 2, 0, 3, 5, 1, 0 }, - { 1, 2, 5, 4, 0, 0 }, - { 0, 2, 1, 3, 5, 4 } - }; - if(blp.outside) - { - if(classify != 7) - throw Exception("Outside with non prisms not yet implemented"); - for(auto i : Range(6)) - vertices[7][i] = i; - } - - el = Element(types[classify]); - for(auto i : Range(el.PNums())) - el.PNums()[i] = nums[vertices[classify][i]]; - } - else // sel.GetType() == QUAD - { - int nums[] = { sel[0], sel[1], sel[2], sel[3], - mapto[sel[0]], mapto[sel[1]], - mapto[sel[2]], mapto[sel[3]] }; - ArrayMem vertices; - switch(classify) - { - case 6: - { - if(blp.outside) - throw Exception("Type 6 quad outside layer is not yet implemented!"); - el = Element(PRISM); - vertices = {0, 1, 5, 3, 2, 6}; - break; - } - case 9: - { - if(blp.outside) - throw Exception("Type 9 quad outside layer is not yet implemented!"); - el = Element(PRISM); - vertices = { 1, 4, 0, 2, 7, 3 }; - break; - } - case 15: - { - vertices = { 0, 1, 2, 3, 4, 5, 6, 7 }; - if(!blp.outside) - { - Swap(vertices[1], vertices[3]); - Swap(vertices[5], vertices[7]); - } - el = Element(HEX); - break; - } - default: - throw Exception("Type " + ToString(classify) + " for quad layer not yet implemented!"); - } - for(auto i : Range(el.PNums())) - el.PNums()[i] = nums[vertices[i]]; - } - el.SetIndex(blp.new_matnrs[layer-1]); - mesh.AddVolumeElement(el); - } - } - - // 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.surfid.Contains(sel.GetIndex())) - { - const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex()); - if(blp.outside && - (!blp.domains[fd.DomainIn()] && !blp.domains[fd.DomainOut()])) - continue; - if(!blp.outside && - (blp.domains[fd.DomainIn()] || blp.domains[fd.DomainOut()])) - continue; - } - for(auto& pnum : sel.PNums()) - if(mapto[pnum].IsValid()) - pnum = mapto[pnum]; - } - - for(ElementIndex ei : Range(ne)) - { - auto& el = mesh[ei]; - // only move the elements on the correct side - if(blp.outside ? blp.domains[el.GetIndex()] : !blp.domains[el.GetIndex()]) - for(auto& pnum : el.PNums()) - if(mapto[pnum].IsValid()) - pnum = mapto[pnum]; - } - - // Lock all the prism points so that the rest of the mesh can be - // optimised without invalidating the entire mesh - // for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++) - for (PointIndex pi = 1; pi <= np; pi++) - if(bndnodes.Test(pi)) mesh.AddLockedPoint(pi); - - // Now, actually pull back the old surface points to create - // the actual boundary layers - PrintMessage(3, "Moving and optimising boundary layer points..."); - - for (PointIndex i = 1; i <= np; i++) - { - if(bndnodes.Test(i)) - { - MeshPoint pointtomove; - pointtomove = mesh.Point(i); - mesh.Point(i).SetPoint(pointtomove + layerht * growthvectors[i]); - } - } - mesh.Compress(); - } - - for(int i=1; i <= mesh.GetNFD(); i++) + for(auto i : Range(1, fd_old+1)) + if(si_map[i] != -1) { - auto& fd = mesh.GetFaceDescriptor(i); - if(blp.surfid.Contains(fd.BCProperty())) - { - if(blp.outside) - fd.SetDomainOut(blp.new_matnrs[blp.new_matnrs.Size()-1]); - else - fd.SetDomainIn(blp.new_matnrs[blp.new_matnrs.Size()-1]); - } + if(mesh.GetFaceDescriptor(mesh.GetNFD()).DomainIn() == new_mat_nr) + mesh.GetFaceDescriptor(i).SetDomainOut(new_mat_nr); + else + mesh.GetFaceDescriptor(i).SetDomainIn(new_mat_nr); } - - PrintMessage(3, "New NP: ", mesh.GetNP()); - PrintMessage(1, "Boundary Layer Generation....Done!"); - } + } } diff --git a/libsrc/meshing/boundarylayer.hpp b/libsrc/meshing/boundarylayer.hpp index 113d9b1f..1d33fab1 100644 --- a/libsrc/meshing/boundarylayer.hpp +++ b/libsrc/meshing/boundarylayer.hpp @@ -14,7 +14,7 @@ public: // parameters by Philippose .. Array surfid; Array heights; - Array new_matnrs; + string new_mat; BitArray domains; bool outside = false; // set the boundary layer on the outside bool grow_edges = false; diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index caee2d35..083ef684 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -1017,9 +1017,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def ("BoundaryLayer", [](Mesh & self, variant boundary, variant thickness, - variant material, + string material, variant domain, bool outside, - bool grow_edges, optional project_boundaries) { BoundaryLayerParameters blp; @@ -1050,6 +1049,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) } } } + blp.new_mat = material; if(project_boundaries.has_value()) { @@ -1070,34 +1070,13 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) 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 + int nr_domains = self.GetNDomains(); + blp.domains.SetSize(nr_domains + 1); // one based blp.domains.Clear(); if(string* pdomain = get_if(&domain); pdomain) { regex pattern(*pdomain); - for(auto i : Range(1, first_new_mat)) + for(auto i : Range(1, nr_domains+1)) if(regex_match(self.GetMaterial(i), pattern)) blp.domains.SetBit(i); } @@ -1106,19 +1085,15 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) 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; + blp.grow_edges = true; 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, py::arg("project_boundaries")=nullopt, + py::arg("project_boundaries")=nullopt, R"delimiter( Add boundary layer to mesh. diff --git a/tests/pytest/test_boundarylayer.py b/tests/pytest/test_boundarylayer.py index 27124613..04c0ee8c 100644 --- a/tests/pytest/test_boundarylayer.py +++ b/tests/pytest/test_boundarylayer.py @@ -18,7 +18,7 @@ def test_boundarylayer(outside, capfd): mesh = unit_cube.GenerateMesh(maxh=0.3) ne_before = mesh.ne layer_surfacenames = ["right", "top", "left", "back", "bottom"] - mesh.BoundaryLayer("|".join(layer_surfacenames), [0.01, 0.02], "layer", outside=outside, grow_edges=True) + mesh.BoundaryLayer("|".join(layer_surfacenames), [0.01, 0.01], "layer", outside=outside) should_ne = ne_before + 2 * GetNSurfaceElements(mesh, layer_surfacenames) assert mesh.ne == should_ne @@ -26,7 +26,7 @@ def test_boundarylayer(outside, capfd): assert not "elements are not matching" in capture.out for side in ["front"]: - mesh.BoundaryLayer(side, [0.001, 0.002], "layer", outside=outside, grow_edges=True) + mesh.BoundaryLayer(side, [0.001, 0.001], "layer", outside=outside) should_ne += 2 * GetNSurfaceElements(mesh, [side]) assert mesh.ne == should_ne capture = capfd.readouterr() @@ -53,8 +53,8 @@ def test_boundarylayer2(outside, version, capfd): geo.CloseSurfaces(top, bot, []) mesh = geo.GenerateMesh() should_ne = mesh.ne + 2 * GetNSurfaceElements(mesh, ["default"], "part") - layersize = 0.05 - mesh.BoundaryLayer("default", [0.5 * layersize, layersize], "part", domains="part", outside=outside, grow_edges=True) + layersize = 0.025 + mesh.BoundaryLayer("default", [layersize, layersize], "part", domains="part", outside=outside) assert mesh.ne == should_ne assert not "elements are not matching" in capfd.readouterr().out import netgen.gui @@ -73,8 +73,7 @@ def test_wrong_orientation(outside): mesh = geo.GenerateMesh() - mesh.BoundaryLayer(".*", 0.1, "air", domains="air", outside=outside, - grow_edges=True) + mesh.BoundaryLayer(".*", 0.1, "air", domains="air", outside=outside) ngs = pytest.importorskip("ngsolve") mesh = ngs.Mesh(mesh) assert ngs.Integrate(1, mesh) == pytest.approx(1.2**3 if outside else 1) @@ -88,8 +87,7 @@ def test_splitted_surface(): geo.Add((brick*slots).mat("slot")) mesh = geo.GenerateMesh() - mesh.BoundaryLayer(".*", [0.001, 0.002], "block", "block", outside=False, - grow_edges=True) + mesh.BoundaryLayer(".*", [0.001, 0.001], "block", "block", outside=False) ngs = pytest.importorskip("ngsolve") mesh = ngs.Mesh(mesh) assert ngs.Integrate(1, mesh) == pytest.approx(1)