From 177ecc74594456fa66a438099c3ed3a722a3f3b4 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Wed, 24 Jun 2020 06:41:06 +0000 Subject: [PATCH] Allow curving of mesh if boundarylayer is flat. If surfnr is larger than nr of surfaces then do linear interpolation for PointInBetween and so on. Some fixes in boundarylayer so that surface numbers are correct. --- libsrc/csg/csgeom.cpp | 1 + libsrc/meshing/boundarylayer.cpp | 132 +++++++++++++------------------ libsrc/meshing/python_mesh.cpp | 18 ++++- 3 files changed, 74 insertions(+), 77 deletions(-) diff --git a/libsrc/csg/csgeom.cpp b/libsrc/csg/csgeom.cpp index 4874f0a9..cc3e4464 100644 --- a/libsrc/csg/csgeom.cpp +++ b/libsrc/csg/csgeom.cpp @@ -130,6 +130,7 @@ namespace netgen Point<3> & newp, EdgePointGeomInfo & newgi) const { Point<3> hnewp = p1+secpoint*(p2-p1); + //(*testout) << "hnewp " << hnewp << " s1 " << surfi1 << " s2 " << surfi2 << endl; if (surfi1 != -1 && surfi2 != -1 && surfi1 != surfi2) { diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index 55960395..3ce5891a 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -159,7 +159,8 @@ namespace netgen 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, + // -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); @@ -271,55 +272,53 @@ namespace netgen if(blp.grow_edges) for(SegmentIndex sei = 0; sei < nseg; sei++) { - PointIndex seg_p1 = mesh[sei][0]; - PointIndex seg_p2 = mesh[sei][1]; + auto& segi = mesh[sei]; // 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)) { - if(blp.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++) + if(blp.surfid.Contains(segi.si)) { - PointIndex segpair_p1 = mesh[sej][1]; - PointIndex segpair_p2 = mesh[sej][0]; + // clear the bit to indicate that this segment has been processed + segsel.Clear(sei); - // 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))) + // Find matching segment pair on other surface + for(SegmentIndex sej = 0; sej < nseg; sej++) { - // clear bit to indicate that processing of this segment is done - segsel.Clear(sej); + 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); - // Only worry about those surfaces which are not in the - // boundary layer list - if(!blp.surfid.Contains(mesh[sej].si)) + // 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(segpair_p1); - auto pnt2_elems = meshtopo.GetVertexSurfaceElements(segpair_p2); + 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() == mesh[sej].si) + 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(segpair_p1) + " to " + ToString(segpair_p2)); + throw Exception("Couldn't find element on other side for " + ToString(segj[0]) + " to " + ToString(segj[1])); const auto& commsel = mesh[pnt_commelem]; auto surfelem_vect = GetSurfaceNormal(mesh, commsel); if(blp.outside) surfelem_vect *= -1; 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; @@ -332,7 +331,7 @@ namespace netgen { 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, + FaceDescriptor fd(-1, get<1>(domains), get<2>(domains), -1); @@ -367,46 +366,37 @@ namespace netgen 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]; + Segment s1 = segi; + Segment s2 = segj; 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); - } + if(blp.surfid.Contains(segj.si)) + s2.si = map_surface_index(segj.si); else { - s2.si = domains_to_surf_index[make_tuple(s2.si, - blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(s2.si).DomainOut())]; + auto side_surf = domains_to_surf_index[make_tuple(s2.si, blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(s2.si).DomainOut())]; + if(blp.outside) + s2.si = side_surf; + else + segj.si = side_surf; } s1.si = map_surface_index(s1.si); - if(create_it) - { - mesh.AddSegment(s1); - mesh.AddSegment(s2); - } + s1.surfnr1 = s1.surfnr2 = s2.surfnr1 = s2.surfnr2 = -1; + mesh.AddSegment(s1); + mesh.AddSegment(s2); + } + // segi[0] = mapto[segi[0]] not working somehow? + mesh[sei][0] = mapto[segi[0]]; + mesh[sei][1] = mapto[segi[1]]; + mesh[sej][0] = mapto[segj[0]]; + mesh[sej][1] = mapto[segj[1]]; } - - // 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 { @@ -457,9 +447,9 @@ namespace netgen { for(SurfaceElementIndex si = 0; si < nse; si++) { - if(blp.surfid.Contains(mesh[si].GetIndex())) + const auto& sel = mesh[si]; + if(blp.surfid.Contains(sel.GetIndex())) { - const auto& sel = mesh[si]; Element2d newel = sel; newel.SetIndex(map_surface_index(sel.GetIndex())); mesh.AddSurfaceElement(newel); @@ -548,36 +538,28 @@ namespace netgen for(SurfaceElementIndex sei : Range(nse)) { auto& sel = mesh[sei]; - bool to_move = blp.surfid.Contains(sel.GetIndex()); - if(blp.domains.Size()) + if(!blp.surfid.Contains(sel.GetIndex())) { - if(blp.outside) - to_move |= blp.domains[mesh.GetFaceDescriptor(sel.GetIndex()).DomainIn()]; - else - to_move |= !blp.domains[mesh.GetFaceDescriptor(sel.GetIndex()).DomainIn()]; - } - - if(to_move) - { - for(auto& pnum : sel.PNums()) - // Check (Doublecheck) if the corresponding point has a - // copy available for remapping - if(mapto[pnum].IsValid()) - // Map the surface elements to the new points - pnum = mapto[pnum]; + 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]; - bool to_move = blp.outside ? blp.domains[el.GetIndex()] : !blp.domains[el.GetIndex()]; - if(blp.domains.Size() == 0 || to_move) + // only move the elements on the correct side + if(blp.outside ? blp.domains[el.GetIndex()] : !blp.domains[el.GetIndex()]) for(auto& pnum : el.PNums()) - // Check (Doublecheck) if the corresponding point has a - // copy available for remapping if(mapto[pnum].IsValid()) - // Map the volume elements to the new points pnum = mapto[pnum]; } diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index fdfb2c40..60c719e7 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -949,8 +949,22 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) { 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); + { + auto& fd = self.GetFaceDescriptor(i); + if(regex_match(fd.GetBCName(), pattern)) + { + auto dom_pattern = get_if(&domain); + // only add if adjacent to domain + if(dom_pattern) + { + regex pattern(*dom_pattern); + if(regex_match(self.GetMaterial(fd.DomainIn()), pattern) || (fd.DomainOut() > 0 ? regex_match(self.GetMaterial(fd.DomainOut()), pattern) : false)) + blp.surfid.Append(i); + } + else + blp.surfid.Append(i); + } + } } if(double* pthickness = get_if(&thickness); pthickness)