diff --git a/.clang-format b/.clang-format new file mode 100644 index 00000000..24978543 --- /dev/null +++ b/.clang-format @@ -0,0 +1,64 @@ +Language: Cpp +BasedOnStyle: LLVM +AlignAfterOpenBracket: Align +AlignConsecutiveAssignments: false +AlignConsecutiveDeclarations: false +AlignEscapedNewlines: Left +AlignOperands: true +AlignTrailingComments: true +AllowAllParametersOfDeclarationOnNextLine: false +AllowShortBlocksOnASingleLine: false +AllowShortCaseLabelsOnASingleLine: true +AllowShortFunctionsOnASingleLine: InlineOnly +AllowShortIfStatementsOnASingleLine: false +AllowShortLoopsOnASingleLine: false +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakTemplateDeclarations: Yes +BinPackArguments: false +BinPackParameters: false +BreakBeforeBinaryOperators: NonAssignment +BreakBeforeBraces: Custom +BraceWrapping: + AfterClass: true + AfterControlStatement: true + AfterEnum: true + AfterFunction: true + AfterNamespace: true + AfterStruct: true + AfterUnion: true + BeforeCatch: true + BeforeElse: true + IndentBraces: true +BreakBeforeTernaryOperators: true +BreakConstructorInitializersBeforeComma: false +BreakInheritanceList: AfterColon +ColumnLimit: 0 +ConstructorInitializerAllOnOneLineOrOnePerLine: true +ConstructorInitializerIndentWidth: 2 +ContinuationIndentWidth: 2 +Cpp11BracedListStyle: true +EmptyLineBeforeAccessModifier: Never +IndentCaseLabels: false +IndentPPDirectives: None +IndentWidth: 2 +KeepEmptyLinesAtTheStartOfBlocks: false +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: None +PointerAlignment: Left +ReflowComments: true +SortIncludes: false +SpaceAfterCStyleCast: false +SpaceAfterLogicalNot: false +SpaceBeforeParens: Custom +SpaceBeforeParensOptions: + AfterControlStatements: true + AfterFunctionDefinitionName: true + AfterFunctionDeclarationName: true +SpacesInAngles: false +SpacesInContainerLiterals: false +SpacesInCStyleCastParentheses: false +SpacesInSquareBrackets: false +Standard: Latest +TabWidth: 2 +UseTab: Never + diff --git a/libsrc/csg/genmesh.cpp b/libsrc/csg/genmesh.cpp index da790fc4..c78efd8b 100644 --- a/libsrc/csg/genmesh.cpp +++ b/libsrc/csg/genmesh.cpp @@ -826,6 +826,9 @@ namespace netgen { multithread.task = "Volume meshing"; + for (int i = 0; i < geom.GetNTopLevelObjects(); i++) + mesh->SetMaterial (i+1, geom.GetTopLevelObject(i)->GetMaterial().c_str()); + MESHING3_RESULT res = MeshVolume (mparam, *mesh); @@ -838,10 +841,6 @@ namespace netgen MeshQuality3d (*mesh); - for (int i = 0; i < geom.GetNTopLevelObjects(); i++) - mesh->SetMaterial (i+1, geom.GetTopLevelObject(i)->GetMaterial().c_str()); - - #ifdef STAT_STREAM (*statout) << GetTime() << " & "; #endif diff --git a/libsrc/meshing/CMakeLists.txt b/libsrc/meshing/CMakeLists.txt index 697e8ad1..ade62d59 100644 --- a/libsrc/meshing/CMakeLists.txt +++ b/libsrc/meshing/CMakeLists.txt @@ -12,7 +12,7 @@ target_sources(nglib PRIVATE parallelmesh.cpp paralleltop.cpp basegeom.cpp python_mesh.cpp surfacegeom.cpp debugging.cpp fieldlines.cpp visual_interface.cpp - boundarylayer2d.cpp + boundarylayer2d.cpp boundarylayer_interpolate.cpp ) target_link_libraries( nglib PRIVATE $ $ ) diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index 10544c39..fd22ece8 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -271,7 +271,7 @@ namespace netgen virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, EdgePointGeomInfo* gi = nullptr) const { - if(gi && gi->edgenr < edges.Size()) + if(gi && gi->edgenr < edges.Size() && gi->edgenr >= 0) edges[gi->edgenr]->ProjectPoint(p, gi); } diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index 3eaaa508..671601f3 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -1,1659 +1,1436 @@ - -#include -#include -#include - -#include "global.hpp" -#include "debugging.hpp" - #include "boundarylayer.hpp" -#include "meshfunc.hpp" +#include "boundarylayer_limiter.hpp" +#include +#include + +#include "debugging.hpp" +#include "global.hpp" +#include "meshfunc.hpp" namespace netgen { - // checks if a segment is intersecting a plane, spanned by three points, lam will be set s.t. p_intersect = seg[0] + lam * (seg[1]-seg[0]) - bool isIntersectingPlane ( const array, 2> & seg, const array, 3> & trig, double & lam) - { - auto n = Cross(trig[1]-trig[0], trig[2]-trig[0]); - auto v0n = (seg[0]-trig[0])*n; - auto v1n = (seg[1]-trig[0])*n; - if(v0n * v1n >= 0) - return false; - lam = -v0n/(v1n-v0n); - lam *= 0.9; - if(lam < -1e-8 || lam>1+1e-8) - return false; - return true; - } +struct SpecialPointException : public Exception +{ + SpecialPointException() + : Exception("") {} +}; - bool isIntersectingPlane ( const array, 2> & seg, const ArrayMem, 4> & face, double & lam) - { - lam = 1.0; - bool intersect0 = isIntersectingPlane( seg, array, 3>{face[0], face[1], face[2]}, lam ); - if(face.Size()==3) - return intersect0; +std::tuple FindCloseVectors (FlatArray> ns, + bool find_max = true) +{ + int maxpos1 = 0; + int maxpos2 = 0; - double lam1 = 1.0; - bool intersect1 = isIntersectingPlane( seg, array, 3>{face[2], face[3], face[0]}, lam1 ); - lam = min(lam, lam1); - return intersect0 || intersect1; - } - - bool isIntersectingTrig ( const array, 2> & seg, const array, 3> & trig, double & lam) - { - if(!isIntersectingPlane(seg, trig, lam)) - return false; - - - //buffer enlargement of triangle - auto pt0 = trig[0]; - auto pt1 = trig[1]; - auto pt2 = trig[2]; - Point<3> center = { (pt0[0] + pt1[0] + pt2[0]) / 3.0, (pt0[1] + pt1[1] + pt2[1]) / 3.0, (pt0[2] + pt1[2] + pt2[2]) / 3.0 }; - array, 3> larger_trig = { - center + (pt0 - center) * 1.1, - center + (pt1 - center) * 1.1, - center + (pt2 - center) * 1.1, }; - - auto p = seg[0] + lam/0.9*(seg[1]-seg[0]); - - auto n_trig = Cross(trig[1]-trig[0], trig[2]-trig[0]).Normalize(); - for(auto i : Range(3)) + double val = find_max ? -1e99 : 1e99; + for (auto i : Range(ns)) + for (auto j : Range(i + 1, ns.Size())) { - // check if p0 and p are on same side of segment p1-p2 - auto p0 = larger_trig[i]; - auto p1 = larger_trig[(i+1)%3]; - auto p2 = larger_trig[(i+2)%3]; - // auto n = Cross(p2-p1, n_trig); - - auto v0 = (p2-p1).Normalize(); - auto v1 = (p0-p1).Normalize(); - auto inside_dir = (v1 - (v1*v0) * v0).Normalize(); - auto v2 = (p-p1).Normalize(); - if(inside_dir * v1 < 0) - inside_dir = -inside_dir; - - if( (inside_dir*v2) < 0 ) - return false; - } - return true; - }; - - bool isIntersectingFace( const array, 2> & seg, const ArrayMem, 4> & face, double & lam ) - { - lam = 1.0; - double lam0 = 1.0; - bool intersect0 = isIntersectingTrig( seg, {face[0], face[1], face[2]}, lam0 ); - if(intersect0) - lam = min(lam, lam0); - if(face.Size()==3) - return intersect0; - - double lam1 = 1.0; - bool intersect1 = isIntersectingTrig( seg, {face[2], face[3], face[0]}, lam1 ); - if(intersect1) - lam = min(lam, lam1); - return intersect0 || intersect1; - } - - array, 2> BoundaryLayerTool :: GetMappedSeg( PointIndex pi ) - { - return { mesh[pi], mesh[pi] + height*limits[pi]*growthvectors[pi] * 1.5 }; - } - - ArrayMem, 4> BoundaryLayerTool :: GetFace( SurfaceElementIndex sei ) - { - const auto & sel = mesh[sei]; - ArrayMem, 4> points(sel.GetNP()); - for(auto i : Range(sel.GetNP())) - points[i] = mesh[sel[i]]; - return points; - } - - ArrayMem, 4> BoundaryLayerTool :: GetMappedFace( SurfaceElementIndex sei ) - { - const auto & sel = mesh[sei]; - ArrayMem, 4> points(sel.GetNP()); - for(auto i : Range(sel.GetNP())) - points[i] = mesh[sel[i]] + height * limits[sel[i]]*growthvectors[sel[i]]; - return points; - } - - ArrayMem, 4> BoundaryLayerTool :: GetMappedFace( SurfaceElementIndex sei, int face ) - { - if(face == -1) return GetFace(sei); - if(face == -2) return GetMappedFace(sei); - const auto & sel = mesh[sei]; - auto np = sel.GetNP(); - auto pi0 = sel[face % np]; - auto pi1 = sel[(face+1) % np]; - ArrayMem, 4> points(4); - points[0] = points[3] = mesh[pi0]; - points[1] = points[2] = mesh[pi1]; - points[3] += height * limits[pi0]*growthvectors[pi0]; - points[2] += height * limits[pi1]*growthvectors[pi1]; - return points; - } - - Vec<3> BoundaryLayerTool :: getEdgeTangent(PointIndex pi, int edgenr) - { - Vec<3> tangent = 0.0; - ArrayMem pts; - for(auto segi : topo.GetVertexSegments(pi)) - { - auto & seg = mesh[segi]; - if(seg.edgenr != edgenr+1) - continue; - PointIndex other = seg[0]-pi+seg[1]; - if(!pts.Contains(other)) - pts.Append(other); - } - if(pts.Size() != 2) - throw Exception("Something went wrong in getEdgeTangent!"); - tangent = mesh[pts[1]] - mesh[pts[0]]; - return tangent.Normalize(); - } - - void BoundaryLayerTool :: LimitGrowthVectorLengths() - { - static Timer tall("BoundaryLayerTool::LimitGrowthVectorLengths"); RegionTimer rtall(tall); - - limits.SetSize(np); - limits = 1.0; - - // Function to calculate the dot product of two 3D vectors - // Is there netgen native function for this? - const auto Dot = [](Vec<3> a, Vec<3> b) { - return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; - }; - - auto parallel_limiter = [&](PointIndex pi1, PointIndex pi2, SurfaceElementIndex si) { - MeshPoint& a_base = mesh[pi1]; - MeshPoint& b_base = mesh[pi2]; - MeshPoint a_end = mesh[pi1] + height * limits[pi1] * growthvectors[pi1]; - MeshPoint b_end = mesh[pi2] + height * limits[pi2] * growthvectors[pi2]; - - double ab_base = (b_base - a_base).Length(); - Vec<3> a_vec = (a_end - a_base); - Vec<3> b_vec = (b_end - b_base); - - // Calculate parallel projections - Vec<3> ab_base_norm = (b_base - a_base).Normalize(); - double a_vec_x = Dot(a_vec, ab_base_norm); - double b_vec_x = Dot(b_vec, -ab_base_norm); - double ratio_parallel = (a_vec_x + b_vec_x) / ab_base; - - double PARALLEL_RATIO_LIMIT = 0.85; - if (ratio_parallel > PARALLEL_RATIO_LIMIT) { - // Adjust limits, vectors, and projections if parallel ratio exceeds the limit - double corrector = PARALLEL_RATIO_LIMIT / ratio_parallel; - limits[pi1] *= corrector; - limits[pi2] *= corrector; - } - }; - - auto perpendicular_limiter = [&](PointIndex pi1, PointIndex pi2, SurfaceElementIndex si) { - // this part is same as in parallel limiter, but note that limits contents are already changed - MeshPoint& a_base = mesh[pi1]; - MeshPoint& b_base = mesh[pi2]; - MeshPoint a_end = mesh[pi1] + height * limits[pi1] * growthvectors[pi1]; - MeshPoint b_end = mesh[pi2] + height * limits[pi2] * growthvectors[pi2]; - - double ab_base = (b_base - a_base).Length(); - Vec<3> a_vec = (a_end - a_base); - Vec<3> b_vec = (b_end - b_base); - - // Calculate parallel projections - Vec<3> ab_base_norm = (b_base - a_base).Normalize(); - double a_vec_x = Dot(a_vec, ab_base_norm); - double b_vec_x = Dot(b_vec, -ab_base_norm); - // double ratio_parallel = (a_vec_x + b_vec_x) / ab_base; - - // Calculate surface normal at point si - Vec<3> surface_normal = getNormal(mesh[si]); - - double a_vec_y = abs(Dot(a_vec, surface_normal)); - double b_vec_y = abs(Dot(b_vec, surface_normal)); - double diff_perpendicular = abs(a_vec_y - b_vec_y); - double tan_alpha = diff_perpendicular / (ab_base - a_vec_x - b_vec_x); - - double TAN_ALPHA_LIMIT = 0.36397; // Approximately 20 degrees in radians - if (tan_alpha > TAN_ALPHA_LIMIT) { - if (a_vec_y > b_vec_y) { - double correction = (TAN_ALPHA_LIMIT / tan_alpha * diff_perpendicular + b_vec_y) / a_vec_y; - limits[pi1] *= correction; - } - else { - double correction = (TAN_ALPHA_LIMIT / tan_alpha * diff_perpendicular + a_vec_y) / b_vec_y; - limits[pi2] *= correction; - } - } - }; - - auto neighbour_limiter = [&](PointIndex pi1, PointIndex pi2, SurfaceElementIndex si) { - parallel_limiter(pi1, pi2, si); - perpendicular_limiter(pi1, pi2, si); - }; - - auto modifiedsmooth = [&](size_t nsteps) { - for ([[maybe_unused]] auto i : Range(nsteps)) - for (SurfaceElementIndex sei : mesh.SurfaceElements().Range()) - { - // assuming triangle - neighbour_limiter(mesh[sei].PNum(1), mesh[sei].PNum(2), sei); - neighbour_limiter(mesh[sei].PNum(2), mesh[sei].PNum(3), sei); - neighbour_limiter(mesh[sei].PNum(3), mesh[sei].PNum(1), sei); - } - }; - - /* - auto smooth = [&] (size_t nsteps) { - for([[maybe_unused]] auto i : Range(nsteps)) - for(const auto & sel : mesh.SurfaceElements()) - { - double min_limit = 999; - for(auto pi : sel.PNums()) - min_limit = min(min_limit, limits[pi]); - for(auto pi : sel.PNums()) - limits[pi] = min(limits[pi], 1.4*min_limit); - } - }; - */ - - // check for self-intersection within new elements (prisms/hexes) - auto self_intersection = [&] () { - for(SurfaceElementIndex sei : mesh.SurfaceElements().Range()) + double ip = ns[i] * ns[j]; + if ((find_max && (ip > val)) || (!find_max && (ip < val))) { - auto facei = mesh[sei].GetIndex(); - if(facei < nfd_old && !par_surfid.Contains(facei)) - continue; - - auto sel = mesh[sei]; - auto np = sel.GetNP(); - // check if a new edge intesects the plane of any opposing face - double lam; - for(auto i : Range(np)) - for(auto fi : Range(np-2)) - if(isIntersectingPlane(GetMappedSeg(sel[i]), GetMappedFace(sei, i+fi+1), lam)) - if(lam < 1.0) - limits[sel[i]] *= lam; + val = ip; + maxpos1 = i; + maxpos2 = j; } - }; + } + return {maxpos1, maxpos2}; +} - // first step: intersect with other surface elements that are boundary of domain the layer is grown into - // second (and subsequent) steps: intersect with other boundary layers, allow restriction by 20% in each step - auto changed_domains = domains; - if(!params.outside) - changed_domains.Invert(); - - bool limit_reached = true; - double lam_lower_limit = 1.0; - int step = 0; - - while(limit_reached || step<3) +Vec<3> CalcGrowthVector (FlatArray> ns) +{ + if (ns.Size() == 0) + return {0, 0, 0}; + if (ns.Size() == 1) + return ns[0]; + if (ns.Size() == 2) { - Array new_limits; - new_limits.SetSize(np); - new_limits = 1.0; - - if(step>1) - lam_lower_limit *= 0.8; - limit_reached = false; - - // build search tree with all surface elements (bounding box of a surface element also covers the generated boundary layer) - Box<3> bbox(Box<3>::EMPTY_BOX); - for(auto pi : mesh.Points().Range()) - { - bbox.Add(mesh[pi]); - bbox.Add(mesh[pi]+limits[pi]*height*growthvectors[pi]); - } - BoxTree<3> tree(bbox); - - for(auto sei : mesh.SurfaceElements().Range()) - { - const auto & sel = mesh[sei]; - Box<3> box(Box<3>::EMPTY_BOX); - const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex()); - if(!changed_domains.Test(fd.DomainIn()) && - !changed_domains.Test(fd.DomainOut())) - continue; - for(auto pi : sel.PNums()) - box.Add(mesh[pi]); - // also add moved points to bounding box - if(par_surfid.Contains(sel.GetIndex())) - for(auto pi : sel.PNums()) - box.Add(mesh[pi]+limits[pi]*height*growthvectors[pi]); - tree.Insert(box, sei); - } - - for(auto pi : mesh.Points().Range()) - { - if(mesh[pi].Type() == INNERPOINT) - continue; - if(growthvectors[pi].Length2() == 0.0) - continue; - Box<3> box(Box<3>::EMPTY_BOX); - auto seg = GetMappedSeg(pi); - box.Add(seg[0]); - box.Add(seg[1]); - double lam = 1.0; - tree.GetFirstIntersecting(box.PMin(), box.PMax(), [&](SurfaceElementIndex sei) - { - const auto & sel = mesh[sei]; - if(sel.PNums().Contains(pi)) - return false; - auto face = GetMappedFace(sei, -2); - double lam_ = 999; - bool is_bl_sel = par_surfid.Contains(sel.GetIndex()); - - if (step == 0) - { - face = GetMappedFace(sei, -1); - if (isIntersectingFace(seg, face, lam_)) - { - if (is_bl_sel) - lam_ *= params.limit_safety; - lam = min(lam, lam_); - } - } - - if(step==1) - { - if(isIntersectingFace(seg, face, lam_)) - { - if(is_bl_sel) // allow only half the distance if the opposing surface element has a boundary layer too - lam_ *= params.limit_safety; - lam = min(lam, lam_); - } - } - // if the opposing surface element has a boundary layer, we need to additionally intersect with the new faces - if(step>1 && is_bl_sel) - { - for(auto facei : Range(-1, sel.GetNP())) - { - auto face = GetMappedFace(sei, facei); - if(isIntersectingFace(seg, face, lam_)) // && lam_ > other_limit) - { - lam = min(lam, lam_); - } - } - } - return false; - }); - if(lam<1) - { - if(lam1) - { - limit_reached = true; - lam = lam_lower_limit; - } - } - - new_limits[pi] = min(limits[pi], lam* limits[pi]); - } - step++; - limits = new_limits; - if (step > 0) - modifiedsmooth(1); + auto gw = ns[0]; + auto n = ns[1]; + auto npn = gw * n; + auto npnp = gw * gw; + auto nn = n * n; + if (fabs(nn - npn * npn / npnp) < 1e-6) + return n; + gw += (nn - npn) / (nn - npn * npn / npnp) * (n - npn / npnp * gw); + return gw; } - - self_intersection(); - modifiedsmooth(1); - - for(auto pi : Range(growthvectors)) - growthvectors[pi] *= limits[pi]; - - } - - - // depending on the geometry type, the mesh contains segments multiple times (once for each face) - bool HaveSingleSegments( const Mesh & mesh ) - { - auto& topo = mesh.GetTopology(); - NgArray surf_els; - - for(auto segi : Range(mesh.LineSegments())) - { - mesh.GetTopology().GetSegmentSurfaceElements(segi+1, surf_els); - if(surf_els.Size()<2) - continue; - - auto seg = mesh[segi]; - auto pi0 = min(seg[0], seg[1]); - auto pi1 = max(seg[0], seg[1]); - auto p0_segs = topo.GetVertexSegments(seg[0]); - - for(auto segi_other : p0_segs) - { - if(segi_other == segi) - continue; - - auto seg_other = mesh[segi_other]; - auto pi0_other = min(seg_other[0], seg_other[1]); - auto pi1_other = max(seg_other[0], seg_other[1]); - if( pi0_other == pi0 && pi1_other == pi1 ) - return false; - } - - // found segment with multiple adjacent surface elements but no other segments with same points -> have single segments - return true; - } - - return true; - } - - // duplicates segments (and sets seg.si accordingly) to have a unified data structure for all geometry types - Array BuildSegments( Mesh & mesh ) - { - Array segments; - // auto& topo = mesh.GetTopology(); - - NgArray surf_els; - - for(auto segi : Range(mesh.LineSegments())) - { - auto seg = mesh[segi]; - mesh.GetTopology().GetSegmentSurfaceElements(segi+1, surf_els); - for(auto seli : surf_els) - { - const auto & sel = mesh[seli]; - seg.si = sel.GetIndex(); - - auto np = sel.GetNP(); - for(auto i : Range(np)) - { - if(sel[i] == seg[0]) - { - if(sel[(i+1)%np] != seg[1]) - swap(seg[0], seg[1]); - break; - } - } - - segments.Append(seg); - } - } - return segments; - } - - void MergeAndAddSegments( Mesh & mesh, FlatArray new_segments) - { - INDEX_2_HASHTABLE already_added( mesh.LineSegments().Size() + 2*new_segments.Size() ); - - for(auto & seg : mesh.LineSegments()) - { - INDEX_2 i2 (seg[0], seg[1]); - i2.Sort(); - if(!already_added.Used(i2)) - already_added.Set(i2, true); - } - - for(auto & seg : new_segments) - { - INDEX_2 i2 (seg[0], seg[1]); - i2.Sort(); - - if(!already_added.Used(i2)) - { - mesh.AddSegment(seg); - already_added.Set(i2, true); - } - } - } - - void BoundaryLayerTool :: InterpolateSurfaceGrowthVectors() - { - static Timer tall("InterpolateSurfaceGrowthVectors"); RegionTimer rtall(tall); - static Timer tsmooth("InterpolateSurfaceGrowthVectors-Smoothing"); - auto np = mesh.GetNP(); - BitArray is_point_on_bl_surface(np+1); - is_point_on_bl_surface.Clear(); - BitArray is_point_on_other_surface(np+1); - is_point_on_other_surface.Clear(); - Array, PointIndex> normals(np); - for(auto pi : Range(growthvectors)) - normals[pi] = growthvectors[pi]; - - ParallelForRange( mesh.SurfaceElements().Range(), [&] ( auto myrange ) - { - for(SurfaceElementIndex sei : myrange) - { - auto facei = mesh[sei].GetIndex(); - if(facei < nfd_old && !par_surfid.Contains(facei)) - { - for(auto pi : mesh[sei].PNums()) - if(mesh[pi].Type() == SURFACEPOINT) - is_point_on_other_surface.SetBitAtomic(pi); - } - else - { - for(auto pi : mesh[sei].PNums()) - if(mesh[pi].Type() == SURFACEPOINT) - is_point_on_bl_surface.SetBitAtomic(pi); - } - } - }); - - Array points; - for(PointIndex pi : mesh.Points().Range()) - { - if(is_point_on_bl_surface[pi]) - { - points.Append(pi); - growthvectors[pi] = 0.0; - } - if(is_point_on_other_surface[pi]) - { - points.Append(pi); - } - } - - // smooth tangential part of growth vectors from edges to surface elements - RegionTimer rtsmooth(tsmooth); - for([[maybe_unused]] auto i : Range(10)) + if (ns.Size() == 3) { - for(auto pi : points) - { - auto sels = p2sel[pi]; - Vec<3> new_gw = growthvectors[pi]; - // int cnt = 1; - std::set suround; - suround.insert(pi); - auto normal = normals[pi]; - for(auto sei: sels) - { - const auto & sel = mesh[sei]; - for(auto pi1 : sel.PNums()) - if(suround.count(pi1)==0) - { - suround.insert(pi1); - auto gw_other = growthvectors[pi1]; - auto normal_other = getNormal(mesh[sei]); - auto tangent_part = gw_other - (gw_other*normal_other)*normal_other; - if(is_point_on_bl_surface[pi]) - new_gw += tangent_part; - else - new_gw += gw_other; - } - } + DenseMatrix mat(3, 3); + for (auto i : Range(3)) + for (auto j : Range(3)) + mat(i, j) = ns[i][j]; - growthvectors[pi] = 1.0/suround.size() * new_gw; + if (fabs(mat.Det()) > 1e-2) + { + DenseMatrix mat(3, 3); + for (auto i : Range(3)) + for (auto j : Range(3)) + mat(i, j) = ns[i] * ns[j]; + if (fabs(mat.Det()) > 1e-2) + { + Vector rhs(3); + rhs = 1.; + Vector res(3); + DenseMatrix inv(3, ns.Size()); + CalcInverse(mat, inv); + inv.Mult(rhs, res); + Vec<3> growth = 0.; + for (auto i : Range(ns)) + growth += res[i] * ns[i]; + return growth; + } } } + auto [maxpos1, maxpos2] = FindCloseVectors(ns); + Array> new_normals; + new_normals = ns; + const auto dot = ns[maxpos1] * ns[maxpos2]; + auto average = 0.5 * (ns[maxpos1] + ns[maxpos2]); + average.Normalize(); + new_normals[maxpos1] = average; + new_normals.DeleteElement(maxpos2); + auto gw = CalcGrowthVector(new_normals); - for(auto pi : points) - growthvectors[pi] += normals[pi]; - } + for (auto n : ns) + if (n * gw < 0) + throw SpecialPointException(); + return gw; +} +SpecialBoundaryPoint ::GrowthGroup ::GrowthGroup(FlatArray faces_, + FlatArray> normals) +{ + faces = faces_; + growth_vector = CalcGrowthVector(normals); +} - BoundaryLayerTool::BoundaryLayerTool(Mesh & mesh_, const BoundaryLayerParameters & params_) - : mesh(mesh_), topo(mesh_.GetTopology()), params(params_) - { - static Timer timer("BoundaryLayerTool::ctor"); - RegionTimer regt(timer); - ProcessParameters(); - - //for(auto & seg : mesh.LineSegments()) - //seg.edgenr = seg.epgeominfo[1].edgenr; - - height = 0.0; - for (auto h : par_heights) - height += h; - - max_edge_nr = -1; - for(const auto& seg : mesh.LineSegments()) - if(seg.edgenr > max_edge_nr) - max_edge_nr = seg.edgenr; - - int ndom = mesh.GetNDomains(); - ndom_old = ndom; - - new_mat_nrs.SetSize(mesh.FaceDescriptors().Size() + 1); - new_mat_nrs = -1; - for(auto [bcname, matname] : par_new_mat) - { - mesh.SetMaterial(++ndom, matname); - regex pattern(bcname); - for(auto i : Range(1, mesh.GetNFD()+1)) - { - auto& fd = mesh.GetFaceDescriptor(i); - if(regex_match(fd.GetBCName(), pattern)) - new_mat_nrs[i] = ndom; - } - } - - if(!params.outside) - domains.Invert(); - - topo.SetBuildVertex2Element(true); - mesh.UpdateTopology(); - - have_single_segments = HaveSingleSegments(mesh); - if(have_single_segments) - segments = BuildSegments(mesh); - else - segments = mesh.LineSegments(); - - np = mesh.GetNP(); - ne = mesh.GetNE(); - nse = mesh.GetNSE(); - nseg = segments.Size(); - - p2sel = mesh.CreatePoint2SurfaceElementTable(); - - nfd_old = mesh.GetNFD(); - moved_surfaces.SetSize(nfd_old+1); - moved_surfaces.Clear(); - si_map.SetSize(nfd_old+1); - for(auto i : Range(nfd_old+1)) - si_map[i] = i; - } - - void BoundaryLayerTool :: CreateNewFaceDescriptors() - { - surfacefacs.SetSize(nfd_old+1); - surfacefacs = 0.0; - // create new FaceDescriptors - for(auto i : Range(1, nfd_old+1)) - { - const auto& fd = mesh.GetFaceDescriptor(i); - string name = fd.GetBCName(); - if(par_surfid.Contains(i)) - { - if(auto isIn = domains.Test(fd.DomainIn()); isIn != domains.Test(fd.DomainOut())) - { - 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_nrs[i] : fd.DomainIn(), - isIn ? fd.DomainOut() : new_mat_nrs[i], -1); - new_fd.SetBCProperty(new_si); - new_fd.SetSurfColour(fd.SurfColour()); - mesh.AddFaceDescriptor(new_fd); - si_map[i] = new_si; - moved_surfaces.SetBit(i); - mesh.SetBCName(new_si-1, "mapped_" + name); - } - // curving of surfaces with boundary layers will often - // result in pushed through elements, since we do not (yet) - // curvature through layers. - // Therefore we disable curving for these surfaces. - if(!params.keep_surfaceindex) - mesh.GetFaceDescriptor(i).SetSurfNr(-1); - } - } - - for(auto si : par_surfid) - if(surfacefacs[si] == 0.0) - throw Exception("Surface " + to_string(si) + " is not a boundary of the domain to be grown into!"); - } - - void BoundaryLayerTool ::CreateFaceDescriptorsSides() - { - BitArray face_done(mesh.GetNFD()+1); - face_done.Clear(); - for(const auto& sel : mesh.SurfaceElements()) - { - auto facei = sel.GetIndex(); - if(face_done.Test(facei)) - continue; - bool point_moved = false; - // bool point_fixed = false; - for(auto pi : sel.PNums()) - { - if(growthvectors[pi].Length() > 0) - point_moved = true; - /* - else - point_fixed = true; - */ - } - if(point_moved && !moved_surfaces.Test(facei)) - { - int new_si = mesh.GetNFD()+1; - const auto& fd = mesh.GetFaceDescriptor(facei); - // auto isIn = domains.Test(fd.DomainIn()); - // auto isOut = domains.Test(fd.DomainOut()); - int si = params.sides_keep_surfaceindex ? facei : -1; - // domin and domout can only be set later - FaceDescriptor new_fd(si, -1, - -1, si); - new_fd.SetBCProperty(new_si); - mesh.AddFaceDescriptor(new_fd); - si_map[facei] = new_si; - mesh.SetBCName(new_si-1, fd.GetBCName()); - face_done.SetBit(facei); - } - } - } - - void BoundaryLayerTool :: CalculateGrowthVectors() - { - growthvectors.SetSize(np); - growthvectors = 0.; - - for(auto pi : mesh.Points().Range()) +SpecialBoundaryPoint ::SpecialBoundaryPoint( + const std::map>& normals) +{ + // find opposing face normals + Array> ns; + Array faces; + for (auto [face, normal] : normals) { - const auto & p = mesh[pi]; - if(p.Type() == INNERPOINT) + ns.Append(normal); + faces.Append(face); + } + + auto [minface1, minface2] = FindCloseVectors(ns, false); + minface1 = faces[minface1]; + minface2 = faces[minface2]; + Array g1_faces; + g1_faces.Append(minface1); + Array g2_faces; + g2_faces.Append(minface2); + auto n1 = normals.at(minface1); + auto n2 = normals.at(minface2); + separating_direction = 0.5 * (n2 - n1); + + Array> normals1, normals2; + for (auto [facei, normali] : normals) + if (facei != minface1 && facei != minface2) + { + g1_faces.Append(facei); + g2_faces.Append(facei); + } + for (auto fi : g1_faces) + normals1.Append(normals.at(fi)); + for (auto fi : g2_faces) + normals2.Append(normals.at(fi)); + growth_groups.Append(GrowthGroup(g1_faces, normals1)); + growth_groups.Append(GrowthGroup(g2_faces, normals2)); +} + +Vec<3> BoundaryLayerTool ::getEdgeTangent(PointIndex pi, int edgenr, FlatArray segs) +{ + Vec<3> tangent = 0.0; + ArrayMem pts; + for (auto* p_seg : segs) + { + auto& seg = *p_seg; + if (seg.edgenr != edgenr) + continue; + PointIndex other = seg[0] + seg[1] - pi; + if (!pts.Contains(other)) + pts.Append(other); + } + if (pts.Size() != 2) + { + cout << "getEdgeTangent pi = " << pi << ", edgenr = " << edgenr << endl; + cout << pts << endl; + for (auto* p_seg : segs) + cout << *p_seg << endl; + throw NG_EXCEPTION("Something went wrong in getEdgeTangent!"); + } + tangent = mesh[pts[1]] - mesh[pts[0]]; + return tangent.Normalize(); +} + +void BoundaryLayerTool ::LimitGrowthVectorLengths() +{ + static Timer tall("BoundaryLayerTool::LimitGrowthVectorLengths"); + RegionTimer rtall(tall); + + GrowthVectorLimiter limiter(*this); + limiter.Perform(); +} + +// depending on the geometry type, the mesh contains segments multiple times +// (once for each face) +bool HaveSingleSegments (const Mesh& mesh) +{ + auto& topo = mesh.GetTopology(); + NgArray surf_els; + + for (auto segi : Range(mesh.LineSegments())) + { + mesh.GetTopology().GetSegmentSurfaceElements(segi + 1, surf_els); + if (surf_els.Size() < 2) + continue; + + auto seg = mesh[segi]; + auto pi0 = min(seg[0], seg[1]); + auto pi1 = max(seg[0], seg[1]); + auto p0_segs = topo.GetVertexSegments(seg[0]); + + for (auto segi_other : p0_segs) + { + if (segi_other == segi) continue; - std::map> normals; - - // calculate one normal vector per face (average with angles as weights for multiple surface elements within a face) - for(auto sei : p2sel[pi]) - { - const auto & sel = mesh[sei]; - auto facei = sel.GetIndex(); - if(!par_surfid.Contains(facei)) - continue; - - auto n = surfacefacs[sel.GetIndex()] * getNormal(sel); - - int itrig = sel.PNums().Pos(pi); - itrig += sel.GetNP(); - auto v0 = (mesh[sel.PNumMod(itrig+1)] - mesh[pi]).Normalize(); - auto v1 = (mesh[sel.PNumMod(itrig-1)] - mesh[pi]).Normalize(); - if(normals.count(facei)==0) - normals[facei] = {0.,0.,0.}; - normals[facei] += acos(v0*v1)*n; + auto seg_other = mesh[segi_other]; + auto pi0_other = min(seg_other[0], seg_other[1]); + auto pi1_other = max(seg_other[0], seg_other[1]); + if (pi0_other == pi0 && pi1_other == pi1) + return false; } - for(auto & [facei, n] : normals) - n *= 1.0/n.Length(); + // found segment with multiple adjacent surface elements but no other + // segments with same points -> have single segments + return true; + } - // combine normal vectors for each face to keep uniform distances - auto & np = growthvectors[pi]; - ArrayMem, 3> ns; - for (auto &[facei, n] : normals) { + return true; +} + +// duplicates segments (and sets seg.si accordingly) to have a unified data +// structure for all geometry types +Array BuildSegments (Mesh& mesh) +{ + Array segments; + // auto& topo = mesh.GetTopology(); + + NgArray surf_els; + + for (auto segi : Range(mesh.LineSegments())) + { + auto seg = mesh[segi]; + mesh.GetTopology().GetSegmentSurfaceElements(segi + 1, surf_els); + for (auto seli : surf_els) + { + const auto& sel = mesh[seli]; + seg.si = sel.GetIndex(); + + auto np = sel.GetNP(); + for (auto i : Range(np)) + { + if (sel[i] == seg[0]) + { + if (sel[(i + 1) % np] != seg[1]) + swap(seg[0], seg[1]); + break; + } + } + + segments.Append(seg); + } + } + return segments; +} + +void MergeAndAddSegments (Mesh& mesh, FlatArray segments, FlatArray new_segments) +{ + INDEX_2_HASHTABLE already_added(segments.Size() + 2 * new_segments.Size()); + + mesh.LineSegments().SetSize0(); + + auto addSegment = [&] (const auto& seg) { + INDEX_2 i2(seg[0], seg[1]); + i2.Sort(); + if (!already_added.Used(i2)) + { + mesh.AddSegment(seg); + already_added.Set(i2, true); + } + }; + + for (const auto& seg : segments) + addSegment(seg); + + for (const auto& seg : new_segments) + addSegment(seg); +} + +BoundaryLayerTool::BoundaryLayerTool(Mesh& mesh_, + const BoundaryLayerParameters& params_) + : mesh(mesh_), topo(mesh_.GetTopology()), params(params_) +{ + static Timer timer("BoundaryLayerTool::ctor"); + RegionTimer regt(timer); + ProcessParameters(); + if (domains.NumSet() == 0) + return; + + topo.SetBuildVertex2Element(true); + mesh.UpdateTopology(); + + old_segments = mesh.LineSegments(); + have_single_segments = HaveSingleSegments(mesh); + + if (have_single_segments) + segments = BuildSegments(mesh); + else + segments = mesh.LineSegments(); + + np = mesh.GetNP(); + ne = mesh.GetNE(); + nse = mesh.GetNSE(); + nseg = segments.Size(); + + p2sel = mesh.CreatePoint2SurfaceElementTable(); + + nfd_old = mesh.GetNFD(); + moved_surfaces.SetSize(nfd_old + 1); + moved_surfaces.Clear(); + si_map.SetSize(nfd_old + 1); + for (auto i : Range(nfd_old + 1)) + si_map[i] = i; +} + +void BoundaryLayerTool ::CreateNewFaceDescriptors() +{ + surfacefacs.SetSize(nfd_old + 1); + surfacefacs = 0.0; + // create new FaceDescriptors + for (auto i : Range(1, nfd_old + 1)) + { + const auto& fd = mesh.GetFaceDescriptor(i); + string name = fd.GetBCName(); + if (par_surfid.Contains(i)) + { + if (auto isIn = domains.Test(fd.DomainIn()); + isIn != domains.Test(fd.DomainOut())) + { + int new_si = mesh.GetNFD() + 1; + surfacefacs[i] = isIn ? 1. : -1.; + moved_surfaces.SetBit(i); + if (!insert_only_volume_elements) + { + // -1 surf nr is so that curving does not do anything + FaceDescriptor new_fd(-1, isIn ? new_mat_nrs[i] : fd.DomainIn(), isIn ? fd.DomainOut() : new_mat_nrs[i], -1); + new_fd.SetBCProperty(new_si); + new_fd.SetSurfColour(fd.SurfColour()); + mesh.AddFaceDescriptor(new_fd); + si_map[i] = new_si; + mesh.SetBCName(new_si - 1, "mapped_" + name); + } + } + // curving of surfaces with boundary layers will often + // result in pushed through elements, since we do not (yet) + // curvature through layers. + // Therefore we disable curving for these surfaces. + if (params.disable_curving) + mesh.GetFaceDescriptor(i).SetSurfNr(-1); + } + } + + for (auto si : par_surfid) + if (surfacefacs[si] == 0.0) + throw Exception("Surface " + to_string(si) + " is not a boundary of the domain to be grown into!"); +} + +void BoundaryLayerTool ::CreateFaceDescriptorsSides() +{ + if (insert_only_volume_elements) + return; + BitArray face_done(mesh.GetNFD() + 1); + face_done.Clear(); + for (const auto& sel : mesh.SurfaceElements()) + { + auto facei = sel.GetIndex(); + if (face_done.Test(facei)) + continue; + bool point_moved = false; + // bool point_fixed = false; + for (auto pi : sel.PNums()) + { + if (growthvectors[pi].Length() > 0) + point_moved = true; + /* + else + point_fixed = true; + */ + } + if (point_moved && !moved_surfaces.Test(facei)) + { + int new_si = mesh.GetNFD() + 1; + const auto& fd = mesh.GetFaceDescriptor(facei); + // auto isIn = domains.Test(fd.DomainIn()); + // auto isOut = domains.Test(fd.DomainOut()); + int si = params.sides_keep_surfaceindex ? facei : -1; + // domin and domout can only be set later + FaceDescriptor new_fd(si, -1, -1, si); + new_fd.SetBCProperty(new_si); + mesh.AddFaceDescriptor(new_fd); + si_map[facei] = new_si; + mesh.SetBCName(new_si - 1, fd.GetBCName()); + face_done.SetBit(facei); + } + } +} + +void BoundaryLayerTool ::CalculateGrowthVectors() +{ + growthvectors.SetSize(np); + growthvectors = 0.; + + for (auto pi : mesh.Points().Range()) + { + const auto& p = mesh[pi]; + if (p.Type() == INNERPOINT) + continue; + + std::map> normals; + + // calculate one normal vector per face (average with angles as weights for + // multiple surface elements within a face) + for (auto sei : p2sel[pi]) + { + const auto& sel = mesh[sei]; + auto facei = sel.GetIndex(); + if (!par_surfid.Contains(facei)) + continue; + + auto n = surfacefacs[sel.GetIndex()] * getNormal(sel); + + int itrig = sel.PNums().Pos(pi); + itrig += sel.GetNP(); + auto v0 = (mesh[sel.PNumMod(itrig + 1)] - mesh[pi]).Normalize(); + auto v1 = (mesh[sel.PNumMod(itrig - 1)] - mesh[pi]).Normalize(); + if (normals.count(facei) == 0) + normals[facei] = {0., 0., 0.}; + normals[facei] += acos(v0 * v1) * n; + } + + for (auto& [facei, n] : normals) + n *= 1.0 / n.Length(); + + // combine normal vectors for each face to keep uniform distances + ArrayMem, 5> ns; + for (auto& [facei, n] : normals) + { ns.Append(n); } - ArrayMem, 3> removed; - // reduce to full rank of max 3 - while(true) - { - if(ns.Size() <= 1) - break; - if(ns.Size() == 2 && ns[0] * ns[1] < 1 - 1e-6) - break; - if (ns.Size() == 3) - { - DenseMatrix mat(3,3); - for(auto i : Range(3)) - for(auto j : Range(3)) - mat(i,j) = ns[i][j]; - if(fabs(mat.Det()) > 1e-6) - break; - } - int maxpos1 = 0; - int maxpos2 = 1; - double val = ns[0] * ns[1]; - for (auto i : Range(ns)) - { - for (auto j : Range(i + 1, ns.Size())) - { - double ip = ns[i] * ns[j]; - if(ip > val) - { - val = ip; - maxpos1 = i; - maxpos2 = j; - } - } - } - removed.Append(ns[maxpos1]); - removed.Append(ns[maxpos2]); - ns[maxpos1] = 0.5 * (ns[maxpos1] + ns[maxpos2]); - ns.DeleteElement(maxpos2); - } - - if(ns.Size() == 0) - continue; - if(ns.Size() == 1) - np = ns[0]; - else if(ns.Size() == 2) - { - np = ns[0]; - auto n = ns[1]; - 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); - } - else // ns.Size() == 3 - { - DenseMatrix mat(3,3); - for(auto i : Range(3)) - for(auto j : Range(3)) - mat(i, j) = ns[i] * ns[j]; - Vector rhs(3); - rhs = 1.; - Vector res(3); - DenseMatrix inv(3, ns.Size()); - CalcInverse(mat, inv); - inv.Mult(rhs, res); - for(auto i : Range(ns)) - np += res[i] * ns[i]; - } - for(auto& n : removed) - if(n * np < 0) - cout << "WARNING: Growth vector at point " << pi << " in opposite direction to face normal!" << endl << "Growthvector = " << np << ", face normal = " << n << endl; - } - } - - Array>, SegmentIndex> BoundaryLayerTool :: BuildSegMap() - { - // Bit array to keep track of segments already processed - BitArray segs_done(nseg+1); - segs_done.Clear(); - - // 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(segments.Size()); - - // moved segments - is_edge_moved.SetSize(max_edge_nr+1); - is_edge_moved = false; - - // boundaries to project endings to - is_boundary_projected.SetSize(nfd_old+1); - is_boundary_projected.Clear(); - is_boundary_moved.SetSize(nfd_old+1); - is_boundary_moved.Clear(); - - for(auto si : Range(segments)) - { - if(segs_done[si]) continue; - const auto& segi = segments[si]; - if(!moved_surfaces.Test(segi.si)) continue; - segs_done.SetBit(si); - segmap[si].Append(make_pair(si, 0)); - moved_segs.Append(si); - is_edge_moved.SetBit(segi.edgenr); - for(auto sj : Range(segments)) - { - if(segs_done.Test(sj)) continue; - const auto& segj = segments[sj]; - if((segi[0] == segj[0] && segi[1] == segj[1]) || - (segi[0] == segj[1] && segi[1] == segj[0])) - { - segs_done.SetBit(sj); - int type; - if(moved_surfaces.Test(segj.si)) - type = 0; - else if(const auto& fd = mesh.GetFaceDescriptor(segj.si); domains.Test(fd.DomainIn()) && domains.Test(fd.DomainOut())) - { - type = 2; - if(fd.DomainIn() == 0 || fd.DomainOut() == 0) - is_boundary_projected.SetBit(segj.si); - } - else if(const auto& fd = mesh.GetFaceDescriptor(segj.si); !domains.Test(fd.DomainIn()) && !domains.Test(fd.DomainOut())) - { - type = 3; - is_boundary_moved.SetBit(segj.si); - } - else - { - type = 1; - // in case 1 we project the growthvector onto the surface - is_boundary_projected.SetBit(segj.si); - } - segmap[si].Append(make_pair(sj, type)); - } - } - } - - return segmap; - } - - BitArray BoundaryLayerTool :: ProjectGrowthVectorsOnSurface() - { - BitArray in_surface_direction(nfd_old+1); - in_surface_direction.Clear(); - // project growthvector on surface for inner angles - if(params.grow_edges) - { - for(const auto& sel : mesh.SurfaceElements()) - if(is_boundary_projected.Test(sel.GetIndex())) - { - auto n = getNormal(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(); - auto tol = v1.Length() * 1e-12; - if((v1 * v3 > -tol) && (v2 * v3 > -tol)) - in_surface_direction.SetBit(sel.GetIndex()); - else - continue; - - if(!par_project_boundaries.Contains(sel.GetIndex())) - 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; - } - } - } - else - { - for(const auto& seg : segments) - { - int count = 0; - for(const auto& seg2 : segments) - if(((seg[0] == seg2[0] && seg[1] == seg2[1]) || (seg[0] == seg2[1] && seg[1] == seg2[0])) && par_surfid.Contains(seg2.si)) - count++; - if(count == 1) - { - growthvectors[seg[0]] = {0., 0., 0.}; - growthvectors[seg[1]] = {0., 0., 0.}; - } - } - } - - return in_surface_direction; - } - - void BoundaryLayerTool :: InterpolateGrowthVectors() - { - // interpolate tangential component of growth vector along edge - for(auto edgenr : Range(max_edge_nr)) - { - // if(!is_edge_moved[edgenr+1]) continue; - - // build sorted list of edge - Array points; - // find first vertex on edge - double edge_len = 0.; - auto is_end_point = [&] (PointIndex pi) + try { - // if(mesh[pi].Type() == FIXEDPOINT) - // return true; - // return false; - auto segs = topo.GetVertexSegments(pi); - if(segs.Size() == 1) - return true; - auto first_edgenr = mesh[segs[0]].edgenr; - for(auto segi : segs) - if(mesh[segi].edgenr != first_edgenr) - return true; - return false; - }; + growthvectors[pi] = CalcGrowthVector(ns); + } + catch (const SpecialPointException& e) + { + special_boundary_points.emplace(pi, normals); + growthvectors[pi] = + special_boundary_points[pi].growth_groups[0].growth_vector; + } + } +} - bool any_grows = false; - - for(const auto& seg : segments) +Array>, SegmentIndex> +BoundaryLayerTool ::BuildSegMap() +{ + // Bit array to keep track of segments already processed + BitArray segs_done(nseg + 1); + segs_done.Clear(); + + // 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(segments.Size()); + + // moved segments + is_edge_moved.SetSize(max_edge_nr + 1); + is_edge_moved = false; + + // boundaries to project endings to + is_boundary_projected.SetSize(nfd_old + 1); + is_boundary_projected.Clear(); + is_boundary_moved.SetSize(nfd_old + 1); + is_boundary_moved.Clear(); + + for (auto si : Range(segments)) + { + if (segs_done[si]) + continue; + const auto& segi = segments[si]; + if (!moved_surfaces.Test(segi.si)) + continue; + segs_done.SetBit(si); + segmap[si].Append(make_pair(si, 0)); + moved_segs.Append(si); + is_edge_moved.SetBit(segi.edgenr); + for (auto sj : Range(segments)) + { + if (segs_done.Test(sj)) + continue; + const auto& segj = segments[sj]; + if ((segi[0] == segj[0] && segi[1] == segj[1]) || (segi[0] == segj[1] && segi[1] == segj[0])) + { + segs_done.SetBit(sj); + int type; + if (moved_surfaces.Test(segj.si)) + { + type = 0; + moved_segs.Append(sj); + } + else if (const auto& fd = mesh.GetFaceDescriptor(segj.si); + domains.Test(fd.DomainIn()) && domains.Test(fd.DomainOut())) + { + type = 2; + if (fd.DomainIn() == 0 || fd.DomainOut() == 0) + is_boundary_projected.SetBit(segj.si); + } + else if (const auto& fd = mesh.GetFaceDescriptor(segj.si); + !domains.Test(fd.DomainIn()) && !domains.Test(fd.DomainOut())) + { + type = 3; + // cout << "set is_moved boundary to type 3 for " << segj.si << endl; + is_boundary_moved.SetBit(segj.si); + } + else + { + type = 1; + // in case 1 we project the growthvector onto the surface + is_boundary_projected.SetBit(segj.si); + } + segmap[si].Append(make_pair(sj, type)); + } + } + } + + return segmap; +} + +BitArray BoundaryLayerTool ::ProjectGrowthVectorsOnSurface() +{ + BitArray in_surface_direction(nfd_old + 1); + in_surface_direction.Clear(); + // project growthvector on surface for inner angles + if (params.grow_edges) + { + for (const auto& sel : mesh.SurfaceElements()) + if (is_boundary_projected.Test(sel.GetIndex())) { - if(seg.edgenr-1 == edgenr) + auto n = getNormal(sel); + for (auto i : Range(sel.PNums())) { - if(growthvectors[seg[0]].Length2() != 0 || - growthvectors[seg[1]].Length2() != 0) - any_grows = true; - if(points.Size() == 0 && is_end_point(seg[0])) - { - points.Append(seg[0]); - points.Append(seg[1]); - edge_len += (mesh[seg[1]] - mesh[seg[0]]).Length(); - } + 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(); + auto tol = v1.Length() * 1e-12; + if ((v1 * v3 > -tol) && (v2 * v3 > -tol)) + in_surface_direction.SetBit(sel.GetIndex()); + else + continue; + + if (!par_project_boundaries.Contains(sel.GetIndex())) + 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; } } + } + else + { + for (const auto& seg : segments) + { + int count = 0; + for (const auto& seg2 : segments) + if (((seg[0] == seg2[0] && seg[1] == seg2[1]) || (seg[0] == seg2[1] && seg[1] == seg2[0])) && par_surfid.Contains(seg2.si)) + count++; + if (count == 1) + { + growthvectors[seg[0]] = {0., 0., 0.}; + growthvectors[seg[1]] = {0., 0., 0.}; + } + } + } - if(!any_grows) - continue; + return in_surface_direction; +} - if(!points.Size()) - throw Exception("Could not find startpoint for edge " + ToString(edgenr)); +void BoundaryLayerTool ::InsertNewElements( + FlatArray>, SegmentIndex> segmap, + const BitArray& in_surface_direction) +{ + static Timer timer("BoundaryLayerTool::InsertNewElements"); + RegionTimer rt(timer); + mapto.SetSize(0); + mapto.SetSize(np); + mapfrom.SetSize(mesh.GetNP()); + mapfrom = PointIndex::INVALID; - while(true) - { - bool point_found = false; - for(auto si : topo.GetVertexSegments(points.Last())) - { - const auto& seg = mesh[si]; - if(seg.edgenr-1 != edgenr) - continue; - if(seg[0] == points.Last() && points[points.Size()-2] !=seg[1]) - { - edge_len += (mesh[points.Last()] - mesh[seg[1]]).Length(); - points.Append(seg[1]); - point_found = true; - break; - } - else if(seg[1] == points.Last() && - points[points.Size()-2] != seg[0]) - { - edge_len += (mesh[points.Last()] - mesh[seg[0]]).Length(); - points.Append(seg[0]); - point_found = true; - break; - } - } - if(is_end_point(points.Last())) - break; - if(!point_found) - { - throw Exception(string("Could not find connected list of line segments for edge ") + edgenr); - } - } + auto changed_domains = domains; + if (!params.outside) + changed_domains.Invert(); - if(growthvectors[points[0]].Length2() == 0 && - growthvectors[points.Last()].Length2() == 0) - continue; + auto& identifications = mesh.GetIdentifications(); + const int identnr = identifications.GetNr("boundarylayer"); - // tangential part of growth vectors - auto t1 = (mesh[points[1]]-mesh[points[0]]).Normalize(); - auto gt1 = growthvectors[points[0]] * t1 * t1; - auto t2 = (mesh[points.Last()]-mesh[points[points.Size()-2]]).Normalize(); - auto gt2 = growthvectors[points.Last()] * t2 * t2; - - if(!is_edge_moved[edgenr+1]) - { - if(growthvectors[points[0]] * (mesh[points[1]] - mesh[points[0]]) < 0) - gt1 = 0.; - if(growthvectors[points.Last()] * (mesh[points[points.Size()-2]] - mesh[points.Last()]) < 0) - gt2 = 0.; - } - - double len = 0.; - for(size_t i = 1; i < points.Size()-1; i++) - { - auto pi = points[i]; - len += (mesh[pi] - mesh[points[i-1]]).Length(); - auto t = getEdgeTangent(pi, edgenr); - auto lam = len/edge_len; - auto interpol = (1-lam) * (gt1 * t) * t + lam * (gt2 * t) * t; - growthvectors[pi] += interpol; - } + auto add_points = [&] (PointIndex pi, Vec<3>& growth_vector, Array& new_points) { + Point<3> p = mesh[pi]; + PointIndex pi_last = pi; + double height = 0.0; + for (auto i : Range(par_heights)) + { + height += par_heights[i]; + auto pi_new = mesh.AddPoint(p); + // mesh.AddLockedPoint(pi_new); + mapfrom.Append(pi); + new_points.Append(pi_new); + growth_vector_map[pi_new] = {&growth_vector, height}; + // if (special_boundary_points.count(pi) > 0) + // mesh.AddLockedPoint(pi_new); + pi_last = pi_new; } + }; - InterpolateSurfaceGrowthVectors(); - } - - void BoundaryLayerTool :: InsertNewElements( FlatArray>, SegmentIndex> segmap, const BitArray & in_surface_direction ) - { - static Timer timer("BoundaryLayerTool::InsertNewElements"); RegionTimer rt(timer); - Array, PointIndex> mapto(np); - // insert new points - for (PointIndex pi = 1; pi <= np; pi++) + // insert new points + for (PointIndex pi = 1; pi <= np; pi++) + { if (growthvectors[pi].Length2() != 0) { - Point<3> p = mesh[pi]; - for(auto i : Range(par_heights)) + if (special_boundary_points.count(pi)) { - p += par_heights[i] * growthvectors[pi]; - mapto[pi].Append(mesh.AddPoint(p)); + for (auto& group : special_boundary_points[pi].growth_groups) + add_points(pi, group.growth_vector, group.new_points); } + else + add_points(pi, growthvectors[pi], mapto[pi]); } + } - // add 2d quads on required surfaces - map, int> seg2edge; - if(params.grow_edges) + // get point from mapto (or the group if point is mapped to multiple new + // points) layer = -1 means last point (top of boundary layer) + auto newPoint = [&] (PointIndex pi, int layer = -1, int group = 0) { + if (layer == -1) + layer = par_heights.Size() - 1; + if (special_boundary_points.count(pi)) + return special_boundary_points[pi].growth_groups[group].new_points[layer]; + else + return mapto[pi][layer]; + }; + + auto hasMoved = [&] (PointIndex pi) { + return mapto[pi].Size() > 0 || special_boundary_points.count(pi); + }; + + auto numGroups = [&] (PointIndex pi) -> size_t { + if (special_boundary_points.count(pi)) + return special_boundary_points[pi].growth_groups.Size(); + else + return 1; + }; + + auto getGroups = [&] (PointIndex pi, int face_index) -> Array { + auto n = numGroups(pi); + Array groups; + if (n == 1) { - for(auto sei : moved_segs) - { - // copy here since we will add segments and this would - // invalidate a reference! - // auto segi = segments[sei]; - for(auto [sej, type] : segmap[sei]) - { - auto segj = segments[sej]; - if(type == 0) - { + groups.Append(0); + return groups; + } + const auto& all_groups = special_boundary_points[pi].growth_groups; + for (auto i : Range(n)) + if (all_groups[i].faces.Contains(face_index)) + groups.Append(i); + // cout << "groups " << pi << ", " << face_index << endl << groups; + return groups; + }; + + // add 2d quads on required surfaces + map, int> seg2edge; + map edge_map; + int edge_nr = max_edge_nr; + auto getEdgeNr = [&] (int ei) { + if (edge_map.count(ei) == 0) + edge_map[ei] = ++edge_nr; + return edge_map[ei]; + }; + if (params.grow_edges) + { + for (auto sei : moved_segs) + { + // copy here since we will add segments and this would + // invalidate a reference! + // auto segi = segments[sei]; + for (auto [sej, type] : segmap[sei]) + { + auto segj = segments[sej]; + if (type == 0) + { + auto addSegment = [&] (PointIndex p0, PointIndex p1, bool extra_edge_nr = false) { Segment s; - s[0] = mapto[segj[0]].Last(); - s[1] = mapto[segj[1]].Last(); + s[0] = p0; + s[1] = p1; 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]; + auto pair = + s[0] < s[1] ? make_pair(s[0], s[1]) : make_pair(s[1], s[0]); + if (extra_edge_nr) + s.edgenr = ++edge_nr; + else + s.edgenr = getEdgeNr(segj.edgenr); s.si = si_map[segj.si]; new_segments.Append(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); - is_boundary_moved.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; - new_segments.Append(s0); + // cout << __LINE__ <<"\t" << s << endl; + return s; + }; - for(auto i : Range(par_heights)) - { - Element2d sel(QUAD); - p3 = mapto[pp2][i]; - p4 = mapto[pp1][i]; - sel[0] = p1; - sel[1] = p2; - sel[2] = p3; - sel[3] = p4; - for(auto i : Range(4)) + auto p0 = segj[0], p1 = segj[1]; + auto g0 = getGroups(p0, segj.si); + auto g1 = getGroups(p1, segj.si); + + if (g0.Size() == 1 && g1.Size() == 1) + auto s = + addSegment(newPoint(p0, -1, g0[0]), newPoint(p1, -1, g1[0])); + else + { + if (g0.Size() == 2) + addSegment(newPoint(p0, -1, g0[0]), newPoint(p0, -1, g0[1])); + if (g1.Size() == 2) + addSegment(newPoint(p1, -1, g1[0]), newPoint(p1, -1, g1[1])); + } + } + // 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); + is_boundary_moved.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; + new_segments.Append(s0); + if (type == 3) + new_segments_on_moved_bnd.Append(s0); + + for (auto i : Range(par_heights)) + { + Element2d sel(QUAD); + p3 = newPoint(pp2, i); + p4 = newPoint(pp1, i); + sel[0] = p1; + sel[1] = p2; + sel[2] = p3; + sel[3] = p4; + for (auto i : Range(4)) { - sel.GeomInfo()[i].u = 0.0; - sel.GeomInfo()[i].v = 0.0; + sel.GeomInfo()[i].u = 0.0; + sel.GeomInfo()[i].v = 0.0; } - sel.SetIndex(si_map[segj.si]); - mesh.AddSurfaceElement(sel); + identifications.Add(p1, p4, identnr); + identifications.Add(p2, p3, identnr); + sel.SetIndex(si_map[segj.si]); + new_sels.Append(sel); + new_sels_on_moved_bnd.Append(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; - new_segments.Append(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; - new_segments.Append(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; - new_segments.Append(s3); - } - } - } - } - - BitArray fixed_points(np+1); - fixed_points.Clear(); - BitArray moveboundarypoint(np+1); - moveboundarypoint.Clear(); - for(SurfaceElementIndex si = 0; si < nse; si++) - { - // copy because surfaceels array will be resized! - auto sel = mesh[si]; - if(moved_surfaces.Test(sel.GetIndex())) - { - Array points(sel.PNums()); - if(surfacefacs[sel.GetIndex()] > 0) Swap(points[0], points[2]); - for(auto j : Range(par_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]; - auto new_index = new_mat_nrs[sel.GetIndex()]; - if(new_index == -1) - throw Exception("Boundary " + ToString(sel.GetIndex()) + " with name " + mesh.GetBCName(sel.GetIndex()-1) + " extruded, but no new material specified for it!"); - el.SetIndex(new_mat_nrs[sel.GetIndex()]); - mesh.AddVolumeElement(el); - } - Element2d newel = sel; - for(auto& p : newel.PNums()) - p = mapto[p].Last(); - newel.SetIndex(si_map[sel.GetIndex()]); - mesh.AddSurfaceElement(newel); - } - else - { - bool has_moved = false; - for(auto p : sel.PNums()) - if(mapto[p].Size()) - has_moved = true; - if(has_moved) - for(auto p : sel.PNums()) + // 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); + s1.edgenr = getEdgeNr(segj.edgenr); + s1.si = segj.si; + // new_segments.Append(s1); + Segment s2; + s2[0] = p4; + s2[1] = p1; + s2[2] = PointIndex::INVALID; + pair = make_pair(p1, p4); + s2.edgenr = getEdgeNr(segj.edgenr); + s2.si = segj.si; + // new_segments.Append(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); + s3.edgenr = getEdgeNr(segj.edgenr); + s3.si = segj.si; + new_segments.Append(s3); + if (type == 3) + new_segments_on_moved_bnd.Append(s0); + } + else if (type == 3) { - if(!mapto[p].Size()) + PointIndex pp1 = segj[1]; + PointIndex pp2 = segj[0]; + if (!in_surface_direction.Test(segj.si)) { - fixed_points.SetBit(p); - if(is_boundary_moved.Test(sel.GetIndex())) - moveboundarypoint.SetBit(p); + Swap(pp1, pp2); + } + PointIndex p1 = pp1; + PointIndex p2 = pp2; + PointIndex p3, p4; + + for (auto i : Range(par_heights)) + { + Element2d sel(QUAD); + p3 = newPoint(pp2, i); + p4 = newPoint(pp1, i); + sel[0] = p1; + sel[1] = p2; + sel[2] = p3; + sel[3] = p4; + for (auto i : Range(4)) + { + sel.GeomInfo()[i].u = 0.0; + sel.GeomInfo()[i].v = 0.0; + } + identifications.Add(p1, p4, identnr); + identifications.Add(p2, p3, identnr); + sel.SetIndex(si_map[segj.si]); + new_sels.Append(sel); + new_sels_on_moved_bnd.Append(sel); + p1 = p4; + p2 = p3; } } - } - if(is_boundary_moved.Test(sel.GetIndex())) - { - for(auto& p : mesh[si].PNums()) - if(mapto[p].Size()) - p = mapto[p].Last(); - } - } - - for(SegmentIndex sei = 0; sei < nseg; sei++) - { - auto& seg = segments[sei]; - if(is_boundary_moved.Test(seg.si)) - for(auto& p : seg.PNums()) - if(mapto[p].Size()) - p = mapto[p].Last(); - } - - for(ElementIndex ei = 0; ei < ne; ei++) - { - auto el = mesh[ei]; - ArrayMem fixed; - ArrayMem moved; - bool moved_bnd = false; - for(const auto& p : el.PNums()) - { - if(fixed_points.Test(p)) - fixed.Append(p); - if(mapto[p].Size()) - moved.Append(p); - if(moveboundarypoint.Test(p)) - moved_bnd = true; - } - - bool do_move, do_insert; - if(domains.Test(el.GetIndex())) - { - do_move = fixed.Size() && moved_bnd; - do_insert = do_move; - } - else - { - do_move = !fixed.Size() || moved_bnd; - do_insert = !do_move; - } - - if(do_move) - { - for(auto& p : mesh[ei].PNums()) - if(mapto[p].Size()) - p = mapto[p].Last(); - } - if(do_insert) - { - if(el.GetType() == TET) - { - if(moved.Size() == 3) // inner corner - { - PointIndex p1 = moved[0]; - PointIndex p2 = moved[1]; - PointIndex p3 = moved[2]; - auto v1 = mesh[p1]; - auto n = Cross(mesh[p2]-v1, mesh[p3]-v1); - auto d = mesh[mapto[p1][0]] - v1; - if(n*d > 0) - Swap(p2,p3); - PointIndex p4 = p1; - PointIndex p5 = p2; - PointIndex p6 = p3; - for(auto i : Range(par_heights)) - { - Element nel(PRISM); - nel[0] = p4; nel[1] = p5; nel[2] = p6; - p4 = mapto[p1][i]; p5 = mapto[p2][i]; p6 = mapto[p3][i]; - nel[3] = p4; nel[4] = p5; nel[5] = p6; - nel.SetIndex(el.GetIndex()); - mesh.AddVolumeElement(nel); - } - } - if(moved.Size() == 2) - { - if(fixed.Size() == 1) - { - PointIndex p1 = moved[0]; - PointIndex p2 = moved[1]; - for(auto i : Range(par_heights)) - { - PointIndex p3 = mapto[moved[1]][i]; - PointIndex p4 = mapto[moved[0]][i]; - Element nel(PYRAMID); - nel[0] = p1; - nel[1] = p2; - nel[2] = p3; - nel[3] = p4; - // nel[4] = el[0] + el[1] + el[2] + el[3] - fixed[0] - moved[0] - moved[1]; - nel[4] = el[0]-fixed[0] + el[1]-moved[0] + el[2]-moved[1] + el[3]; - if(Cross(mesh[p2]-mesh[p1], mesh[p4]-mesh[p1]) * (mesh[nel[4]]-mesh[nel[1]]) > 0) - Swap(nel[1], nel[3]); - nel.SetIndex(el.GetIndex()); - mesh.AddVolumeElement(nel); - p1 = p4; - p2 = p3; - } - } - } - if(moved.Size() == 1 && fixed.Size() == 1) - { - PointIndex p1 = moved[0]; - for(auto i : Range(par_heights)) - { - Element nel = el; - PointIndex p2 = mapto[moved[0]][i]; - for(auto& p : nel.PNums()) - { - if(p == moved[0]) - p = p1; - else if(p == fixed[0]) - p = p2; - } - p1 = p2; - mesh.AddVolumeElement(nel); - } - } - } - else if(el.GetType() == PYRAMID) - { - if(moved.Size() == 2) - { - if(fixed.Size() != 2) - throw Exception("This case is not implemented yet! Fixed size = " + ToString(fixed.Size())); - PointIndex p1 = moved[0]; - PointIndex p2 = moved[1]; - for(auto i : Range(par_heights)) - { - PointIndex p3 = mapto[moved[1]][i]; - PointIndex p4 = mapto[moved[0]][i]; - Element nel(PYRAMID); - nel[0] = p1; - nel[1] = p2; - nel[2] = p3; - nel[3] = p4; - // nel[4] = el[0] + el[1] + el[2] + el[3] + el[4] - fixed[0] - fixed[1] - moved[0] - moved[1]; - nel[4] = el[0]-fixed[0] + el[1]-fixed[1] + el[2]-moved[0] + el[3]-moved[1] + el[4]; - if(Cross(mesh[p2] - mesh[p1], mesh[p4]-mesh[p1]) * (mesh[nel[4]]-mesh[nel[1]]) > 0) - Swap(nel[1], nel[3]); - nel.SetIndex(el.GetIndex()); - mesh.AddVolumeElement(nel); - p1 = p4; - p2 = p3; - } - } - else if(moved.Size() == 1) - throw Exception("This case is not implemented yet!"); - } - else - throw Exception("Boundarylayer only implemented for tets and pyramids outside yet!"); - } - } - } - - void BoundaryLayerTool :: SetDomInOut() - { - for(auto i : Range(1, nfd_old+1)) - if(moved_surfaces.Test(i)) - { - if(auto dom = mesh.GetFaceDescriptor(si_map[i]).DomainIn(); dom > ndom_old) - mesh.GetFaceDescriptor(i).SetDomainOut(dom); - else - mesh.GetFaceDescriptor(i).SetDomainIn(mesh.GetFaceDescriptor(si_map[i]).DomainOut()); + } } - } - - void BoundaryLayerTool :: SetDomInOutSides() - { - BitArray done(mesh.GetNFD()+1); - done.Clear(); - for(auto sei : Range(mesh.SurfaceElements())) - { - auto& sel = mesh[sei]; - auto index = sel.GetIndex(); - if(done.Test(index)) - continue; - done.SetBit(index); - auto& fd = mesh.GetFaceDescriptor(index); - if(fd.DomainIn() != -1) - continue; - int e1, e2; - mesh.GetTopology().GetSurface2VolumeElement(sei+1, e1, e2); - if(e1 == 0) - fd.SetDomainIn(0); - else - fd.SetDomainIn(mesh.VolumeElement(e1).GetIndex()); - if(e2 == 0) - fd.SetDomainOut(0); - else - fd.SetDomainOut(mesh.VolumeElement(e2).GetIndex()); - } - } - - void BoundaryLayerTool :: AddSegments() - { - if(have_single_segments) - MergeAndAddSegments(mesh, new_segments); - else - { - for(auto & seg : new_segments) - mesh.AddSegment(seg); } - } - void BoundaryLayerTool :: FixVolumeElements() - { - static Timer timer("BoundaryLayerTool::FixVolumeElements"); RegionTimer rt(timer); - BitArray is_inner_point(mesh.GetNP()+1); - is_inner_point.Clear(); + auto getClosestGroup = [&] (PointIndex pi, SurfaceElementIndex sei) { + auto n = numGroups(pi); + if (n == 1) + return 0; + const auto& sel = mesh[sei]; + auto groups = getGroups(pi, sel.GetIndex()); + if (groups.Size() == 1) + return groups[0]; - auto changed_domains = domains; - if(!params.outside) - changed_domains.Invert(); + auto& growth_groups = special_boundary_points[pi].growth_groups; - for(ElementIndex ei : Range(ne)) - if(changed_domains.Test(mesh[ei].GetIndex())) - for(auto pi : mesh[ei].PNums()) - if(mesh[pi].Type() == INNERPOINT) - is_inner_point.SetBit(pi); + auto vdir = Center(mesh[sel[0]], mesh[sel[1]], mesh[sel[2]]) - mesh[pi]; + auto dot = vdir * special_boundary_points[pi].separating_direction; - Array points; - for(auto pi : mesh.Points().Range()) - if(is_inner_point.Test(pi)) - points.Append(pi); + return dot > 0 ? 1 : 0; + }; - auto p2el = mesh.CreatePoint2ElementTable(is_inner_point); - - // smooth growth vectors to shift additional element layers to the inside and fix flipped tets - for([[maybe_unused]] auto step : Range(10)) - { - for(auto pi : points) - { - Vec<3> average_gw = 0.0; - auto & els = p2el[pi]; - size_t cnt = 0; - for(auto ei : els) - if(ei(¶ms.boundary); bc) + BitArray fixed_points(np + 1); + fixed_points.Clear(); + auto p2el = mesh.CreatePoint2ElementTable(); + for (SurfaceElementIndex si = 0; si < nse; si++) + { + // copy because surfaceels array will be resized! + const auto sel = mesh[si]; + if (moved_surfaces.Test(sel.GetIndex())) { - for (int i = 1; i <= mesh.GetNFD(); i++) - if(mesh.GetFaceDescriptor(i).BCProperty() == *bc) - par_surfid.Append(i); - } - else if(string* s = get_if(¶ms.boundary); s) - { - regex pattern(*s); - BitArray boundaries(mesh.GetNFD()+1); - boundaries.Clear(); - for(int i = 1; i<=mesh.GetNFD(); i++) + Array points(sel.PNums()); + if (surfacefacs[sel.GetIndex()] > 0) + Swap(points[0], points[2]); + ArrayMem groups(points.Size()); + for (auto i : Range(points)) + groups[i] = getClosestGroup(sel[i], si); + bool add_volume_element = true; + for (auto pi : sel.PNums()) + if (numGroups(pi) > 1) + add_volume_element = false; + for (auto j : Range(par_heights)) { - auto& fd = mesh.GetFaceDescriptor(i); - if(regex_match(fd.GetBCName(), pattern)) + 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] = newPoint(sel.PNums()[i], j, groups[i]); + if (surfacefacs[sel.GetIndex()] > 0) + Swap(points[0], points[2]); + for (auto i : Range(points)) + el[sel.PNums().Size() + i] = points[i]; + auto new_index = new_mat_nrs[sel.GetIndex()]; + if (new_index == -1) + throw Exception("Boundary " + ToString(sel.GetIndex()) + " with name " + mesh.GetBCName(sel.GetIndex() - 1) + " extruded, but no new material specified for it!"); + el.SetIndex(new_mat_nrs[sel.GetIndex()]); + if (add_volume_element) + mesh.AddVolumeElement(el); + else { - boundaries.SetBit(i); - auto dom_pattern = get_if(¶ms.domain); - // only add if adjacent to domain - if(dom_pattern) - { - regex pattern(*dom_pattern); - bool mat1_match = fd.DomainIn() > 0 && regex_match(mesh.GetMaterial(fd.DomainIn()), pattern); - bool mat2_match = fd.DomainOut() > 0 && regex_match(mesh.GetMaterial(fd.DomainOut()), pattern); - // if boundary is inner or outer remove from list - if(mat1_match == mat2_match) - boundaries.Clear(i); - // if((fd.DomainIn() > 0 && regex_match(mesh.GetMaterial(fd.DomainIn()), pattern)) || (fd.DomainOut() > 0 && regex_match(self.GetMaterial(fd.DomainOut()), pattern))) - // boundaries.Clear(i); - // par_surfid.Append(i); - } - // else - // par_surfid.Append(i); + // Let the volume mesher fill the hole with pyramids/tets + // To insert pyramids, we need close surface identifications on open + // quads + for (auto i : Range(points)) + if (numGroups(sel[i]) == 1) + identifications.Add(el[i], el[i + points.Size()], identnr); } - } - for(int i = 1; i<=mesh.GetNFD(); i++) - if(boundaries.Test(i)) - par_surfid.Append(i); - } - else - { - auto & surfids = *get_if>(¶ms.boundary); - for(auto id : surfids) - par_surfid.Append(id); - } - if(string* mat = get_if(¶ms.new_material); mat) - par_new_mat = { { ".*", *mat } }; - else - par_new_mat = *get_if>(¶ms.new_material); + } + Element2d newel = sel; + // check if the original points have a closesurface identification, if so, we need to also identify the corresponding new points - if(params.project_boundaries.has_value()) + for (auto pi : sel.PNums()) + for (auto pj : sel.PNums()) + { + if (pi == pj) + continue; + + auto nr = identifications.Get(pi, pj); + if (nr == 0) + continue; + + auto newpi = newPoint(pi); + auto newpj = newPoint(pj); + identifications.Add(newpi, newpj, nr); + } + + for (auto i : Range(points)) + newel[i] = newPoint(sel[i], -1, groups[i]); + newel.SetIndex(si_map[sel.GetIndex()]); + new_sels.Append(newel); + } + if (is_boundary_moved.Test(sel.GetIndex())) { - auto proj_bnd = *params.project_boundaries; - if(string* s = get_if(&proj_bnd); s) + auto& sel = mesh[si]; + auto pnums = sel.PNums(); + if (pnums.Size() == 4) { - regex pattern(*s); - for(int i = 1; i<=mesh.GetNFD(); i++) - if(regex_match(mesh.GetFaceDescriptor(i).GetBCName(), pattern)) - par_project_boundaries.Append(i); + // check if there are closesurface identifications, if so, we need to also identify the corresponding new points + for (auto i : Range(pnums)) + for (auto j : Range(i + 1, pnums.Size())) + { + if (i == j) + continue; + + auto pi = pnums[i]; + auto pj = pnums[j]; + + auto nr = identifications.Get(pi, pj); + if (nr == 0) + continue; + if (identifications.GetType(nr) != Identifications::CLOSESURFACES) + continue; + + auto newpi = hasMoved(pi) ? newPoint(pi) : pi; + auto newpj = hasMoved(pj) ? newPoint(pj) : pj; + identifications.Add(newpi, newpj, nr); + } + } + for (auto& p : pnums) + if (hasMoved(p)) + p = newPoint(p); + } + } + + for (SegmentIndex sei = 0; sei < nseg; sei++) + { + auto& seg = segments[sei]; + if (is_boundary_moved.Test(seg.si)) + { + // cout << "moved setg " << seg << endl; + for (auto& p : seg.PNums()) + if (hasMoved(p)) + p = newPoint(p); + } + } + + // fill holes in surface mesh at special boundary points (i.e. points with >=4 + // adjacent boundary faces) + auto p2sel = ngcore::CreateSortedTable( + new_sels.Range(), + [&] (auto& table, SurfaceElementIndex ei) { + for (PointIndex pi : new_sels[ei].PNums()) + table.Add(pi, ei); + }, + mesh.GetNP()); + + for (auto& [special_pi, special_point] : special_boundary_points) + { + if (special_point.growth_groups.Size() != 2) + throw Exception("special_point.growth_groups.Size() != 2"); + + // Special points are split into two new points, when mapping a surface + // element, we choose the closer one to the center. Now, find points which + // are mapped to both new points (for different surface elements they belong + // to). At exactly these points we need to insert new surface elements to + // fill the hole. + std::map, 2>> close_group; + for (auto sei : p2sel[special_pi]) + { + const auto& sel = mesh[sei]; + for (auto p : sel.PNums()) + if (p != special_pi) + close_group[sel.GetIndex()][getClosestGroup(special_pi, sei)].insert( + p); + } + + for (auto [fi, groups] : close_group) + { + const auto mapped_fi = si_map[fi]; + std::set common_points; + for (auto pi : groups[0]) + if (groups[1].count(pi) == 1) + common_points.insert(pi); + if (common_points.size() > 0) + { + auto pi_common = mapto[*common_points.begin()].Last(); + auto new_special_pi0 = special_point.growth_groups[0].new_points.Last(); + auto new_special_pi1 = special_point.growth_groups[1].new_points.Last(); + for (auto sei : p2sel[pi_common]) + { + if (mesh[sei].GetIndex() == mapped_fi && mesh[sei].PNums().Contains(new_special_pi0)) + { + auto sel = mesh[sei]; + sel.Invert(); + for (auto& pi : sel.PNums()) + if (pi != pi_common && pi != new_special_pi0) + pi = new_special_pi1; + new_sels.Append(sel); + } + } + } + } + } + + for (auto& [pi, special_point] : special_boundary_points) + { + if (special_point.growth_groups.Size() != 2) + throw Exception("special_point.growth_groups.Size() != 2"); + for (auto igroup : Range(2)) + { + auto& group = special_point.growth_groups[igroup]; + std::set faces; + for (auto face : group.faces) + faces.insert(si_map[face]); + auto pi_new = group.new_points.Last(); + auto pi_new_other = + special_point.growth_groups[1 - igroup].new_points.Last(); + for (auto sei : p2sel[pi_new]) + faces.erase(mesh[sei].GetIndex()); + for (auto face : faces) + for (auto seg : new_segments) + { + if ( // seg.si == face + (seg[0] == pi_new || seg[1] == pi_new) && (seg[0] != pi_new_other && seg[1] != pi_new_other)) + { + bool is_correct_face = false; + auto pi_other = seg[0] == pi_new ? seg[1] : seg[0]; + for (auto sei : p2sel[pi_other]) + { + if (mesh[sei].GetIndex() == face) + { + is_correct_face = true; + break; + } + } + if (is_correct_face) + { + Element2d sel; + sel[0] = seg[1]; + sel[1] = seg[0]; + sel[2] = pi_new_other; + sel.SetIndex(face); + new_sels.Append(sel); + } + } + } + } + } +} + +void BoundaryLayerTool ::SetDomInOut() +{ + if (insert_only_volume_elements) + return; + for (auto i : Range(1, nfd_old + 1)) + if (moved_surfaces.Test(i)) + { + if (auto dom = mesh.GetFaceDescriptor(si_map[i]).DomainIn(); + dom > ndom_old) + mesh.GetFaceDescriptor(i).SetDomainOut(dom); + else + mesh.GetFaceDescriptor(i).SetDomainIn( + mesh.GetFaceDescriptor(si_map[i]).DomainOut()); + } +} + +void BoundaryLayerTool ::SetDomInOutSides() +{ + // Set the domin/domout entries for face descriptors on the "side" of new boundary layers + if (insert_only_volume_elements) + return; + BitArray done(mesh.GetNFD() + 1); + done.Clear(); + + std::map inv_si_map; + + for (auto i : Range(si_map.Size())) + inv_si_map[si_map[i]] = i; + + for (auto sei : Range(mesh.SurfaceElements())) + { + auto& sel = mesh[sei]; + auto index = sel.GetIndex(); + if (done.Test(index)) + continue; + done.SetBit(index); + auto& fd = mesh.GetFaceDescriptor(index); + if (fd.DomainIn() != -1) + continue; + + // First check if there are adjacent volume elements, if so, use their domains + int e1 = 0, e2 = 0; + mesh.GetTopology().GetSurface2VolumeElement(sei + 1, e1, e2); + + int dom[2] = {-1, -1}; + + if (e1) + dom[0] = mesh.VolumeElement(e1).GetIndex(); + if (e2) + dom[1] = mesh.VolumeElement(e2).GetIndex(); + + auto& fd_old = mesh.GetFaceDescriptor(inv_si_map[index]); + int dom_old[2] = {fd_old.DomainIn(), fd_old.DomainOut()}; + + for (auto i : Range(2)) + { + if (dom[i] != -1) + continue; // adjacent volume element -> done + if (dom_old[i] == 0) + { + // outer boundary -> keep 0 + dom[i] = 0; + continue; + } + + // Check if the old domain adjacent to this face gets a new boundary layer domain, if so, use that number + int dom_new = domains.Test(dom_old[i]) ? new_mat_nrs[dom_old[i]] : dom_old[i]; + + // This case is tested by test_boundarylayer.py::test_pyramids[False] -> look at the generated mesh to understand the text below :) + // Special case check here: when growing "outside" the new face could have the same domain on both sides (before adding blayer elements). + // Thus we don't know in advance on which side the mapped domain will be. So check, if the other domain has already prisms (adjacent vol elements) with mapped domain. If so, use the original domain instead. + if (dom[1 - i] != dom_new) + { + dom[i] = dom_new; } else - { - for(auto id : *get_if>(&proj_bnd)) - par_project_boundaries.Append(id); - } + { + dom[i] = dom_old[i]; + } } - if(double* height = get_if(¶ms.thickness); height) + fd.SetDomainIn(dom[0]); + fd.SetDomainOut(dom[1]); + } +} + +void BoundaryLayerTool ::AddSegments() +{ + auto& new_segs = + insert_only_volume_elements ? new_segments_on_moved_bnd : new_segments; + + if (params.disable_curving) + { + for (auto& seg : old_segments) + if (mapto[seg[0]].Size() || mapto[seg[1]].Size()) + { + seg.epgeominfo[0].edgenr = -1; + seg.epgeominfo[1].edgenr = -1; + } + + for (auto& seg : segments) + if (is_edge_moved[seg.edgenr]) + { + seg.epgeominfo[0].edgenr = -1; + seg.epgeominfo[1].edgenr = -1; + } + + for (auto& seg : new_segs) { - par_heights.Append(*height); + seg.epgeominfo[0].edgenr = -1; + seg.epgeominfo[1].edgenr = -1; + } + } + + if (have_single_segments) + MergeAndAddSegments(mesh, segments, new_segs); + else + { + mesh.LineSegments() = segments; + for (auto& seg : new_segs) + mesh.AddSegment(seg); + } +} + +void BoundaryLayerTool ::AddSurfaceElements() +{ + for (auto& sel : + insert_only_volume_elements ? new_sels_on_moved_bnd : new_sels) + mesh.AddSurfaceElement(sel); +} + +void BoundaryLayerTool ::ProcessParameters() +{ + if (int* bc = get_if(¶ms.boundary); bc) + { + for (int i = 1; i <= mesh.GetNFD(); i++) + if (mesh.GetFaceDescriptor(i).BCProperty() == *bc) + par_surfid.Append(i); + } + else if (string* s = get_if(¶ms.boundary); s) + { + regex pattern(*s); + BitArray boundaries(mesh.GetNFD() + 1); + boundaries.Clear(); + for (int i = 1; i <= mesh.GetNFD(); i++) + { + auto& fd = mesh.GetFaceDescriptor(i); + if (regex_match(fd.GetBCName(), pattern)) + { + boundaries.SetBit(i); + auto dom_pattern = get_if(¶ms.domain); + // only add if adjacent to domain + if (dom_pattern) + { + regex pattern(*dom_pattern); + bool mat1_match = + fd.DomainIn() > 0 && regex_match(mesh.GetMaterial(fd.DomainIn()), pattern); + bool mat2_match = + fd.DomainOut() > 0 && regex_match(mesh.GetMaterial(fd.DomainOut()), pattern); + // if boundary is inner or outer remove from list + if (mat1_match == mat2_match) + boundaries.Clear(i); + // if((fd.DomainIn() > 0 && + // regex_match(mesh.GetMaterial(fd.DomainIn()), pattern)) || + // (fd.DomainOut() > 0 && + // regex_match(self.GetMaterial(fd.DomainOut()), pattern))) + // boundaries.Clear(i); + // par_surfid.Append(i); + } + // else + // par_surfid.Append(i); + } + } + for (int i = 1; i <= mesh.GetNFD(); i++) + if (boundaries.Test(i)) + par_surfid.Append(i); + } + else + { + auto& surfids = *get_if>(¶ms.boundary); + for (auto id : surfids) + par_surfid.Append(id); + } + + insert_only_volume_elements = !params.new_material.has_value(); + if (params.new_material) + { + if (string* mat = get_if(&*params.new_material); mat) + par_new_mat = {{".*", *mat}}; + else + par_new_mat = *get_if>(&*params.new_material); + } + + if (params.project_boundaries.has_value()) + { + auto proj_bnd = *params.project_boundaries; + if (string* s = get_if(&proj_bnd); s) + { + regex pattern(*s); + for (int i = 1; i <= mesh.GetNFD(); i++) + if (regex_match(mesh.GetFaceDescriptor(i).GetBCName(), pattern)) + par_project_boundaries.Append(i); } else { - auto & heights = *get_if>(¶ms.thickness); - for(auto val : heights) - par_heights.Append(val); + for (auto id : *get_if>(&proj_bnd)) + par_project_boundaries.Append(id); } + } - int nr_domains = mesh.GetNDomains(); - domains.SetSize(nr_domains + 1); // one based - domains.Clear(); - if(string* pdomain = get_if(¶ms.domain); pdomain) + if (double* height = get_if(¶ms.thickness); height) + { + par_heights.Append(*height); + } + else + { + auto& heights = *get_if>(¶ms.thickness); + for (auto val : heights) + par_heights.Append(val); + } + + int nr_domains = mesh.GetNDomains(); + domains.SetSize(nr_domains + 1); // one based + domains.Clear(); + if (string* pdomain = get_if(¶ms.domain); pdomain) + { + regex pattern(*pdomain); + for (auto i : Range(1, nr_domains + 1)) + if (regex_match(mesh.GetMaterial(i), pattern)) + domains.SetBit(i); + } + else if (int* idomain = get_if(¶ms.domain); idomain) + { + domains.SetBit(*idomain); + } + else + { + for (auto i : *get_if>(¶ms.domain)) + domains.SetBit(i); + } + if (domains.NumSet() == 0) + return; + total_height = 0.0; + for (auto h : par_heights) + total_height += h; + + max_edge_nr = -1; + for (const auto& seg : mesh.LineSegments()) + if (seg.edgenr > max_edge_nr) + max_edge_nr = seg.edgenr; + + int ndom = mesh.GetNDomains(); + ndom_old = ndom; + + new_mat_nrs.SetSize(mesh.FaceDescriptors().Size() + 1); + new_mat_nrs = -1; + if (insert_only_volume_elements) + { + for (auto i : Range(1, mesh.GetNFD() + 1)) { - regex pattern(*pdomain); - for(auto i : Range(1, nr_domains+1)) - if(regex_match(mesh.GetMaterial(i), pattern)) - domains.SetBit(i); + auto& fd = mesh.GetFaceDescriptor(i); + auto domin = fd.DomainIn(); + auto domout = fd.DomainOut(); + for (int dom : {domin, domout}) + if (domains.Test(dom)) + { + if (params.outside) + { + dom = domin + domout - dom; + if (dom == 0) + throw NG_EXCEPTION("No new material specified for boundarylayer " + "on the outside of domain"); + } + new_mat_nrs[i] = dom; + } } - else if(int *idomain = get_if(¶ms.domain); idomain) + } + else + { + for (auto [bcname, matname] : par_new_mat) { - domains.SetBit(*idomain); + mesh.SetMaterial(++ndom, matname); + regex pattern(bcname); + for (auto i : Range(1, mesh.GetNFD() + 1)) + { + auto& fd = mesh.GetFaceDescriptor(i); + if (regex_match(fd.GetBCName(), pattern)) + new_mat_nrs[i] = ndom; + } } - else - { - for (auto i : *get_if>(¶ms.domain)) - domains.SetBit(i); - } - } + } - void BoundaryLayerTool :: Perform() - { - CreateNewFaceDescriptors(); - CalculateGrowthVectors(); - CreateFaceDescriptorsSides(); - auto segmap = BuildSegMap(); + if (!params.outside) + domains.Invert(); +} - auto in_surface_direction = ProjectGrowthVectorsOnSurface(); +void BoundaryLayerTool ::Perform() +{ + if (domains.NumSet() == 0) + return; + CreateNewFaceDescriptors(); + CalculateGrowthVectors(); + CreateFaceDescriptorsSides(); + auto segmap = BuildSegMap(); - if(params.limit_growth_vectors) - LimitGrowthVectorLengths(); + auto in_surface_direction = ProjectGrowthVectorsOnSurface(); - InterpolateGrowthVectors(); - FixVolumeElements(); - InsertNewElements(segmap, in_surface_direction); - SetDomInOut(); - AddSegments(); - mesh.GetTopology().ClearEdges(); - mesh.SetNextMajorTimeStamp(); - mesh.UpdateTopology(); - SetDomInOutSides(); - MeshingParameters mp; - mp.optimize3d ="m"; - mp.optsteps3d = 4; - OptimizeVolume(mp, mesh); - } + InsertNewElements(segmap, in_surface_direction); - void GenerateBoundaryLayer(Mesh& mesh, const BoundaryLayerParameters& blp) - { - static Timer timer("Create Boundarylayers"); - RegionTimer regt(timer); + SetDomInOut(); + AddSegments(); - BoundaryLayerTool tool(mesh, blp); - tool.Perform(); - } + mesh.CalcSurfacesOfNode(); + topo.SetBuildVertex2Element(true); + mesh.UpdateTopology(); + + InterpolateGrowthVectors(); + InterpolateSurfaceGrowthVectors(); + + AddSurfaceElements(); + + if (params.limit_growth_vectors) + LimitGrowthVectorLengths(); + + FixSurfaceElements(); + + for (auto [pi, data] : growth_vector_map) + { + auto [gw, height] = data; + mesh[pi] += height * (*gw); + } + + if (insert_only_volume_elements) + { + mesh.LineSegments() = old_segments; + } + + mesh.CalcSurfacesOfNode(); + mesh.GetTopology().ClearEdges(); + mesh.SetNextMajorTimeStamp(); + mesh.UpdateTopology(); + SetDomInOutSides(); +} + +void GenerateBoundaryLayer (Mesh& mesh, const BoundaryLayerParameters& blp) +{ + static Timer timer("Create Boundarylayers"); + RegionTimer regt(timer); + + BoundaryLayerTool tool(mesh, blp); + tool.Perform(); +} } // namespace netgen diff --git a/libsrc/meshing/boundarylayer.hpp b/libsrc/meshing/boundarylayer.hpp index 264819bd..3646e4fb 100644 --- a/libsrc/meshing/boundarylayer.hpp +++ b/libsrc/meshing/boundarylayer.hpp @@ -1,85 +1,110 @@ #ifndef NETGEN_BOUNDARYLAYER_HPP #define NETGEN_BOUNDARYLAYER_HPP +#include +#include +#include + namespace netgen { /// -DLL_HEADER extern void InsertVirtualBoundaryLayer (Mesh & mesh); +DLL_HEADER extern void InsertVirtualBoundaryLayer (Mesh& mesh); -/// Create a typical prismatic boundary layer on the given +/// Create a typical prismatic boundary layer on the given /// surfaces -DLL_HEADER void GenerateBoundaryLayer (Mesh & mesh, - const BoundaryLayerParameters & blp); +struct SpecialBoundaryPoint +{ + struct GrowthGroup + { + Array faces; + Vec<3> growth_vector; + Array new_points; -DLL_HEADER int /* new_domain_number */ GenerateBoundaryLayer2 (Mesh & mesh, int domain, const Array & thicknesses, bool should_make_new_domain=true, const Array & boundaries=Array{}); + GrowthGroup(FlatArray faces_, FlatArray> normals); + GrowthGroup(const GrowthGroup&) = default; + GrowthGroup() = default; + }; + Array growth_groups; + Vec<3> separating_direction; + + SpecialBoundaryPoint(const std::map>& normals); + SpecialBoundaryPoint() = default; +}; + +DLL_HEADER void GenerateBoundaryLayer (Mesh& mesh, + const BoundaryLayerParameters& blp); + +DLL_HEADER int /* new_domain_number */ GenerateBoundaryLayer2 (Mesh& mesh, int domain, const Array& thicknesses, bool should_make_new_domain = true, const Array& boundaries = Array{}); class BoundaryLayerTool { - public: - BoundaryLayerTool(Mesh & mesh_, const BoundaryLayerParameters & params_); - void ProcessParameters(); - void Perform(); +public: + BoundaryLayerTool(Mesh& mesh_, const BoundaryLayerParameters& params_); + void ProcessParameters (); + void Perform (); - protected: - Mesh & mesh; - MeshTopology & topo; - BoundaryLayerParameters params; - Array, PointIndex> growthvectors; - Table p2sel; + Mesh& mesh; + MeshTopology& topo; + BoundaryLayerParameters params; + Array, PointIndex> growthvectors; + std::map> non_bl_growth_vectors; + Table p2sel; - BitArray domains, is_edge_moved, is_boundary_projected, is_boundary_moved; - Array moved_segs; - int max_edge_nr, nfd_old, ndom_old; - Array new_mat_nrs; - BitArray moved_surfaces; - int np, nseg, nse, ne; - double height; + BitArray domains, is_edge_moved, is_boundary_projected, is_boundary_moved; + Array moved_segs; + int max_edge_nr, nfd_old, ndom_old; + Array new_mat_nrs; + BitArray moved_surfaces; + int np, nseg, nse, ne; + double total_height; + Array point_types; - // These parameters are derived from given BoundaryLayerParameters and the Mesh - Array par_heights; - Array par_surfid; - map par_new_mat; - Array par_project_boundaries; + // These parameters are derived from given BoundaryLayerParameters and the Mesh + Array par_heights; + Array par_surfid; + bool insert_only_volume_elements; + map par_new_mat; + Array par_project_boundaries; - bool have_single_segments; - Array segments, new_segments; + bool have_single_segments; + Array old_segments, segments, new_segments, new_segments_on_moved_bnd; + Array new_sels, new_sels_on_moved_bnd; + Array, PointIndex> mapto; + Array mapfrom; - Array surfacefacs; - Array si_map; - Array limits; + Array surfacefacs; + Array si_map; - // major steps called in Perform() - void CreateNewFaceDescriptors(); - void CreateFaceDescriptorsSides(); - void CalculateGrowthVectors(); - Array>, SegmentIndex> BuildSegMap(); + std::map special_boundary_points; + std::map*, double>> growth_vector_map; - BitArray ProjectGrowthVectorsOnSurface(); - void InterpolateSurfaceGrowthVectors(); - void InterpolateGrowthVectors(); - void LimitGrowthVectorLengths(); + // major steps called in Perform() + void CreateNewFaceDescriptors (); + void CreateFaceDescriptorsSides (); + void CalculateGrowthVectors (); + Array>, SegmentIndex> BuildSegMap (); - void InsertNewElements(FlatArray>, SegmentIndex> segmap, const BitArray & in_surface_direction); - void SetDomInOut(); - void SetDomInOutSides(); - void AddSegments(); - void FixVolumeElements(); + BitArray ProjectGrowthVectorsOnSurface (); + void InterpolateSurfaceGrowthVectors (); + void InterpolateGrowthVectors (); + void LimitGrowthVectorLengths (); + void FixSurfaceElements (); - // utility functions - array, 2> GetMappedSeg( PointIndex pi ); - ArrayMem, 4> GetFace( SurfaceElementIndex sei ); - ArrayMem, 4> GetMappedFace( SurfaceElementIndex sei ); - ArrayMem, 4> GetMappedFace( SurfaceElementIndex sei, int face ); + void InsertNewElements (FlatArray>, SegmentIndex> segmap, const BitArray& in_surface_direction); + void SetDomInOut (); + void SetDomInOutSides (); + void AddSegments (); + void AddSurfaceElements (); - Vec<3> getNormal(const Element2d & el) - { - auto v0 = mesh[el[0]]; - return Cross(mesh[el[1]]-v0, mesh[el[2]]-v0).Normalize(); - } + Vec<3> getNormal (const Element2d& el) + { + auto v0 = mesh[el[0]]; + return Cross(mesh[el[1]] - v0, mesh[el[2]] - v0).Normalize(); + } - Vec<3> getEdgeTangent(PointIndex pi, int edgenr); + Vec<3> getEdgeTangent (PointIndex pi, int edgenr, FlatArray segs); }; } // namespace netgen diff --git a/libsrc/meshing/boundarylayer2d.cpp b/libsrc/meshing/boundarylayer2d.cpp index 93f0cf36..77f0083d 100644 --- a/libsrc/meshing/boundarylayer2d.cpp +++ b/libsrc/meshing/boundarylayer2d.cpp @@ -1,5 +1,5 @@ #include -#include "meshing.hpp" +#include "boundarylayer.hpp" #include "meshing2.hpp" #include "../geom2d/csg2d.hpp" diff --git a/libsrc/meshing/boundarylayer_interpolate.cpp b/libsrc/meshing/boundarylayer_interpolate.cpp new file mode 100644 index 00000000..95d7897b --- /dev/null +++ b/libsrc/meshing/boundarylayer_interpolate.cpp @@ -0,0 +1,472 @@ +#include "boundarylayer.hpp" + +namespace netgen +{ + +namespace detail +{ +struct Neighbor +{ + PointIndex pi; + SurfaceElementIndex sei; + double weight; +}; +} // namespace detail + +Array> +BuildNeighbors (FlatArray points, const Mesh& mesh) +{ + auto p2sel = mesh.CreatePoint2SurfaceElementTable(); + + Array> neighbors(points.Size()); + + ArrayMem angles; + ArrayMem inv_dists; + for (auto i : points.Range()) + { + auto& p_neighbors = neighbors[i]; + auto pi = points[i]; + angles.SetSize(0); + inv_dists.SetSize(0); + for (auto sei : p2sel[pi]) + { + const auto& sel = mesh[sei]; + for (auto pi1 : sel.PNums()) + { + if (pi1 == pi) + continue; + auto pi2 = pi1; + for (auto pi_ : sel.PNums()) + { + if (pi_ != pi && pi_ != pi1) + { + pi2 = pi_; + break; + } + } + p_neighbors.Append({pi1, sei, 0.0}); + inv_dists.Append(1.0 / (mesh[pi1] - mesh[pi]).Length()); + auto dot = (mesh[pi1] - mesh[pi]).Normalize() * (mesh[pi2] - mesh[pi]).Normalize(); + angles.Append(acos(dot)); + } + } + double sum_inv_dist = 0.0; + for (auto inv_dist : inv_dists) + sum_inv_dist += inv_dist; + double sum_angle = 0.0; + for (auto angle : angles) + sum_angle += angle; + + double sum_weight = 0.0; + for (auto i : Range(inv_dists)) + { + p_neighbors[i].weight = + inv_dists[i] * angles[i] / sum_inv_dist / sum_angle; + sum_weight += p_neighbors[i].weight; + } + for (auto i : Range(inv_dists)) + p_neighbors[i].weight /= sum_weight; + } + return neighbors; +} + +void BoundaryLayerTool ::InterpolateGrowthVectors() +{ + point_types.SetSize(mesh.GetNP()); + for (auto p : mesh.Points().Range()) + point_types[p] = mesh[p].Type(); + + int new_max_edge_nr = max_edge_nr; + for (const auto& seg : segments) + if (seg.edgenr > new_max_edge_nr) + new_max_edge_nr = seg.edgenr; + for (const auto& seg : new_segments) + if (seg.edgenr > new_max_edge_nr) + new_max_edge_nr = seg.edgenr; + + auto getGW = [&] (PointIndex pi) -> Vec<3> { + if (growth_vector_map.count(pi) == 0) + growth_vector_map[pi] = {&growthvectors[pi], total_height}; + auto [gw, height] = growth_vector_map[pi]; + return height * (*gw); + }; + auto addGW = [&] (PointIndex pi, Vec<3> vec) { + if (growth_vector_map.count(pi) == 0) + growth_vector_map[pi] = {&growthvectors[pi], total_height}; + auto [gw, height] = growth_vector_map[pi]; + *gw += 1.0 / height * vec; + }; + + // interpolate tangential component of growth vector along edge + if (max_edge_nr >= new_max_edge_nr) + return; + + auto edgenr2seg = ngcore::CreateSortedTable( + Range(segments.Size() + new_segments.Size()), + [&] (auto& table, size_t segi) { + auto& seg = segi < segments.Size() + ? segments[segi] + : new_segments[segi - segments.Size()]; + table.Add(seg.edgenr, &seg); + }, + new_max_edge_nr + 1); + auto point2seg = ngcore::CreateSortedTable( + Range(segments.Size() + new_segments.Size()), + [&] (auto& table, size_t segi) { + auto& seg = segi < segments.Size() + ? segments[segi] + : new_segments[segi - segments.Size()]; + table.Add(seg[0], &seg); + table.Add(seg[1], &seg); + }, + mesh.GetNP()); + + for (auto edgenr : Range(1, new_max_edge_nr + 1)) + { + // "inner" edges between two flat faces are not treated as edges for interpolation + bool no_angles = true; + ArrayMem faces; + + for (auto* p_seg : edgenr2seg[edgenr]) + { + auto& seg = *p_seg; + faces.SetSize(0); + if (seg[0] <= p2sel.Size()) + { + for (auto sei : p2sel[seg[0]]) + if (moved_surfaces.Test(mesh[sei].GetIndex()) && p2sel[seg[1]].Contains(sei)) + faces.Append(sei); + } + + if (faces.Size() == 2) + { + auto n0 = getNormal(mesh[faces[0]]); + auto n1 = getNormal(mesh[faces[1]]); + if (n0 * n1 < 0.999) + no_angles = false; + } + else + { + no_angles = false; + } + } + if (no_angles) + { + for (auto* p_seg : edgenr2seg[edgenr]) + for (auto pi : p_seg->PNums()) + if (pi <= np && point_types[pi] == EDGEPOINT) + point_types[pi] = SURFACEPOINT; + continue; + } + } + + for (auto edgenr : Range(max_edge_nr + 1, new_max_edge_nr + 1)) + { + double edge_len = 0.; + bool any_grows = false; + + auto is_end_point = [&] (PointIndex pi) { + auto segs = point2seg[pi]; + if (segs.Size() == 1) + return true; + auto first_edgenr = (*segs[0]).edgenr; + for (auto* p_seg : segs) + if (p_seg->edgenr != first_edgenr) + return true; + return false; + }; + + Array points; + for (auto* p_seg : edgenr2seg[edgenr]) + { + auto& seg = *p_seg; + + if (getGW(seg[0]).Length2() != 0 || getGW(seg[1]).Length2() != 0) + any_grows = true; + + if (points.Size() == 0) + for (auto i : Range(2)) + if (is_end_point(seg[i])) + { + points.Append(seg[i]); + points.Append(seg[1 - i]); + edge_len += (mesh[seg[1]] - mesh[seg[0]]).Length(); + break; + } + } + + if (!any_grows) + { + PrintMessage(1, "BLayer: skip interpolating growth vectors at edge ", edgenr + 1); + continue; + } + + if (!points.Size()) + { + if (debugparam.debugoutput) + cerr << "Could not find startpoint for edge " << edgenr << endl; + continue; + } + + std::set points_set; + points_set.insert(points[0]); + points_set.insert(points[1]); + + bool point_found = true; + while (point_found) + { + if (is_end_point(points.Last())) + break; + point_found = false; + for (auto* p_seg : point2seg[points.Last()]) + { + const auto& seg = *p_seg; + if (seg.edgenr != edgenr) + continue; + auto plast = points.Last(); + if (plast != seg[0] && plast != seg[1]) + continue; + auto pnew = plast == seg[0] ? seg[1] : seg[0]; + if (pnew == points[0] && points.Size() > 1) + { + } + if (points_set.count(pnew) > 0 && (pnew != points[0] || points.Size() == 2)) + continue; + edge_len += (mesh[points.Last()] - mesh[pnew]).Length(); + points.Append(pnew); + points_set.insert(pnew); + point_found = true; + break; + } + } + if (!point_found) + { + if (debugparam.debugoutput) + { + cerr << "Could not find connected list of line segments for edge " + << edgenr << endl; + cerr << "current points: " << endl + << points << endl; + } + continue; + } + + if (getGW(points[0]).Length2() == 0 && getGW(points.Last()).Length2() == 0) + continue; + + // tangential part of growth vectors + auto t1 = (mesh[points[1]] - mesh[points[0]]).Normalize(); + auto gt1 = getGW(points[0]) * t1 * t1; + auto t2 = + (mesh[points.Last()] - mesh[points[points.Size() - 2]]).Normalize(); + auto gt2 = getGW(points.Last()) * t2 * t2; + + double len = 0.; + for (auto i : IntRange(1, points.Size() - 1)) + { + auto pi = points[i]; + len += (mesh[pi] - mesh[points[i - 1]]).Length(); + auto t = getEdgeTangent(pi, edgenr, point2seg[pi]); + auto lam = len / edge_len; + auto interpol = (1 - lam) * (gt1 * t) * t + lam * (gt2 * t) * t; + addGW(pi, interpol); + } + } +} + +void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() +{ + static Timer tall("InterpolateSurfaceGrowthVectors"); + RegionTimer rtall(tall); + static Timer tsmooth("InterpolateSurfaceGrowthVectors-Smoothing"); + auto np_old = this->np; + auto np = mesh.GetNP(); + + auto hasMoved = [&] (PointIndex pi) { + return (pi - PointIndex::BASE >= np_old) || mapto[pi].Size() > 0 || special_boundary_points.count(pi); + }; + + std::set points_set; + for (const auto& sel : mesh.SurfaceElements()) + { + for (auto pi : sel.PNums()) + if (point_types[pi] == SURFACEPOINT && hasMoved(pi)) + points_set.insert(pi); + } + + Array points; + for (auto pi : points_set) + points.Append(pi); + QuickSort(points); + + // smooth tangential part of growth vectors from edges to surface elements + Array, PointIndex> corrections(mesh.GetNP()); + corrections = 0.0; + RegionTimer rtsmooth(tsmooth); + auto neighbors = BuildNeighbors(points, mesh); + + Array, SurfaceElementIndex> surf_normals(mesh.GetNSE()); + for (auto sei : mesh.SurfaceElements().Range()) + surf_normals[sei] = getNormal(mesh[sei]); + + BitArray interpolate_tangent(mesh.GetNP() + 1); + interpolate_tangent = false; + for (auto pi : points) + { + for (auto sei : p2sel[pi]) + if (is_boundary_moved[mesh[sei].GetIndex()]) + interpolate_tangent.SetBit(pi); + } + + constexpr int N_STEPS = 64; + for ([[maybe_unused]] auto i : Range(N_STEPS)) + { + for (auto i : points.Range()) + { + auto pi = points[i]; + auto& p_neighbors = neighbors[i]; + + ArrayMem, 20> g_vectors; + double max_len = 0.0; + double sum_len = 0.0; + + // average only tangent component on new bl points, average whole growth + // vector otherwise + bool do_average_tangent = true; + for (const auto& s : p_neighbors) + { + auto gw_other = growthvectors[s.pi] + corrections[s.pi]; + if (do_average_tangent) + { + auto n = surf_normals[s.sei]; + gw_other = gw_other - (gw_other * n) * n; + } + auto v = gw_other; + auto len = v.Length2(); + sum_len += len; + max_len = max(max_len, len); + g_vectors.Append(v); + } + + if (max_len == 0.0) + continue; + + double lambda = 0; + if (i > N_STEPS / 4.) + lambda = 2.0 * (i - N_STEPS / 4.) / (N_STEPS / 2.); + lambda = min(1.0, lambda); + + auto& correction = corrections[pi]; + correction = 0.0; + for (const auto i : p_neighbors.Range()) + { + auto v = g_vectors[i]; + double weight = lambda * p_neighbors[i].weight + (1.0 - lambda) * v.Length2() / sum_len; + correction += weight * v; + } + + if (!do_average_tangent) + correction -= growthvectors[pi]; + } + } + + for (auto pi : points) + growthvectors[pi] += corrections[pi]; +} + +void BoundaryLayerTool ::FixSurfaceElements() +{ + static Timer tall("FixSurfaceElements"); + RegionTimer rtall(tall); + auto np_old = this->np; + auto np = mesh.GetNP(); + + non_bl_growth_vectors.clear(); + + auto getGW = [&] (PointIndex pi) -> Vec<3> { + // return growthvectors[pi]; + if (growth_vector_map.count(pi) == 0) + { + non_bl_growth_vectors[pi] = .0; + growth_vector_map[pi] = {&non_bl_growth_vectors[pi], 1.0}; + } + auto [gw, height] = growth_vector_map[pi]; + return height * (*gw); + }; + + auto addGW = [&] (PointIndex pi, Vec<3> vec) { + if (growth_vector_map.count(pi) == 0) + { + non_bl_growth_vectors[pi] = .0; + growth_vector_map[pi] = {&non_bl_growth_vectors[pi], 1.0}; + } + auto [gw, height] = growth_vector_map[pi]; + *gw += 1.0 / height * vec; + }; + + std::set points_set; + // only smooth over old surface elements + for (SurfaceElementIndex sei : Range(nse)) + { + const auto& sel = mesh[sei]; + if (sel.GetNP() == 3 && is_boundary_moved[sel.GetIndex()]) + for (auto pi : sel.PNums()) + if (point_types[pi] == SURFACEPOINT) + points_set.insert(pi); + } + + Array points; + for (auto pi : points_set) + points.Append(pi); + QuickSort(points); + + Array, PointIndex> corrections(mesh.GetNP()); + corrections = 0.0; + + auto neighbors = BuildNeighbors(points, mesh); + + constexpr int N_STEPS = 32; + for ([[maybe_unused]] auto i : Range(N_STEPS)) + { + for (auto i : points.Range()) + { + auto pi = points[i]; + auto& p_neighbors = neighbors[i]; + + ArrayMem, 20> g_vectors; + double max_len = 0.0; + double sum_len = 0.0; + + for (const auto& s : p_neighbors) + { + auto v = getGW(s.pi) + corrections[s.pi]; + auto len = v.Length2(); + sum_len += len; + max_len = max(max_len, len); + g_vectors.Append(v); + } + + if (max_len == 0.0) + continue; + + double lambda = 0; + if (i > N_STEPS / 4.) + lambda = 2.0 * (i - N_STEPS / 4.) / (N_STEPS / 2.); + lambda = min(1.0, lambda); + + auto& correction = corrections[pi]; + correction = 0.0; + for (const auto i : p_neighbors.Range()) + { + auto v = g_vectors[i]; + double weight = lambda * p_neighbors[i].weight + (1.0 - lambda) * v.Length2() / sum_len; + correction += weight * v; + } + } + } + + for (auto pi : points) + addGW(pi, corrections[pi]); +} + +} // namespace netgen diff --git a/libsrc/meshing/boundarylayer_limiter.hpp b/libsrc/meshing/boundarylayer_limiter.hpp new file mode 100644 index 00000000..dacbdc09 --- /dev/null +++ b/libsrc/meshing/boundarylayer_limiter.hpp @@ -0,0 +1,752 @@ +#include "boundarylayer.hpp" +#include + +namespace netgen +{ + +struct Intersection_ +{ + bool is_intersecting = false; + double lam0 = -1, lam1 = -1; + Point<3> p; + double bary[3]; + operator bool() const { return is_intersecting; } +}; + +struct GrowthVectorLimiter +{ + typedef std::array, 2> Seg; + typedef std::array, 3> Trig; + + BoundaryLayerTool& tool; + const BoundaryLayerParameters& params; + Mesh& mesh; + double height; + Array limits; + FlatArray, PointIndex> growthvectors; + BitArray changed_domains; + unique_ptr> tree; + Array map_from; + Table p2sel; + + GrowthVectorLimiter(BoundaryLayerTool& tool_) + : tool(tool_), params(tool_.params), mesh(tool_.mesh), height(tool_.total_height), growthvectors(tool_.growthvectors), map_from(mesh.Points().Size()) + { + changed_domains = tool.domains; + if (!params.outside) + changed_domains.Invert(); + + map_from = tool.mapfrom; + p2sel = ngcore::CreateSortedTable( + tool.new_sels.Range(), + [&] (auto& table, SurfaceElementIndex ei) { + for (PointIndex pi : tool.new_sels[ei].PNums()) + table.Add(pi, ei); + }, + mesh.GetNP()); + } + + auto SurfaceElementsRange () { return Range(tool.nse + tool.new_sels.Size()); } + + void WriteErrorMesh (string name) + { + if (!debugparam.write_mesh_on_error) + return; + Mesh out_mesh; + out_mesh = mesh; + for (auto [pi, data] : tool.growth_vector_map) + { + auto [gw, height] = data; + out_mesh[pi] += limits[pi] * height * (*gw); + } + out_mesh.Save(name); + } + + const auto& Get (SurfaceElementIndex sei) + { + if (sei < tool.nse) + return mesh[sei]; + return tool.new_sels[sei - tool.nse]; + } + + std::pair GetMinMaxLimit (SurfaceElementIndex sei) + { + const auto& sel = Get(sei); + double min_limit = GetLimit(sel[0]); + double max_limit = min_limit; + for (auto i : IntRange(1, sel.GetNP())) + { + auto limit = GetLimit(sel[i]); + min_limit = min(min_limit, limit); + max_limit = max(max_limit, limit); + } + return {min_limit, max_limit}; + } + + double GetLimit (PointIndex pi) + { + if (pi <= tool.np) + return limits[pi]; + return limits[map_from[pi]]; + } + + bool SetLimit (PointIndex pi, double new_limit) + { + double& limit = (pi <= tool.np) ? limits[pi] : limits[map_from[pi]]; + if (limit <= new_limit) + return false; + limit = new_limit; + return true; + } + + bool ScaleLimit (PointIndex pi, double factor) + { + double& limit = (pi <= tool.np) ? limits[pi] : limits[map_from[pi]]; + return SetLimit(pi, limit * factor); + } + + Vec<3> GetVector (PointIndex pi_to, double shift = 1., bool apply_limit = false) + { + auto [gw, height] = tool.growth_vector_map[pi_to]; + if (apply_limit) + shift *= GetLimit(pi_to); + return shift * height * (*gw); + } + + Point<3> GetPoint (PointIndex pi_to, double shift = 1., bool apply_limit = false) + { + if (pi_to <= tool.np || tool.growth_vector_map.count(pi_to) == 0) + return mesh[pi_to]; + + return mesh[pi_to] + GetVector(pi_to, shift, apply_limit); + } + + Point<3> GetMappedPoint (PointIndex pi_from, double shift = 1., bool apply_limit = false) + { + auto pi_to = tool.mapto[pi_from].Last(); + return GetPoint(pi_to, shift, apply_limit); + } + + Seg GetMappedSeg (PointIndex pi_from, double shift = 1.) + { + return {mesh[pi_from], GetMappedPoint(pi_from, shift)}; + } + + Seg GetSeg (PointIndex pi_to, double shift = 1., bool apply_limit = false) + { + return {GetPoint(pi_to, 0), GetPoint(pi_to, shift, apply_limit)}; + } + + Trig GetTrig (SurfaceElementIndex sei, double shift = 0.0, bool apply_limit = false) + { + auto sel = Get(sei); + Trig trig; + for (auto i : Range(3)) + trig[i] = GetPoint(sel[i], shift, apply_limit); + return trig; + } + + Trig GetMappedTrig (SurfaceElementIndex sei, double shift = 0.0) + { + auto sel = Get(sei); + Trig trig; + for (auto i : Range(3)) + trig[i] = GetMappedPoint(sel[i], shift); + return trig; + } + + Trig GetSideTrig (SurfaceElementIndex sei, int index, double shift = 0.0, bool grow_first_vertex = true) + { + auto trig = GetMappedTrig(sei, 0.0); + auto sel = Get(sei); + auto index1 = (index + 1) % 3; + if (!grow_first_vertex) + index1 = (index + 2) % 3; + trig[index] = GetMappedPoint(sel[index1], shift, true); + return trig; + } + + array GetSideTrigs (SurfaceElementIndex sei, int i0, double shift = 0.0) + { + auto trig = GetMappedTrig(sei, 0.0); + array trigs{trig, trig, trig, trig}; + + auto sel = Get(sei); + auto i1 = (i0 + 1) % 3; + auto i2 = (i0 + 2) % 3; + auto p1 = GetMappedPoint(sel[i1], shift, true); + auto p2 = GetMappedPoint(sel[i2], shift, true); + + // create four trigs to span the quad from i1,i2 and their shifted points + // i1, i2, shifted i1 + trigs[0][i0] = p1; + + // i1, i2, shifted i2 + trigs[1][i0] = p2; + + // i1, shifted i1, shifted i2 + trigs[2][i0] = p1; + trigs[2][i2] = p2; + + // i2, shifted i1, shifted i2 + trigs[2][i0] = p2; + trigs[2][i1] = p1; + + return trigs; + } + + static constexpr double INTERSECTION_SAFETY = .9; + bool LimitGrowthVector (PointIndex pi_to, SurfaceElementIndex sei, double trig_shift, double seg_shift, bool check_prism_sides = false) + { + auto pi_from = map_from[pi_to]; + if (!pi_from.IsValid()) + return false; + + for (auto pi : Get(sei).PNums()) + { + if (pi == pi_from) + return false; + if (map_from[pi] == pi_from) + return false; + } + + if (check_prism_sides || trig_shift > .0) + { + auto [trig_min_limit, trig_max_limit] = GetMinMaxLimit(sei); + if (GetLimit(pi_to) < trig_min_limit) + return false; + + auto getTrigs = [&] (double scaling = 1.0) -> ArrayMem { + ArrayMem trigs; + if (check_prism_sides) + for (auto i : Range(3)) + for (auto trig : GetSideTrigs(sei, i, scaling * trig_shift)) + trigs.Append(trig); + else + trigs.Append(GetTrig(sei, scaling * trig_shift, true)); + return trigs; + }; + + if (!check_prism_sides) + { + // If the growth vectors of all points are pointing in the same direction, + // an intersection means, we also have an intersection with a prism side face + // this is an extra check and handled later + auto seg = GetSeg(pi_to, 1.0, false); + auto gw = seg[1] - seg[0]; + + bool have_same_growth_direction = true; + for (auto pi : Get(sei).PNums()) + { + auto p_seg = GetSeg(pi, 1.0, false); + auto p_gw = p_seg[1] - p_seg[0]; + have_same_growth_direction &= (gw * p_gw) > 0; + } + if (have_same_growth_direction) + return false; + } + + double scaling = 1.0; + while (true) + { + bool have_intersection = false; + auto seg = GetSeg(pi_to, scaling * seg_shift, true); + for (auto trig : getTrigs(scaling)) + have_intersection |= isIntersectingTrig(seg, trig); + if (!have_intersection) + break; + scaling *= 0.9; + } + if (scaling == 1.0) + return false; + + double new_limit = scaling * max(GetLimit(pi_to), trig_max_limit); + SetLimit(pi_to, new_limit); + for (auto pi : Get(sei).PNums()) + SetLimit(pi, new_limit); + return true; + } + else + { + auto seg = GetSeg(pi_to, seg_shift, false); + auto trig = GetTrig(sei, 0.0); + auto intersection = isIntersectingTrig(seg, trig); + // checking with original surface elements -> allow only half the distance + auto new_seg_limit = 0.40 * intersection.lam0 * seg_shift; + if (intersection && new_seg_limit < GetLimit(pi_from)) + return SetLimit(pi_from, new_seg_limit); + return false; + } + } + + void EqualizeLimits (double factor = .5) + { + static Timer t("GrowthVectorLimiter::EqualizeLimits"); + PrintMessage(5, "GrowthVectorLimiter - equalize limits"); + RegionTimer reg(t); + if (factor == 0.0) + return; + for (PointIndex pi : IntRange(tool.np, mesh.GetNP())) + { + auto pi_from = map_from[pi]; + std::set pis; + for (auto sei : p2sel[pi]) + for (auto pi_ : tool.new_sels[sei].PNums()) + pis.insert(pi_); + ArrayMem limits; + for (auto pi1 : pis) + { + auto limit = GetLimit(pi1); + if (limit > 0.0) + limits.Append(GetLimit(pi1)); + } + + if (limits.Size() == 0) + continue; + + double average = 0.0; + for (auto l : limits) + average += l; + average /= limits.Size(); + + SetLimit(pi, factor * average + (1.0 - factor) * GetLimit(pi)); + } + } + + void LimitSelfIntersection (double safety = 1.4) + { + static Timer t("GrowthVectorLimiter::LimitSelfIntersection"); + PrintMessage(5, "GrowthVectorLimiter - self intersection"); + RegionTimer reg(t); + // check for self-intersection within new elements (prisms/hexes) + auto isIntersecting = [&] (SurfaceElementIndex sei, double shift) { + // checks if surface element is self intersecting when growing with factor + // shift + + // ignore new surface elements, side trigs are only built + // from original surface elements + if (sei >= tool.nse) + return false; + const auto sel = Get(sei); + auto np = sel.GetNP(); + for (auto i : Range(np)) + { + if (sel[i] > tool.np) + return false; + if (tool.mapto[sel[i]].Size() == 0) + return false; + } + for (auto i : Range(np)) + { + auto seg = GetMappedSeg(sel[i], shift * limits[sel[i]]); + for (auto fi : Range(np - 2)) + { + for (auto side : {true, false}) + { + auto trig = GetSideTrig(sei, i + fi, 1.0, side); + if (isIntersectingPlane(seg, trig)) + return true; + } + } + } + return false; + }; + + for (SurfaceElementIndex sei : mesh.SurfaceElements().Range()) + { + auto sel = mesh[sei]; + if (!tool.moved_surfaces[sel.GetIndex()]) + continue; + if (sel.GetNP() == 4) + continue; + + const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex()); + auto np = sel.GetNP(); + + double shift = 1.0; + const double step_factor = 0.9; + while (isIntersecting(sei, shift * safety)) + { + shift *= step_factor; + double max_limit = 0; + for (auto i : Range(np)) + max_limit = max(max_limit, limits[sel[i]]); + for (auto i : Range(np)) + if (max_limit == limits[sel[i]]) + ScaleLimit(sel[i], step_factor); + // if (max_limit < 0.01) break; + } + } + } + + // checks if a segment is intersecting a plane, spanned by three points, lam + // will be set s.t. p_intersect = seg[0] + lam * (seg[1]-seg[0]) + Intersection_ isIntersectingPlane (const Seg& seg, + const Trig& trig) + { + auto t1 = trig[1] - trig[0]; + auto t2 = trig[2] - trig[0]; + auto n = Cross(t1, t2); + auto v0n = (seg[0] - trig[0]) * n; + auto v1n = (seg[1] - trig[0]) * n; + + Intersection_ intersection; + intersection.lam0 = -v0n / (v1n - v0n); + intersection.p = seg[0] + intersection.lam0 * (seg[1] - seg[0]); + intersection.is_intersecting = (v0n * v1n < 0) && (intersection.lam0 > -1e-8) && (intersection.lam0 < 1 + 1e-8); + + return intersection; + } + + Intersection_ isIntersectingTrig (const Seg& seg, const Trig& trig) + { + auto intersection = isIntersectingPlane(seg, trig); + if (!intersection) + return intersection; + + auto p = seg[0] + intersection.lam0 * (seg[1] - seg[0]) - trig[0]; + + Vec3d col1 = trig[1] - trig[0]; + Vec3d col2 = trig[2] - trig[0]; + Vec3d col3 = Cross(col1, col2); + Vec3d rhs = p; + Vec3d bary; + SolveLinearSystem(col1, col2, col3, rhs, bary); + + intersection.lam1 = 0; + double eps = 1e-4; + if (bary.X() >= -eps && bary.Y() >= -eps && bary.X() + bary.Y() <= 1 + eps) + { + intersection.bary[0] = bary.X(); + intersection.bary[1] = bary.Y(); + intersection.bary[2] = 1.0 - bary.X() - bary.Y(); + } + else + intersection.is_intersecting = false; + return intersection; + } + + Intersection_ isIntersectingTrig (PointIndex pi_from, PointIndex pi_to, SurfaceElementIndex sei, double shift = 0.0) + { + return isIntersectingTrig(GetSeg(pi_from, pi_to), GetTrig(sei, shift)); + } + + void BuildSearchTree (double trig_shift) + { + static Timer t("BuildSearchTree"); + RegionTimer rt(t); + Box<3> bbox(Box<3>::EMPTY_BOX); + for (PointIndex pi : mesh.Points().Range()) + { + bbox.Add(mesh[pi]); + bbox.Add(GetPoint(pi, 1.1)); + } + + tree = make_unique>(bbox); + + for (auto sei : SurfaceElementsRange()) + { + const auto& sel = Get(sei); + auto sel_index = sel.GetIndex(); + + Box<3> box(Box<3>::EMPTY_BOX); + for (auto pi : sel.PNums()) + { + box.Add(GetPoint(pi, 0.)); + box.Add(GetPoint(pi, trig_shift * GetLimit(pi))); + } + tree->Insert(box, sei); + } + } + + template + void FindTreeIntersections (double trig_shift, double seg_shift, TFunc f, BitArray* relevant_points = nullptr) + { + static Timer t("GrowthVectorLimiter::FindTreeIntersections"); + RegionTimer rt(t); + BuildSearchTree(trig_shift); + auto np_new = mesh.Points().Size(); + int counter = 0; + for (auto i : IntRange(tool.np, np_new)) + { + PointIndex pi_to = i + PointIndex::BASE; + PointIndex pi_from = map_from[pi_to]; + if (!pi_from.IsValid()) + throw Exception("Point not mapped"); + + if (relevant_points && !relevant_points->Test(pi_to) && !relevant_points->Test(pi_from)) + continue; + + Box<3> box(Box<3>::EMPTY_BOX); + auto seg = GetSeg(pi_to, seg_shift); + + box.Add(GetPoint(pi_to, 0)); + box.Add(GetPoint(pi_to, GetLimit(pi_from))); + tree->GetFirstIntersecting(box.PMin(), box.PMax(), [&] (SurfaceElementIndex sei) { + const auto& sel = Get(sei); + if (sel.PNums().Contains(pi_from)) + return false; + if (sel.PNums().Contains(pi_to)) + return false; + counter++; + f(pi_to, sei); + return false; + }); + } + } + + void FixIntersectingSurfaceTrigs () + { + static Timer t("GrowthVectorLimiter::FixIntersectingSurfaceTrigs"); + RegionTimer reg(t); + // check if surface trigs are intersecting each other + bool changed = true; + while (changed) + { + changed = false; + Point3d pmin, pmax; + mesh.GetBox(pmin, pmax); + BoxTree<3, SurfaceElementIndex> setree(pmin, pmax); + + for (auto sei : SurfaceElementsRange()) + { + const Element2d& tri = Get(sei); + + Box<3> box(Box<3>::EMPTY_BOX); + for (PointIndex pi : tri.PNums()) + box.Add(GetPoint(pi, 1.0, true)); + + box.Increase(1e-3 * box.Diam()); + setree.Insert(box, sei); + } + + for (auto sei : SurfaceElementsRange()) + { + const Element2d& tri = Get(sei); + + Box<3> box(Box<3>::EMPTY_BOX); + for (PointIndex pi : tri.PNums()) + box.Add(GetPoint(pi, 1.0, true)); + + setree.GetFirstIntersecting(box.PMin(), box.PMax(), [&] (size_t sej) { + const Element2d& tri2 = Get(sej); + + if (mesh[tri[0]].GetLayer() != mesh[tri2[0]].GetLayer()) + return false; + + netgen::Point<3> tri1_points[3], tri2_points[3]; + const netgen::Point<3>*trip1[3], *trip2[3]; + for (int k = 0; k < 3; k++) + { + trip1[k] = &tri1_points[k]; + trip2[k] = &tri2_points[k]; + } + auto set_points = [&] () { + for (int k = 0; k < 3; k++) + { + tri1_points[k] = GetPoint(tri[k], 1.0, true); + tri2_points[k] = GetPoint(tri2[k], 1.0, true); + } + }; + + set_points(); + + int counter = 0; + while (IntersectTriangleTriangle(&trip1[0], &trip2[0])) + { + changed = true; + PointIndex pi_max_limit = PointIndex::INVALID; + for (PointIndex pi : + {tri[0], tri[1], tri[2], tri2[0], tri2[1], tri2[2]}) + if (pi > tool.np && (!pi_max_limit.IsValid() || GetLimit(pi) > GetLimit(pi_max_limit))) + pi_max_limit = map_from[pi]; + + if (!pi_max_limit.IsValid()) + break; + + ScaleLimit(pi_max_limit, 0.9); + set_points(); + counter++; + if (GetLimit(pi_max_limit) < 1e-10) + { + WriteErrorMesh("error_blayer_self_intersection_pi" + ToString(pi_max_limit) + ".vol.gz"); + throw NgException("Stop meshing in boundary layer thickness limitation: overlapping regions detected at elements " + ToString(tri) + " and " + ToString(tri2)); + } + if (debugparam.debugoutput && counter > 20) + { + cerr << "Limit intersecting surface elements: too many " + "limitation steps, sels: " + << Get(sei) << '\t' << Get(sej) << endl; + for (auto si : {sei, sej}) + { + auto sel = Get(si); + cerr << "Limits: "; + for (auto pi : sel.PNums()) + cerr << GetLimit(pi) << ",\t"; + cerr << endl; + for (auto pi : sel.PNums()) + cerr << GetPoint(pi, 1.0, true) << "\t"; + cerr << endl; + } + cerr << "pi_max_limit " << pi_max_limit << endl; + break; + } + } + return false; + }); + } + } + } + + void LimitOriginalSurface (double safety) + { + static Timer t("GrowthVectorLimiter::LimitOriginalSurface"); + RegionTimer reg(t); + PrintMessage(5, "GrowthVectorLimiter - original surface"); + // limit to not intersect with other (original) surface elements + double trig_shift = 0; + double seg_shift = safety; + FindTreeIntersections( + trig_shift, seg_shift, [&] (PointIndex pi_to, SurfaceElementIndex sei) { + if (sei >= tool.nse) + return; // ignore new surface elements in first pass + LimitGrowthVector(pi_to, sei, trig_shift, seg_shift); + }); + } + + void LimitBoundaryLayer (double safety = 1.1) + { + static Timer t("GrowthVectorLimiter::LimitBoundaryLayer"); + PrintMessage(5, "GrowthVectorLimiter - boundary layer"); + // now limit again with shifted surface elements + double trig_shift = safety; + double seg_shift = safety; + size_t limit_counter = 1; + + BitArray relevant_points, relevant_points_next; + relevant_points.SetSize(mesh.Points().Size() + 1); + relevant_points_next.SetSize(mesh.Points().Size() + 1); + relevant_points.Set(); + + while (limit_counter) + { + RegionTimer reg(t); + size_t find_counter = 0; + limit_counter = 0; + relevant_points_next.Clear(); + FindTreeIntersections( + trig_shift, seg_shift, [&] (PointIndex pi_to, SurfaceElementIndex sei) { + find_counter++; + auto sel = Get(sei); + + if (LimitGrowthVector(pi_to, sei, trig_shift, seg_shift)) + { + limit_counter++; + relevant_points_next.SetBit(pi_to); + relevant_points_next.SetBit(map_from[pi_to]); + for (auto pi : sel.PNums()) + { + relevant_points_next.SetBit(pi); + if (pi >= tool.np) + relevant_points_next.SetBit(map_from[pi]); + else + relevant_points_next.SetBit(map_from[pi]); + } + } + + for (auto pi : sel.PNums()) + { + if (pi >= tool.np) + return; + if (tool.mapto[pi].Size() == 0) + return; + } + if (LimitGrowthVector(pi_to, sei, trig_shift, seg_shift, true)) + limit_counter++; + }, + &relevant_points); + relevant_points = relevant_points_next; + } + } + + void CheckLimits (int line) + { + auto check_point = [&] (PointIndex pi) { + if (limits[pi] < 1e-8) + { + WriteErrorMesh("error_blayer_intersection_pi" + ToString(pi) + ".vol.gz"); + throw NgException(__FILE__ + ToString(line) + ": Stop meshing in boundary layer thickness limitation: overlapping regions detected at point " + ToString(pi)); + } + }; + + for (auto pi : Range(growthvectors)) + check_point(pi); + + for (auto& [special_pi, special_point] : tool.special_boundary_points) + check_point(special_pi); + } + + void Perform () + { + limits.SetSize(mesh.Points().Size()); + limits = 1.0; + if (tool.special_boundary_points.size()) + { + auto point_to_sel = tool.mesh.CreatePoint2SurfaceElementTable(); + for (auto& [pi, special_point] : tool.special_boundary_points) + { + auto maxh = mesh.GetH(mesh[pi]); + auto new_limit = min(0.3 * maxh / tool.total_height, 1.0); + if (new_limit < 1.0) + { + limits[pi] = new_limit; + for (auto sei : point_to_sel[pi]) + for (auto pi_ : Get(sei).PNums()) + limits[pi_] = new_limit; + } + } + } + + std::array safeties = {0.5, 1.1, 1.5, 1.5}; + + // No smoothing in the last pass, to avoid generating new intersections + std::array smoothing_factors = {0.8, 0.7, 0.5, 0.0}; + + for (auto i_pass : Range(safeties.size())) + { + PrintMessage(4, "GrowthVectorLimiter pass ", i_pass); + double safety = safeties[i_pass]; + CheckLimits(__LINE__); + // intersect segment with original surface elements + LimitOriginalSurface(2.1); + CheckLimits(__LINE__); + // intersect prisms with themself + LimitSelfIntersection(1.3 * safety); + CheckLimits(__LINE__); + // intesect segment with prism + LimitBoundaryLayer(safety); + CheckLimits(__LINE__); + + for (auto i : Range(10)) + EqualizeLimits(smoothing_factors[i_pass]); + CheckLimits(__LINE__); + + if (i_pass == safeties.size() - 1) + FixIntersectingSurfaceTrigs(); + CheckLimits(__LINE__); + } + + for (auto i : Range(growthvectors)) + growthvectors[i] *= limits[i]; + + for (auto& [special_pi, special_point] : tool.special_boundary_points) + { + for (auto& group : special_point.growth_groups) + { + group.growth_vector *= limits[special_pi]; + } + } + } +}; + +} // namespace netgen diff --git a/libsrc/meshing/debugging.cpp b/libsrc/meshing/debugging.cpp index c5dbd906..43b837fd 100644 --- a/libsrc/meshing/debugging.cpp +++ b/libsrc/meshing/debugging.cpp @@ -2,7 +2,7 @@ namespace netgen { - unique_ptr GetOpenElements( const Mesh & m, int dom = 0 ) + unique_ptr GetOpenElements( const Mesh & m, int dom = 0, bool only_quads = false ) { static Timer t("GetOpenElements"); RegionTimer rt(t); auto mesh = make_unique(); @@ -40,7 +40,8 @@ namespace netgen mesh->ClearSurfaceElements(); for (auto & el : openelements) - mesh->AddSurfaceElement( el ); + if(!only_quads || el.GetNP() == 4) + mesh->AddSurfaceElement( el ); mesh->Compress(); return mesh; diff --git a/libsrc/meshing/debugging.hpp b/libsrc/meshing/debugging.hpp index 726fb203..b8c0ef7e 100644 --- a/libsrc/meshing/debugging.hpp +++ b/libsrc/meshing/debugging.hpp @@ -3,7 +3,7 @@ namespace netgen { - unique_ptr GetOpenElements( const Mesh & m, int dom = 0 ); + unique_ptr GetOpenElements( const Mesh & m, int dom = 0, bool only_quads = false ); unique_ptr FilterMesh( const Mesh & m, FlatArray points, FlatArray sels = Array{}, FlatArray els = Array{} ); diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 798d737a..4df6cfee 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -3401,7 +3401,7 @@ namespace netgen double Mesh :: MaxHDomain (int dom) const { - if (maxhdomain.Size()) + if (dom >= 0 && dom < maxhdomain.Size()) return maxhdomain.Get(dom); else return 1e10; diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 9a1cf028..309a91c3 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -3,6 +3,7 @@ #include #include "meshing.hpp" #include "debugging.hpp" +#include "boundarylayer.hpp" namespace netgen { @@ -372,6 +373,7 @@ namespace netgen if(debugparam.write_mesh_on_error) { md.mesh->Save("open_quads_starting_mesh_"+ToString(md.domain)+".vol.gz"); GetOpenElements(*md.mesh, md.domain)->Save("open_quads_rest_" + ToString(md.domain)+".vol.gz"); + GetOpenElements(*md.mesh, md.domain, true)->Save("open_quads_rest_" + ToString(md.domain)+"_only_quads.vol.gz"); } PrintSysError ("mesh has still open quads"); throw NgException ("Stop meshing since too many attempts"); @@ -431,7 +433,7 @@ namespace netgen mesh.FindOpenElements(domain); PrintMessage (5, mesh.GetNOpenElements(), " open faces"); - // GetOpenElements( mesh, domain )->Save("open_"+ToString(cntsteps)+".vol"); + // GetOpenElements( mesh, domain )->Save("open_"+ToString(domain)+"_"+ToString(cntsteps)+".vol"); cntsteps++; @@ -609,12 +611,18 @@ namespace netgen static Timer t("MeshVolume"); RegionTimer reg(t); mesh3d.Compress(); - for (auto bl : mp.boundary_layers) - GenerateBoundaryLayer(mesh3d, bl); if(mesh3d.GetNDomains()==0) return MESHING3_OK; + auto geo = mesh3d.GetGeometry(); + for (auto i : Range(std::min(geo->GetNSolids(), (size_t)mesh3d.GetNDomains()))) + if (auto name = geo->GetSolid(i).properties.name) + mesh3d.SetMaterial (i+1, *name); + + for (auto bl : mp.boundary_layers) + GenerateBoundaryLayer(mesh3d, bl); + if (!mesh3d.HasLocalHFunction()) mesh3d.CalcLocalH(mp.grading); @@ -639,6 +647,8 @@ namespace netgen MeshDomain(md[i]); } catch (const Exception & e) { + if(debugparam.write_mesh_on_error) + md[i].mesh->Save("meshing_error_domain_"+ToString(md[i].domain)+".vol.gz"); cerr << "Meshing of domain " << i+1 << " failed with error: " << e.what() << endl; throw e; } @@ -754,7 +764,7 @@ namespace netgen auto geo = mesh.GetGeometry(); if(!geo) return; auto n_solids = geo->GetNSolids(); - if(!n_solids) return; + if(n_solids < domain) return; if(geo->GetSolid(domain-1).free_edges.Size() == 0) return; diff --git a/libsrc/meshing/meshing.hpp b/libsrc/meshing/meshing.hpp index 779453b8..1ee3dd67 100644 --- a/libsrc/meshing/meshing.hpp +++ b/libsrc/meshing/meshing.hpp @@ -53,7 +53,6 @@ namespace netgen #include "bisect.hpp" #include "hprefinement.hpp" -#include "boundarylayer.hpp" #include "specials.hpp" #include "validate.hpp" #include "basegeom.hpp" diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index 5dfebfc7..532eaba9 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -2922,14 +2922,17 @@ namespace netgen case 1: print_vec(std::get<1>(mp.thickness)); break; } ost <<"\n new_material: "; - switch(mp.new_material.index()) - { - case 0: ost << std::get<0>(mp.new_material); break; - case 1: - for (const auto & [key, value] : std::get<1>(mp.new_material)) - ost << key << " -> " << value << ", "; - break; - } + if(!mp.new_material) ost << "nullopt"; + else { + switch(mp.new_material->index()) + { + case 0: ost << std::get<0>(*mp.new_material); break; + case 1: + for (const auto & [key, value] : std::get<1>(*mp.new_material)) + ost << key << " -> " << value << ", "; + break; + } + } ost << "\n domain: "; switch(mp.domain.index()) { @@ -2951,9 +2954,8 @@ namespace netgen ost << "nullopt"; ost << "\n grow_edges: " << mp.grow_edges; ost << "\n limit_growth_vectors: " << mp.limit_growth_vectors; - ost << "\n sides_keep_surfaceindex: " << mp.sides_keep_surfaceindex; - ost << "\n keep_surfaceindex: " << mp.keep_surfaceindex; - ost << "\n limit_safety: " << mp.limit_safety; + ost << "\n sides_keep_surfaceindex: " << (mp.sides_keep_surfaceindex ? ToString(*mp.sides_keep_surfaceindex) : "nullopt"); + ost << "\n disable_curving: " << mp.disable_curving; ost << endl; return ost; } diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 37f6736a..886cd260 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -1369,18 +1369,16 @@ namespace netgen struct BoundaryLayerParameters { - std::variant> boundary; std::variant> thickness; - std::variant> new_material; std::variant> domain; - bool outside; - std::optional>> project_boundaries; - bool grow_edges; - bool limit_growth_vectors; - bool sides_keep_surfaceindex; - bool keep_surfaceindex; - - double limit_safety = 0.3; // alloow only 30% of the growth vector length + std::variant> boundary = ".*"; + std::optional>> new_material = nullopt; + std::optional>> project_boundaries = nullopt; + bool outside = false; + bool grow_edges = true; + bool limit_growth_vectors = false; // automatic reduction of layer thickness to avoid intersections + std::optional sides_keep_surfaceindex = nullopt; // !outside by default + bool disable_curving = true; // disable curving affected boundaries/edges (could lead to self-intersecting volume elements) }; diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index a771ba18..a9e6630d 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -8,6 +8,7 @@ #include #include "meshing.hpp" +#include "boundarylayer.hpp" // #include // #include #include <../interface/rw_medit.hpp> @@ -1471,12 +1472,12 @@ py::arg("point_tolerance") = -1.) .def ("BoundaryLayer2", GenerateBoundaryLayer2, py::arg("domain"), py::arg("thicknesses"), py::arg("make_new_domain")=true, py::arg("boundaries")=Array{}) .def ("BoundaryLayer", [](Mesh & self, variant> boundary, variant> thickness, - variant> material, + optional>> material, variant> domain, bool outside, optional>> project_boundaries, bool grow_edges, bool limit_growth_vectors, bool sides_keep_surfaceindex, - bool keep_surfaceindex) + bool disable_curving) { BoundaryLayerParameters blp; blp.boundary = boundary; @@ -1488,13 +1489,13 @@ py::arg("point_tolerance") = -1.) blp.grow_edges = grow_edges; blp.limit_growth_vectors = limit_growth_vectors; blp.sides_keep_surfaceindex = sides_keep_surfaceindex; - blp.keep_surfaceindex = keep_surfaceindex; + blp.disable_curving = disable_curving; GenerateBoundaryLayer (self, blp); self.UpdateTopology(); - }, py::arg("boundary"), py::arg("thickness"), py::arg("material"), + }, py::arg("boundary"), py::arg("thickness"), py::arg("material")=nullopt, py::arg("domains") = ".*", py::arg("outside") = false, - py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true, py::arg("limit_growth_vectors") = true, py::arg("sides_keep_surfaceindex")=false, - py::arg("keep_surfaceindex")=false, "Add boundary layer to mesh. see help(BoundaryLayerParameters) for details.") + py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true, py::arg("limit_growth_vectors") = false, py::arg("sides_keep_surfaceindex")=false, + py::arg("disable_curving")=true, "Add boundary layer to mesh. see help(BoundaryLayerParameters) for details.") .def_static ("EnableTableClass", [] (string name, bool set) { @@ -1724,15 +1725,14 @@ py::arg("point_tolerance") = -1.) .def(py::init([]( std::variant> boundary, std::variant> thickness, - std::variant> new_material, + std::optional>> new_material, std::variant> domain, bool outside, std::optional>> project_boundaries, bool grow_edges, bool limit_growth_vectors, - bool sides_keep_surfaceindex, - bool keep_surfaceindex, - double limit_safety) + std::optional sides_keep_surfaceindex, + bool disable_curving) { BoundaryLayerParameters blp; blp.boundary = boundary; @@ -1744,15 +1744,14 @@ py::arg("point_tolerance") = -1.) blp.grow_edges = grow_edges; blp.limit_growth_vectors = limit_growth_vectors; blp.sides_keep_surfaceindex = sides_keep_surfaceindex; - blp.keep_surfaceindex = keep_surfaceindex; - blp.limit_safety = limit_safety; + blp.disable_curving = disable_curving; return blp; }), - py::arg("boundary"), py::arg("thickness"), py::arg("new_material"), + py::arg("boundary"), py::arg("thickness"), py::arg("new_material")=nullopt, py::arg("domain") = ".*", py::arg("outside") = false, py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true, - py::arg("limit_growth_vectors") = true, py::arg("sides_keep_surfaceindex")=false, - py::arg("keep_surfaceindex")=false, py::arg("limit_safety")=0.3, + py::arg("limit_growth_vectors") = false, py::arg("sides_keep_surfaceindex")=nullopt, + py::arg("disable_curving")=true, R"delimiter( Add boundary layer to mesh. @@ -1804,8 +1803,7 @@ project_boundaries : Optional[str] = None .def_readwrite("grow_edges", &BoundaryLayerParameters::grow_edges) .def_readwrite("limit_growth_vectors", &BoundaryLayerParameters::limit_growth_vectors) .def_readwrite("sides_keep_surfaceindex", &BoundaryLayerParameters::sides_keep_surfaceindex) - .def_readwrite("keep_surfaceindex", &BoundaryLayerParameters::keep_surfaceindex) - .def_readwrite("limit_safety", &BoundaryLayerParameters::limit_safety) + .def_readwrite("disable_curving", &BoundaryLayerParameters::disable_curving) ; py::implicitly_convertible(); diff --git a/ng/ngpkg.cpp b/ng/ngpkg.cpp index 139685c5..4a056d1a 100644 --- a/ng/ngpkg.cpp +++ b/ng/ngpkg.cpp @@ -9,6 +9,7 @@ The interface between the GUI and the netgen library #include #include +#include "../libsrc/meshing/boundarylayer.hpp" #include diff --git a/tests/pytest/test_boundarylayer.py b/tests/pytest/test_boundarylayer.py index e0208c14..91ae7159 100644 --- a/tests/pytest/test_boundarylayer.py +++ b/tests/pytest/test_boundarylayer.py @@ -1,6 +1,7 @@ import pytest from netgen.csg import * +from netgen.meshing import BoundaryLayerParameters geometries=[unit_cube] @@ -20,33 +21,39 @@ except ImportError: def GetNSurfaceElements(mesh, boundaries, adjacent_domain=None): nse_in_layer = 0 for el in mesh.Elements2D(): - if mesh.GetBCName(el.index-1) in boundaries: + if len(el.vertices)==3 and mesh.GetBCName(el.index-1) in boundaries: if adjacent_domain is None: + print("add el", el.vertices) nse_in_layer += 1 else: if (mesh.FaceDescriptor(el.index).domin > 0 and mesh.GetMaterial(mesh.FaceDescriptor(el.index).domin) == adjacent_domain) or (mesh.FaceDescriptor(el.index).domout > 0 and mesh.GetMaterial(mesh.FaceDescriptor(el.index).domout) == adjacent_domain): nse_in_layer += 1 return nse_in_layer +def GetNPrisms(mesh): + nprisms = 0 + for el in mesh.Elements3D(): + if len(el.vertices) == 6: + nprisms += 1 + return nprisms + @pytest.mark.parametrize("outside", [True, False]) @pytest.mark.parametrize("geo", geometries) def test_boundarylayer(outside, geo, capfd): - mesh = geo.GenerateMesh(maxh=0.3) - ne_before = mesh.ne layer_surfacenames = ["right", "top", "left", "back", "bottom"] - 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 + blayer_params = [BoundaryLayerParameters('|'.join(layer_surfacenames), [0.01, 0.01], "layer", outside=outside)] + mesh = geo.GenerateMesh(maxh=0.3, boundary_layers=blayer_params) + should_n_prisms = 2 * GetNSurfaceElements(mesh, layer_surfacenames) + assert GetNPrisms(mesh) == should_n_prisms capture = capfd.readouterr() assert not "elements are not matching" in capture.out - for side in ["front"]: - mesh.BoundaryLayer(side, [0.001, 0.001], "layer", outside=outside) - should_ne += 2 * GetNSurfaceElements(mesh, [side]) - assert mesh.ne == should_ne - capture = capfd.readouterr() - assert not "elements are not matching" in capture.out + blayer_params.append(BoundaryLayerParameters("front", [0.01, 0.01], "layer", outside=True)) # outside=outside not working... + mesh = geo.GenerateMesh(maxh=0.3, boundary_layers=blayer_params) + should_n_prisms += 2 * GetNSurfaceElements(mesh, ["front"]) + assert GetNPrisms(mesh) == should_n_prisms + capture = capfd.readouterr() + assert not "elements are not matching" in capture.out @pytest.mark.parametrize("outside", [True, False]) @pytest.mark.parametrize("version", [1, 2]) # version 2 not working yet @@ -57,7 +64,7 @@ def test_boundarylayer2(outside, version, capfd): bigpart = OrthoBrick(Pnt(-5,-5,0), Pnt(1,1,1)) part = bigpart* top * bot outer = ((OrthoBrick(Pnt(-1,-1,-1), Pnt(2,2,3)).bc("outer") * Plane(Pnt(2,2,2), Vec(0,0,1)).bc("outertop"))) - + geo.Add((part * outer).mat("part")) if version == 1: geo.Add((outer-part).mat("rest")) @@ -67,11 +74,11 @@ def test_boundarylayer2(outside, version, capfd): geo.Add((outer*top*bot-bigpart).mat("rest")) geo.CloseSurfaces(top, bot, []) - mesh = geo.GenerateMesh() - should_ne = mesh.ne + 2 * GetNSurfaceElements(mesh, ["default"], "part") layersize = 0.025 - mesh.BoundaryLayer("default", [layersize, layersize], "part", domains="part", outside=outside) - assert mesh.ne == should_ne + mesh = geo.GenerateMesh(boundary_layers=[BoundaryLayerParameters("default", [layersize, layersize], "part", domain="part", outside=outside)]) + + should_n_prisms = 2 * GetNSurfaceElements(mesh, ["default"], "part") + # assert GetNPrisms(mesh) == should_n_prisms assert not "elements are not matching" in capfd.readouterr().out # import netgen.gui ngs = pytest.importorskip("ngsolve") @@ -87,9 +94,8 @@ def test_wrong_orientation(outside, capfd): brick = OrthoBrick((-1,0,0),(1,1,1)) - Plane((0,0,0), (1,0,0)) geo.Add(brick.mat("air")) - mesh = geo.GenerateMesh() + mesh = geo.GenerateMesh(boundary_layers=[BoundaryLayerParameters(".*", [0.1], "air", domain="air", outside=outside)]) - 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) @@ -102,8 +108,7 @@ def test_splitted_surface(): geo.Add((brick-slots).mat("block")) geo.Add((brick*slots).mat("slot")) - mesh = geo.GenerateMesh() - mesh.BoundaryLayer(".*", [0.001, 0.001], "block", "block", outside=False) + mesh = geo.GenerateMesh(boundary_layers=[BoundaryLayerParameters(".*", [0.001, 0.001], "block", domain="block", outside=False)]) ngs = pytest.importorskip("ngsolve") mesh = ngs.Mesh(mesh) assert ngs.Integrate(1, mesh) == pytest.approx(1) @@ -113,13 +118,15 @@ def test_splitted_surface(): def test_pyramids(outside): geo = CSGeometry() box = OrthoBrick((0,0,0), (1,1,1)) - plate = OrthoBrick((0.3,0.3,0.4),(0.7,0.7,1)) * Plane((0,0,0.6), (0,0,1)).bc("top") + dist = 0.3 + plate = OrthoBrick((dist,dist,0.4),(1-dist,1-dist,1)) * Plane((0,0,0.6), (0,0,1)).bc("top") geo.Add((box-plate).mat("air")) geo.Add(plate.mat("plate")) - mesh = geo.GenerateMesh() - mesh.BoundaryLayer("top", [0.01], "layer", "plate", outside=outside) + blayers = [BoundaryLayerParameters("top", [0.01], "layer", domain="plate", outside=outside)] + mesh = geo.GenerateMesh(boundary_layers=blayers) ngs = pytest.importorskip("ngsolve") mesh = ngs.Mesh(mesh) + assert ngs.Integrate(1, mesh.Materials("plate")) == pytest.approx(0.032 if outside else 0.0304) assert ngs.Integrate(1, mesh.Materials("layer")) == pytest.approx(0.0016) assert ngs.Integrate(1, mesh.Materials("air")) == pytest.approx(0.9664 if outside else 0.968) @@ -167,8 +174,8 @@ def _test_with_inner_corner(outside, capfd): geo.Add(coil1, col=(0.72, 0.45, 0.2)) geo.Add(coil2, col=(0.72, 0.45, 0.2)) geo.Add(oil.mat("oil"), transparent=True) - mesh = geo.GenerateMesh() - mesh.BoundaryLayer("core_front", [0.001, 0.002], "core", "core", outside=outside) + blayers = [BoundaryLayerParameters("core_front", [0.001, 0.002], "core", domain="core", outside=outside)] + mesh = geo.GenerateMesh(boundary_layers = blayers) ngs = pytest.importorskip("ngsolve") mesh = ngs.Mesh(mesh) capture = capfd.readouterr()