From a76d4079793eb5f24067ae2a7690dcd1948470e3 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Mon, 28 Oct 2019 19:58:35 +0100 Subject: [PATCH] implement functionality to restrict meshsize in base class not yet used in any derived geometry --- libsrc/meshing/basegeom.cpp | 111 +++++++++++++++++++++++++++++++++++- libsrc/meshing/basegeom.hpp | 18 +++++- 2 files changed, 126 insertions(+), 3 deletions(-) diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp index 7f160bac..7b772e5d 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -10,6 +10,81 @@ namespace netgen GeometryRegister :: ~GeometryRegister() { ; } + void GeometryFace :: RestrictHTrig(Mesh& mesh, + const PointGeomInfo& gi0, + const PointGeomInfo& gi1, + const PointGeomInfo& gi2, + const MeshingParameters& mparam, + int depth, double h) const + { + auto p0 = GetPoint(gi0); + auto p1 = GetPoint(gi1); + auto p2 = GetPoint(gi2); + auto longest = (p0-p1).Length(); + int cutedge = 2; + if(auto len = (p0-p2).Length(); len > longest) + { + longest = len; + cutedge = 1; + } + if(auto len = (p1-p2).Length(); len > longest) + { + longest = len; + cutedge = 0; + } + PointGeomInfo gi_mid; + gi_mid.u = (gi0.u + gi1.u + gi2.u)/3; + gi_mid.v = (gi0.v + gi1.v + gi2.v)/3; + + if(depth % 3 == 0) + { + double curvature = 0.; + curvature = max2(max2(curvature, GetCurvature(gi_mid)), + max3(GetCurvature(gi0), GetCurvature(gi1), + GetCurvature(gi2))); + if(curvature < 1e-3) + return; + double kappa = curvature * mparam.curvaturesafety; + h = mparam.maxh * kappa < 1 ? mparam.maxh : 1./kappa; + if(h < 1e-4 * longest) + return; + } + + if(h < longest && depth < 10) + { + if(cutedge == 0) + { + PointGeomInfo gi_m; + gi_m.u = 0.5 * (gi1.u + gi2.u); + gi_m.v = 0.5 * (gi1.v + gi2.v); + RestrictHTrig(mesh, gi_m, gi2, gi0, mparam, depth+1, h); + RestrictHTrig(mesh, gi_m, gi0, gi1, mparam, depth+1, h); + } + else if(cutedge == 1) + { + PointGeomInfo gi_m; + gi_m.u = 0.5 * (gi0.u + gi2.u); + gi_m.v = 0.5 * (gi0.v + gi2.v); + RestrictHTrig(mesh, gi_m, gi1, gi2, mparam, depth+1, h); + RestrictHTrig(mesh, gi_m, gi0, gi1, mparam, depth+1, h); + } + else if(cutedge == 2) + { + PointGeomInfo gi_m; + gi_m.u = 0.5 * (gi0.u + gi1.u); + gi_m.v = 0.5 * (gi0.v + gi1.v); + RestrictHTrig(mesh, gi_m, gi1, gi2, mparam, depth+1, h); + RestrictHTrig(mesh, gi_m, gi2, gi0, mparam, depth+1, h); + } + } + else + { + auto pmid = GetPoint(gi_mid); + for(const auto& p : {p0, p1, p2, pmid}) + mesh.RestrictLocalH(p, h); + } + } + void NetgenGeometry :: Analyse(Mesh& mesh, const MeshingParameters& mparam) const { @@ -21,7 +96,41 @@ namespace netgen mparam.grading); if(mparam.uselocalh) - RestrictLocalMeshsize(mesh, mparam); + { + double eps = 1e-12 * bounding_box.Diam(); + + // restrict meshsize on edges + for(const auto & edge : edges) + { + auto length = edge->GetLength(); + // skip very short edges + if(length < eps) + continue; + static constexpr int npts = 20; + // restrict mesh size based on edge length + for(auto i : Range(npts+1)) + mesh.RestrictLocalH(edge->GetPoint(double(i)/npts), length/mparam.segmentsperedge); + + // restrict mesh size based on edge curvature + double t = 0.; + auto p_old = edge->GetPoint(t); + while(t < 1.-eps) + { + t += edge->CalcStep(t, 1./mparam.curvaturesafety); + if(t < 1.) + { + auto p = edge->GetPoint(t); + auto dist = (p-p_old).Length(); + mesh.RestrictLocalH(p, dist); + p_old = p; + } + } + } + + // restrict meshsize on faces + for(const auto& face : faces) + face->RestrictH(mesh, mparam); + } mesh.LoadLocalMeshSize(mparam.meshsizefilename); } diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index 3835937c..6ef03aef 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -28,6 +28,8 @@ namespace netgen virtual const GeometryVertex& GetEndVertex() const = 0; virtual double GetLength() const = 0; virtual Point<3> GetPoint(double t) const = 0; + // Calculate parameter step respecting edges sag value + virtual double CalcStep(double t, double sag) const = 0; virtual bool OrientedLikeGlobal() const = 0; virtual size_t GetHash() const = 0; }; @@ -42,10 +44,24 @@ namespace netgen // Project point using geo info. Fast if point is close to // parametrization in geo info. virtual bool ProjectPointGI(Point<3>& p, PointGeomInfo& gi) const =0; + virtual Point<3> GetPoint(const PointGeomInfo& gi) const = 0; virtual void CalcEdgePointGI(const GeometryEdge& edge, double t, EdgePointGeomInfo& egi) const = 0; virtual Box<3> GetBoundingBox() const = 0; + + // Get curvature in point from local coordinates in PointGeomInfo + virtual double GetCurvature(const PointGeomInfo& gi) const = 0; + + virtual void RestrictH(Mesh& mesh, const MeshingParameters& mparam) const = 0; + + protected: + void RestrictHTrig(Mesh& mesh, + const PointGeomInfo& gi0, + const PointGeomInfo& gi1, + const PointGeomInfo& gi2, + const MeshingParameters& mparam, + int depth = 0, double h = 0.) const; }; class DLL_HEADER NetgenGeometry @@ -76,8 +92,6 @@ namespace netgen virtual Mesh::GEOM_TYPE GetGeomType() const { return Mesh::NO_GEOM; } virtual void Analyse(Mesh& mesh, const MeshingParameters& mparam) const; - virtual void RestrictLocalMeshsize(Mesh& mesh, - const MeshingParameters& mparam) const {} virtual void FindEdges(Mesh& mesh, const MeshingParameters& mparam) const; virtual void MeshSurface(Mesh& mesh, const MeshingParameters& mparam) const; virtual void OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const;