From 36440970fbeb57cdcacf69374336d12b0ce70ee7 Mon Sep 17 00:00:00 2001 From: "mhochsteger@cerbsim.com" Date: Tue, 1 Mar 2022 20:18:05 +0100 Subject: [PATCH] boundarylayer - some more fixes on growth vector interpolation --- libsrc/meshing/boundarylayer.cpp | 160 +++++++++++++++++++++---------- libsrc/meshing/delaunay2d.cpp | 61 +++++++----- libsrc/meshing/delaunay2d.hpp | 2 +- libsrc/meshing/topology.hpp | 1 + 4 files changed, 147 insertions(+), 77 deletions(-) diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index cd9c3597..a8728763 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -180,12 +180,12 @@ namespace netgen } } - void InterpolateSurfaceGrowthVectors(const Mesh & mesh, const BoundaryLayerParameters& blp, int fd_old, FlatArray, PointIndex> growthvectors) + void InterpolateSurfaceGrowthVectors(const Mesh & mesh, const BoundaryLayerParameters& blp, int fd_old, FlatArray, PointIndex> growthvectors, const Table & p2sel) { + auto np = mesh.GetNP(); // interpolate growth vectors at inner surface points from surrounding edge points - Array, PointIndex> delaunay_points(mesh.GetNP()); - Array p2face(mesh.GetNP()); - p2face = 0; + Array, PointIndex> delaunay_points(np); + Array pmap(np); // maps duplicated points Array surface_els; Array edge_points; @@ -195,33 +195,75 @@ namespace netgen if(!blp.surfid.Contains(facei)) continue; - p2face = 0; - edge_points.SetSize(0); surface_points.SetSize(0); surface_els.SetSize(0); + delaunay_points.SetSize(np); + pmap.SetSize(np); mesh.GetSurfaceElementsOfFace (facei, surface_els); Box<2> bbox ( Box<2>::EMPTY_BOX ); for(auto sei : surface_els) { - const auto & sel = mesh[sei]; + auto sel = mesh[sei]; for (auto i : Range(sel.GetNP())) { - auto pi = sel[i]; - if(p2face[pi] != 0) + auto & gi = sel.GeomInfo()[i]; + Point<2> p = {gi.u, gi.v}; + bbox.Add(p); + } + } + + BoxTree<2> tree(bbox); + + for(auto pi : mesh.Points().Range()) + { + auto n_surf_els = p2sel[pi].Size(); + + bool has_relevant_sel = false; + for(auto sei : p2sel[pi]) + if(mesh[sei].GetIndex()==facei) + { + has_relevant_sel = true; + break; + } + if(!has_relevant_sel) + continue; + + if(mesh[pi].Type() <= EDGEPOINT) + edge_points.Append(pi); + else + surface_points.Append(pi); + + // the same point might have different uv coordinates (closed edges for instance) + // duplicate these points for the delaunay tree + bool inserted = false; + for(auto sei : p2sel[pi]) + { + auto sel = mesh[sei]; + if(sel.GetIndex()!=facei) continue; - p2face[pi] = facei; + PointGeomInfo gi = sel.GeomInfo()[sel.PNums().Pos(pi)]; + Point<2> p = {gi.u, gi.v}; + bool found = false; + tree.GetFirstIntersecting( p, p, [&] (const auto pi_found) { return found = true; }); + if(found) + continue; - if(mesh[pi].Type() <= EDGEPOINT) - edge_points.Append(pi); + auto pi_new = pi; + if(inserted) + { + pi_new = delaunay_points.Append(p); + pmap.Append(pi); + edge_points.Append(pi_new); + } else - surface_points.Append(pi); - - auto & gi = sel.GeomInfo()[i]; - // TODO: project to plane if u,v not available? - delaunay_points[pi] = {gi.u, gi.v}; - bbox.Add(delaunay_points[pi]); + { + delaunay_points[pi] = p; + pmap[pi] = pi; + } + tree.Insert(p, pi_new); + inserted = true; } } @@ -231,10 +273,7 @@ namespace netgen DelaunayMesh dmesh( delaunay_points, bbox ); for(auto pi : edge_points) - { - p2face[pi] = 0; dmesh.AddPoint(pi); - } std::map weights; for(auto pi : surface_points) @@ -245,15 +284,13 @@ namespace netgen auto n = 1./v.Length() * v; for(auto & [pi_other, weight] : weights) { - auto t = weight * growthvectors[pi_other]; + // interpolate only the tangential part of the growth vector + auto t = weight * growthvectors[pmap[pi_other]]; t -= (t * n) * n; v += t; } } - } - - } void GenerateBoundaryLayer(Mesh& mesh, const BoundaryLayerParameters& blp) @@ -277,21 +314,17 @@ namespace netgen if(!blp.outside) domains.Invert(); + auto& meshtopo = mesh.GetTopology(); + meshtopo.SetBuildVertex2Element(true); mesh.UpdateTopology(); bool have_single_segments = HaveSingleSegments(mesh); Array segments; if(have_single_segments) - { segments = BuildSegments(mesh); - mesh.SetNextMajorTimeStamp(); - mesh.UpdateTopology(); - } else segments = mesh.LineSegments(); - auto& meshtopo = mesh.GetTopology(); - int np = mesh.GetNP(); int ne = mesh.GetNE(); int nse = mesh.GetNSE(); @@ -339,24 +372,50 @@ namespace netgen } } - // mark points for remapping - for(const auto& sel : mesh.SurfaceElements()) - { - auto n = surfacefacs[sel.GetIndex()] * getSurfaceNormal(sel); - if(n.Length2() != 0.) - { - for(auto pi : sel.PNums()) - { - auto & np = growthvectors[pi]; - if(np.Length() == 0) { np = n; continue; } - auto npn = np * n; - auto npnp = np * np; - auto nn = n * n; - if(nn-npn*npn/npnp == 0) { np = n; continue; } - np += (nn - npn)/(nn - npn*npn/npnp) * (n - npn/npnp * np); - } - } - } + const auto & p2sel = mesh.CreatePoint2SurfaceElementTable(); + + 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(!blp.surfid.Contains(facei)) + continue; + + auto n = surfacefacs[sel.GetIndex()] * getSurfaceNormal(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 + auto & np = growthvectors[pi]; + for(auto & [facei, n] : normals) + { + if(np.Length() == 0) { np = n; continue; } + auto npn = np * n; + auto npnp = np * np; + auto nn = n * n; + if(nn-npn*npn/npnp == 0) { np = n; continue; } + np += (nn - npn)/(nn - npn*npn/npnp) * (n - npn/npnp * np); + } + } // Bit array to keep track of segments already processed BitArray segs_done(nseg+1); @@ -551,7 +610,7 @@ namespace netgen } if(interpolate_growth_vectors) - InterpolateSurfaceGrowthVectors(mesh, blp, fd_old, growthvectors); + InterpolateSurfaceGrowthVectors(mesh, blp, fd_old, growthvectors, p2sel); // insert new points for (PointIndex pi = 1; pi <= np; pi++) @@ -889,6 +948,7 @@ namespace netgen } mesh.GetTopology().ClearEdges(); mesh.UpdateTopology(); + mesh.SetGeometry(nullptr); } void AddDirection( Vec<3> & a, Vec<3> b ) diff --git a/libsrc/meshing/delaunay2d.cpp b/libsrc/meshing/delaunay2d.cpp index 02aba5de..07732695 100644 --- a/libsrc/meshing/delaunay2d.cpp +++ b/libsrc/meshing/delaunay2d.cpp @@ -190,26 +190,7 @@ namespace netgen if(definitive_overlapping_trig==-1) { - Mesh m; - m.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 0)); - for(auto pi : points.Range()) - m.AddPoint(P3(points[pi])); - - for (DelaunayTrig & trig : trigs) - { - if (trig[0] < 0) continue; - - Vec<3> n = Cross (P3(points[trig[1]])-P3(points[trig[0]]), - P3(points[trig[2]])-P3(points[trig[0]])); - if (n(2) < 0) Swap (trig[1], trig[2]); - - Element2d el(trig[0], trig[1], trig[2]); - el.SetIndex (1); - m.AddSurfaceElement (el); - } - m.Compress(); - m.AddPoint(P3(points[pi_new])); - m.Save("error.vol.gz"); + // GetMesh(pi_new)->Save("error.vol.gz"); throw Exception("point not in any circle "+ ToString(pi_new)); } @@ -299,15 +280,18 @@ namespace netgen auto pi_last = *points.Range().end()-3; for(auto edge : edges) { + auto v0 = points[edge[0]] - p; + auto v1 = points[edge[1]] - p; + v0.Normalize(); + v1.Normalize(); + double angle = acos(v0*v1); for(PointIndex pi : {edge[0], edge[1]}) { - if(pi>=pi_last) + if(pi>=pi_last) continue; - if(weights.count(pi)) - continue; - double weight = 1.0/(eps+Dist(p, points[pi])); - sum += weight; - weights[pi] = weight; + double weight = angle/(eps+Dist(p, points[pi])); + sum += weight; + weights[pi] += weight; } } double isum = 1.0/sum; @@ -336,6 +320,31 @@ namespace netgen tree->DeleteElement (j); } + unique_ptr DelaunayMesh::GetMesh(PointIndex pi_new) + { + auto mesh = make_unique(); + Mesh & m = *mesh; + m.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 0)); + for(auto pi : points.Range()) + m.AddPoint(P3(points[pi])); + + for (DelaunayTrig & trig : trigs) + { + if (trig[0] < 0) continue; + + Vec<3> n = Cross (P3(points[trig[1]])-P3(points[trig[0]]), + P3(points[trig[2]])-P3(points[trig[0]])); + if (n(2) < 0) Swap (trig[1], trig[2]); + + Element2d el(trig[0], trig[1], trig[2]); + el.SetIndex (1); + m.AddSurfaceElement (el); + } + m.Compress(); + m.AddPoint(P3(points[pi_new])); + return mesh; + } + ostream & operator<< (ostream & ost, DelaunayTrig trig) { ost << trig[0] << "-" << trig[1] << "-" << trig[2] << endl; diff --git a/libsrc/meshing/delaunay2d.hpp b/libsrc/meshing/delaunay2d.hpp index dbc7524a..40f3b000 100644 --- a/libsrc/meshing/delaunay2d.hpp +++ b/libsrc/meshing/delaunay2d.hpp @@ -68,7 +68,7 @@ namespace netgen void CalcWeights( PointIndex pi_new, std::map & weights ); void AddPoint( PointIndex pi_new ); Array & GetElements() { return trigs; } - + unique_ptr GetMesh(PointIndex pi_new); // for debugging purposes }; } // namespace netgen diff --git a/libsrc/meshing/topology.hpp b/libsrc/meshing/topology.hpp index 9447c545..ce753a4a 100644 --- a/libsrc/meshing/topology.hpp +++ b/libsrc/meshing/topology.hpp @@ -186,6 +186,7 @@ public: { return vert2element[vnr]; } void GetVertexSurfaceElements( int vnr, Array& elements ) const; + const auto & GetVertexSurfaceElements( ) const { return vert2surfelement; } FlatArray GetVertexSurfaceElements(PointIndex vnr) const { return vert2surfelement[vnr]; }