From 016b1692e2b10ebcc8ac24755c5fdd40a7bb0711 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Tue, 1 Mar 2022 13:23:06 +0100 Subject: [PATCH 1/5] fix point type of geo vertices (FIXEDPOINT) -> locked points --- libsrc/meshing/basegeom.cpp | 2 ++ libsrc/meshing/boundarylayer.cpp | 12 ++++++++++-- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp index ca1b7651..2b221de6 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -529,6 +529,8 @@ namespace netgen tree.Insert(mesh[pi], pi); vert2meshpt[vert->GetHash()] = pi; mesh[pi].Singularity(vert->properties.hpref); + mesh[pi].SetType(FIXEDPOINT); + mesh.AddLockedPoint(pi); if(vert->properties.name) { diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index 8ae125ab..4c1bd8c0 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -276,7 +276,11 @@ namespace netgen bool have_single_segments = HaveSingleSegments(mesh); Array segments; if(have_single_segments) + { segments = BuildSegments(mesh); + mesh.SetNextMajorTimeStamp(); + mesh.UpdateTopology(); + } else segments = mesh.LineSegments(); @@ -492,25 +496,29 @@ namespace netgen if(seg.edgenr-1 == edgenr && mesh[seg[0]].Type() == FIXEDPOINT) { points.Append(seg[0]); + points.Append(seg[1]); break; } } while(true) { + bool point_found = false; for(auto si : meshtopo.GetVertexSegments(points.Last())) { const auto& seg = mesh[si]; if(seg.edgenr-1 != edgenr) continue; - if(seg[0] == points.Last() && - (points.Size() < 2 || points[points.Size()-2] !=seg[1])) + if(seg[0] == points.Last() && points[points.Size()-2] !=seg[1]) { points.Append(seg[1]); + point_found = true; break; } } if(mesh[points.Last()].Type() == FIXEDPOINT) break; + if(!point_found) + throw Exception(string("Could not find connected list of line segements for edge ") + edgenr); } // tangential part of growth vectors From 13a0b78e264f89d807d3a0b454ad467301747897 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Tue, 1 Mar 2022 14:34:18 +0100 Subject: [PATCH 2/5] interpolate only tangential part of growth vector --- libsrc/meshing/boundarylayer.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index 4c1bd8c0..4d1276db 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -241,8 +241,13 @@ namespace netgen { dmesh.AddPoint(pi, &weights); auto & v = growthvectors[pi]; + auto n = 1./v.Length() * v; for(auto & [pi_other, weight] : weights) - v += weight * growthvectors[pi_other]; + { + auto t = weight * growthvectors[pi_other]; + t -= (t * n) * t; + v += t; + } } } @@ -341,8 +346,6 @@ namespace netgen { for(auto pi : sel.PNums()) { - if(interpolate_growth_vectors && mesh[pi].Type() >= SURFACEPOINT) - continue; auto & np = growthvectors[pi]; if(np.Length() == 0) { np = n; continue; } auto npn = np * n; @@ -436,8 +439,6 @@ namespace netgen for(auto i : Range(sel.PNums())) { auto pi = sel.PNums()[i]; - if(interpolate_growth_vectors && mesh[pi].Type() >= SURFACEPOINT) - continue; if(growthvectors[pi].Length2() == 0.) continue; auto next = sel.PNums()[(i+1)%sel.GetNV()]; From e8c9d8e1fc6da94a5aece638f2cc864727ea6aee Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Tue, 1 Mar 2022 14:56:01 +0100 Subject: [PATCH 3/5] really interpolate only tangential part... --- libsrc/meshing/boundarylayer.cpp | 5 +++-- libsrc/meshing/delaunay2d.cpp | 5 +---- libsrc/meshing/delaunay2d.hpp | 2 +- 3 files changed, 5 insertions(+), 7 deletions(-) diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index 4d1276db..cd9c3597 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -239,13 +239,14 @@ namespace netgen std::map weights; for(auto pi : surface_points) { - dmesh.AddPoint(pi, &weights); + dmesh.CalcIntersecting(pi); + dmesh.CalcWeights(pi, weights); auto & v = growthvectors[pi]; auto n = 1./v.Length() * v; for(auto & [pi_other, weight] : weights) { auto t = weight * growthvectors[pi_other]; - t -= (t * n) * t; + t -= (t * n) * n; v += t; } } diff --git a/libsrc/meshing/delaunay2d.cpp b/libsrc/meshing/delaunay2d.cpp index 558dfcef..02aba5de 100644 --- a/libsrc/meshing/delaunay2d.cpp +++ b/libsrc/meshing/delaunay2d.cpp @@ -315,15 +315,12 @@ namespace netgen weight *= isum; } - void DelaunayMesh::AddPoint( PointIndex pi_new, std::map * weights ) + void DelaunayMesh::AddPoint( PointIndex pi_new) { static Timer t("AddPoint"); RegionTimer reg(t); CalcIntersecting(pi_new); - if(weights) - CalcWeights(pi_new, *weights); - for (int j : intersecting) { UnsetNeighbours(j); diff --git a/libsrc/meshing/delaunay2d.hpp b/libsrc/meshing/delaunay2d.hpp index 3e451f13..dbc7524a 100644 --- a/libsrc/meshing/delaunay2d.hpp +++ b/libsrc/meshing/delaunay2d.hpp @@ -66,7 +66,7 @@ namespace netgen void CalcIntersecting( PointIndex pi_new ); void CalcWeights( PointIndex pi_new, std::map & weights ); - void AddPoint( PointIndex pi_new, std::map * weights = nullptr ); + void AddPoint( PointIndex pi_new ); Array & GetElements() { return trigs; } }; From 36440970fbeb57cdcacf69374336d12b0ce70ee7 Mon Sep 17 00:00:00 2001 From: "mhochsteger@cerbsim.com" Date: Tue, 1 Mar 2022 20:18:05 +0100 Subject: [PATCH 4/5] 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]; } From 9730a383fdf94d0bf8f1a857d85ed46612a0c123 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Wed, 2 Mar 2022 11:34:02 +0100 Subject: [PATCH 5/5] geo vertices as pointelements not locked points --- libsrc/meshing/basegeom.cpp | 12 ++++-------- libsrc/meshing/meshclass.cpp | 2 ++ libsrc/meshing/meshfunc.cpp | 7 +++++-- 3 files changed, 11 insertions(+), 10 deletions(-) diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp index 2b221de6..ff994594 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -530,15 +530,11 @@ namespace netgen vert2meshpt[vert->GetHash()] = pi; mesh[pi].Singularity(vert->properties.hpref); mesh[pi].SetType(FIXEDPOINT); - mesh.AddLockedPoint(pi); - if(vert->properties.name) - { - Element0d el(pi, pi); - el.name = vert->properties.GetName(); - mesh.SetCD3Name(pi, el.name); - mesh.pointelements.Append (el); - } + Element0d el(pi, pi); + el.name = vert->properties.GetName(); + mesh.SetCD3Name(pi, el.name); + mesh.pointelements.Append (el); } for(auto & vert : vertices) diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index b381f760..89375e47 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -2352,6 +2352,8 @@ namespace netgen for (int i = 0; i < lockedpoints.Size(); i++) points[lockedpoints[i]].SetType(FIXEDPOINT); + for(const auto& pointel : pointelements) + points[pointel.pnum].SetType(FIXEDPOINT); /* for (i = 0; i < openelements.Size(); i++) diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 700137e1..07c2ed58 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -102,7 +102,7 @@ namespace netgen // mark locked/fixed points for each domain TODO: domain bounding box to add only relevant points? for(auto pi : mesh.LockedPoints()) for(auto i : Range(ret)) - ipmap[i][pi] = 1; + ipmap[i][pi] = 2; // add used points to domain mesh, build point mapping for(auto i : Range(ret)) @@ -113,7 +113,10 @@ namespace netgen for(auto pi : Range(ipmap[i])) if(ipmap[i][pi]) { - auto pi_new = m.AddPoint( mesh[pi] ); + const auto& mp = mesh[pi]; + auto pi_new = m.AddPoint( mp, mp.GetLayer(), mp.Type() ); + if(ipmap[i][pi] == 2) + mesh.AddLockedPoint(pi_new); ipmap[i][pi] = pi_new; pmap.Append( pi ); }