From fbf6d92895bb52f37f59d01f90f9e0e4c146d2ed Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Thu, 31 Oct 2019 13:34:40 +0100 Subject: [PATCH 1/5] clean up geometry interface, fix if number of subdivided edges was not correct --- libsrc/meshing/basegeom.cpp | 24 ++++++------- libsrc/meshing/basegeom.hpp | 69 +++++++++++++++++++++++++++++++------ 2 files changed, 71 insertions(+), 22 deletions(-) diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp index 1adb07dd..0e5fa405 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -10,14 +10,6 @@ namespace netgen GeometryRegister :: ~GeometryRegister() { ; } - Array> GeometryEdge :: GetEquidistantPointArray(size_t npoints) const - { - Array> pts(npoints); - for(auto i : Range(npoints)) - pts[i] = GetPoint(double(i)/(npoints-1)); - return pts; - } - void GeometryFace :: RestrictHTrig(Mesh& mesh, const PointGeomInfo& gi0, const PointGeomInfo& gi1, @@ -148,6 +140,8 @@ namespace netgen static Timer t1("MeshEdges"); RegionTimer regt(t1); static Timer tdivide("Divide Edges"); static Timer tdivedgesections("Divide edge sections"); + const char* savetask = multithread.task; + multithread.task = "Mesh Edges"; // create face descriptors and set bc names mesh.SetNBCNames(faces.Size()); @@ -196,11 +190,15 @@ namespace netgen double hvalue[divide_edge_sections+1]; hvalue[0] = 0; + Point<3> old_pt = edge->GetPoint(0.); // calc local h for edge tdivedgesections.Start(); - auto edgepts = edge->GetEquidistantPointArray(divide_edge_sections+1); - for(auto i : Range(edgepts.Size()-1)) - hvalue[i+1] = hvalue[i] + 1./mesh.GetH(edgepts[i+1]) * (edgepts[i+1]-edgepts[i]).Length(); + for(auto i : Range(divide_edge_sections)) + { + auto pt = edge->GetPoint(double(i+1)/divide_edge_sections); + hvalue[i+1] = hvalue[i] + 1./mesh.GetH(pt) * (pt-old_pt).Length(); + old_pt = pt; + } int nsubedges = max2(1, int(floor(hvalue[divide_edge_sections]+0.5))); tdivedgesections.Stop(); mps.SetSize(nsubedges-1); @@ -235,7 +233,7 @@ namespace netgen cout << "CORRECTED" << endl; mps.SetSize (nsubedges-2); params.SetSize (nsubedges); - params[nsubedges] = 1.; + params[nsubedges-1] = 1.; } tdivide.Stop(); // ----------- Add Points to mesh and create segments ----- @@ -292,6 +290,8 @@ namespace netgen } } } + mesh.CalcSurfacesOfNode(); + multithread.task = savetask; } void NetgenGeometry :: MeshSurface(Mesh& mesh, diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index c854b490..b9e7bb16 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -32,7 +32,15 @@ namespace netgen virtual double CalcStep(double t, double sag) const = 0; virtual bool OrientedLikeGlobal() const = 0; virtual size_t GetHash() const = 0; - virtual Array> GetEquidistantPointArray(size_t npoints) const; + virtual void ProjectPoint(Point<3>& p, EdgePointGeomInfo& gi) const = 0; + virtual void PointBetween(const Point<3>& p1, + const Point<3>& p2, + double secpoint, + const EdgePointGeomInfo& gi1, + const EdgePointGeomInfo& gi2, + Point<3>& newp, + EdgePointGeomInfo& newgi) const = 0; + virtual Vec<3> GetTangent(const Point<3>& p) const = 0; }; class GeometryFace @@ -42,9 +50,11 @@ namespace netgen virtual size_t GetNBoundaries() const = 0; virtual Array> GetBoundary(size_t index) const = 0; virtual string GetName() const { return "default"; } + virtual void Project(Point<3>& p) const = 0; // 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 bool CalcPointGeomInfo(const Point<3>& p, PointGeomInfo& gi) const = 0; virtual Point<3> GetPoint(const PointGeomInfo& gi) const = 0; virtual void CalcEdgePointGI(const GeometryEdge& edge, double t, @@ -55,6 +65,21 @@ namespace netgen virtual double GetCurvature(const PointGeomInfo& gi) const = 0; virtual void RestrictH(Mesh& mesh, const MeshingParameters& mparam) const = 0; + virtual Vec<3> GetNormal(const Point<3>& p) const + { + return {0.,0.,1.}; + } + virtual Vec<3> GetNormal(const Point<3>& p, const PointGeomInfo& gi) const + { + return GetNormal(p); + } + virtual void PointBetween(const Point<3>& p1, + const Point<3>& p2, + double secpoint, + const PointGeomInfo& gi1, + const PointGeomInfo& gi2, + Point<3>& newp, + PointGeomInfo& newgi) const = 0; protected: void RestrictHTrig(Mesh& mesh, @@ -100,21 +125,32 @@ namespace netgen virtual void FinalizeMesh(Mesh& mesh) const {} virtual void ProjectPoint (int surfind, Point<3> & p) const - { } - virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p) const { } - virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, EdgePointGeomInfo& gi) const - { ProjectPointEdge(surfind, surfind2, p); } + { + faces[surfind-1]->Project(p); + } - virtual bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const {return false;} + virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p) const + { + throw Exception("In ProjectPointEdge of basegeometry"); + } + virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, EdgePointGeomInfo& gi) const + { + edges[gi.edgenr]->ProjectPoint(p, gi); + } + + virtual bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const + { + return faces[surfind-1]->CalcPointGeomInfo(p3, gi); + } virtual bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const { - throw Exception("ProjectPointGI not overloaded in class" + Demangle(typeid(*this).name())); + return faces[surfind-1]->ProjectPointGI(p, gi); } virtual Vec<3> GetNormal(int surfind, const Point<3> & p) const - { return {0.,0.,1.}; } + { return faces[surfind-1]->GetNormal(p); } virtual Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const - { return GetNormal(surfind, p); } + { return faces[surfind-1]->GetNormal(p, gi); } [[deprecated]] void GetNormal(int surfind, const Point<3> & p, Vec<3> & n) const { @@ -129,6 +165,11 @@ namespace netgen Point<3> & newp, PointGeomInfo & newgi) const { + if(faces.Size()) + { + faces[surfi-1]->PointBetween(p1, p2, secpoint, gi1, gi2, newp, newgi); + return; + } newp = p1 + secpoint * (p2-p1); } @@ -140,13 +181,21 @@ namespace netgen Point<3> & newp, EdgePointGeomInfo & newgi) const { + if(edges.Size()) + { + edges[ap1.edgenr]->PointBetween(p1, p2, secpoint, + ap1, ap2, newp, newgi); + return; + } newp = p1+secpoint*(p2-p1); } virtual Vec<3> GetTangent(const Point<3> & p, int surfi1, int surfi2, const EdgePointGeomInfo & egi) const - { throw Exception("Call GetTangent of " + Demangle(typeid(*this).name())); } + { + return edges[egi.edgenr]->GetTangent(p); + } virtual size_t GetEdgeIndex(const GeometryEdge& edge) const { From 249d78508433a6e03535a3320f6d62cd4aa0bae3 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Thu, 31 Oct 2019 15:17:28 +0100 Subject: [PATCH 2/5] progress bars for find edges, mesh surface,... --- libsrc/meshing/basegeom.cpp | 9 ++++++++- libsrc/meshing/improve3.cpp | 12 ++++++------ libsrc/meshing/meshfunc.cpp | 10 +++++++--- libsrc/meshing/smoothing3.cpp | 10 +++++----- 4 files changed, 26 insertions(+), 15 deletions(-) diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp index 0e5fa405..eae2031e 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -171,6 +171,7 @@ namespace netgen auto boundary = face.GetBoundary(facebndnr); for(auto enr : Range(boundary)) { + multithread.percent = 100. * ((double(enr)/boundary.Size() + facebndnr)/face.GetNBoundaries() + facenr)/faces.Size(); const auto& oriented_edge = *boundary[enr]; auto edgenr = GetEdgeIndex(oriented_edge); const auto& edge = edges[edgenr]; @@ -298,10 +299,13 @@ namespace netgen const MeshingParameters& mparam) const { static Timer t1("Surface Meshing"); RegionTimer regt(t1); + const char* savetask = multithread.task; + multithread.task = "Mesh Surface"; Array glob2loc(mesh.GetNP()); for(auto k : Range(faces)) { + multithread.percent = 100. * k/faces.Size(); const auto& face = *faces[k]; auto bb = face.GetBoundingBox(); bb.Increase(bb.Diam()/10); @@ -354,6 +358,7 @@ namespace netgen mesh.SurfaceElements()[i].SetIndex(k+1); } } + multithread.task = savetask; } void NetgenGeometry :: OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const @@ -366,9 +371,11 @@ namespace netgen auto meshopt = MeshOptimize2d(mesh); for(auto i : Range(mparam.optsteps2d)) { - PrintMessage(2, "Optimization step ", i); + PrintMessage(3, "Optimization step ", i); + int innerstep = 0; for(auto optstep : mparam.optimize2d) { + multithread.percent = 100. * (double(innerstep++)/mparam.optimize2d.size() + i)/mparam.optsteps2d; switch(optstep) { case 's': diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp index 7138e0ab..d9cc3534 100644 --- a/libsrc/meshing/improve3.cpp +++ b/libsrc/meshing/improve3.cpp @@ -183,7 +183,7 @@ void MeshOptimize3d :: CombineImproveSequential (Mesh & mesh, // mesh.CalcSurfacesOfNode (); const char * savetask = multithread.task; - multithread.task = "Combine Improve"; + multithread.task = "Optimize Volume: Combine Improve"; double totalbad = 0; @@ -435,7 +435,7 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh, // mesh.CalcSurfacesOfNode (); const char * savetask = multithread.task; - multithread.task = "Combine Improve"; + multithread.task = "Optimize Volume: Combine Improve"; tbad.Start(); @@ -712,7 +712,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, Array elerrs(ne); const char * savetask = multithread.task; - multithread.task = "Split Improve"; + multithread.task = "Optimize Volume: Split Improve"; PrintMessage (3, "SplitImprove"); (*testout) << "start SplitImprove" << "\n"; @@ -826,7 +826,7 @@ void MeshOptimize3d :: SplitImproveSequential (Mesh & mesh, illegaltet.Clear(); const char * savetask = multithread.task; - multithread.task = "Split Improve"; + multithread.task = "Optimize Volume: Split Improve"; PrintMessage (3, "SplitImprove"); (*testout) << "start SplitImprove" << "\n"; @@ -1121,7 +1121,7 @@ void MeshOptimize3d :: SwapImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal, (*testout) << "\n" << "Start SwapImprove" << endl; const char * savetask = multithread.task; - multithread.task = "Swap Improve"; + multithread.task = "Optimize Volume: Swap Improve"; // mesh.CalcSurfacesOfNode (); @@ -2617,7 +2617,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, (*testout) << "\n" << "Start SwapImprove" << endl; const char * savetask = multithread.task; - multithread.task = "Swap Improve"; + multithread.task = "Optimize Volume: Swap Improve"; INDEX_3_HASHTABLE faces(mesh.GetNOpenElements()/3 + 2); if (goal == OPT_CONFORM) diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 4b92ed8f..f21e0812 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -646,6 +646,8 @@ namespace netgen { static Timer t("OptimizeVolume"); RegionTimer reg(t); RegionTaskManager rtm(mp.parallel_meshing ? mp.nthreads : 0); + const char* savetask = multithread.task; + multithread.task = "Optimize Volume"; int i; @@ -663,7 +665,7 @@ namespace netgen */ mesh3d.CalcSurfacesOfNode(); - for (i = 1; i <= mp.optsteps3d; i++) + for (auto i : Range(mp.optsteps3d)) { if (multithread.terminate) break; @@ -672,12 +674,13 @@ namespace netgen // teterrpow = mp.opterrpow; // for (size_t j = 1; j <= strlen(mp.optimize3d); j++) - for (size_t j = 1; j <= mp.optimize3d.length(); j++) + for (auto j : Range(mp.optimize3d.size())) { + multithread.percent = 100.* (double(j)/mp.optimize3d.size() + i)/mp.optsteps3d; if (multithread.terminate) break; - switch (mp.optimize3d[j-1]) + switch (mp.optimize3d[j]) { case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break; case 'd': optmesh.SplitImprove(mesh3d); break; @@ -698,6 +701,7 @@ namespace netgen MeshQuality3d (mesh3d); } + multithread.task = savetask; return MESHING3_OK; } diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp index de0b6cb5..e84a18ba 100644 --- a/libsrc/meshing/smoothing3.cpp +++ b/libsrc/meshing/smoothing3.cpp @@ -1164,7 +1164,7 @@ void Mesh :: ImproveMesh (const CSG eometry & geometry, OPTIMIZEGOAL goal) } const char * savetask = multithread.task; - multithread.task = "Smooth Mesh"; + multithread.task = "Optimize Volume: Smooth Mesh"; TABLE surfelementsonpoint(points.Size()); @@ -1398,7 +1398,7 @@ void Mesh :: ImproveMeshSequential (const MeshingParameters & mp, OPTIMIZEGOAL g const char * savetask = multithread.task; - multithread.task = "Smooth Mesh"; + multithread.task = "Optimize Volume: Smooth Mesh"; for (PointIndex pi : points.Range()) if ( (*this)[pi].Type() == INNERPOINT ) @@ -1524,7 +1524,7 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) } const char * savetask = multithread.task; - multithread.task = "Smooth Mesh"; + multithread.task = "Optimize Volume: Smooth Mesh"; topt.Start(); int counter = 0; @@ -1659,7 +1659,7 @@ void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp, const char * savetask = multithread.task; - multithread.task = "Smooth Mesh Jacobian"; + multithread.task = "Optimize Volume: Smooth Mesh Jacobian"; // for (PointIndex pi = points.Begin(); i < points.End(); pi++) for (PointIndex pi : points.Range()) @@ -1815,7 +1815,7 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp, const char * savetask = multithread.task; - multithread.task = "Smooth Mesh Jacobian"; + multithread.task = "Optimize Volume: Smooth Mesh Jacobian"; // for (PointIndex pi = points.Begin(); pi <= points.End(); pi++) for (PointIndex pi : points.Range()) From 1e3ed047db2a79c8ae9542c2d13d92cef60809a0 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Thu, 31 Oct 2019 15:25:47 +0100 Subject: [PATCH 3/5] progress for analyse geometry --- libsrc/meshing/basegeom.cpp | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp index eae2031e..d9542149 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -98,10 +98,14 @@ namespace netgen if(mparam.uselocalh) { double eps = 1e-12 * bounding_box.Diam(); + const char* savetask = multithread.task; + multithread.task = "Analyse Edges"; // restrict meshsize on edges - for(const auto & edge : edges) + for(auto i : Range(edges)) { + multithread.percent = 100. * i/edges.Size(); + const auto & edge = edges[i]; auto length = edge->GetLength(); // skip very short edges if(length < eps) @@ -127,9 +131,15 @@ namespace netgen } } + multithread.task = "Analyse Faces"; // restrict meshsize on faces - for(const auto& face : faces) - face->RestrictH(mesh, mparam); + for(auto i : Range(faces)) + { + multithread.percent = 100. * i/faces.Size(); + const auto& face = faces[i]; + face->RestrictH(mesh, mparam); + } + multithread.task = savetask; } mesh.LoadLocalMeshSize(mparam.meshsizefilename); } From 6c012675aa84fdbfe271441c78da2103da13ccb6 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Thu, 31 Oct 2019 17:08:29 +0100 Subject: [PATCH 4/5] project point without geominfo returns new geominfo --- libsrc/csg/csgeom.cpp | 3 ++- libsrc/csg/csgeom.hpp | 2 +- libsrc/meshing/basegeom.cpp | 5 ++++- libsrc/meshing/basegeom.hpp | 6 +++--- libsrc/meshing/meshing2.cpp | 2 +- libsrc/occ/occgeom.cpp | 9 ++++++--- libsrc/occ/occgeom.hpp | 2 +- libsrc/stlgeom/stlgeom.cpp | 2 +- libsrc/stlgeom/stlgeom.hpp | 2 +- 9 files changed, 20 insertions(+), 13 deletions(-) diff --git a/libsrc/csg/csgeom.cpp b/libsrc/csg/csgeom.cpp index f8145009..ad3912d4 100644 --- a/libsrc/csg/csgeom.cpp +++ b/libsrc/csg/csgeom.cpp @@ -72,11 +72,12 @@ namespace netgen Clean(); } - void CSGeometry :: ProjectPoint(int surfind, Point<3> & p) const + PointGeomInfo CSGeometry :: ProjectPoint(int surfind, Point<3> & p) const { Point<3> hp = p; GetSurface(surfind)->Project (hp); p = hp; + return PointGeomInfo(); } bool CSGeometry :: ProjectPointGI(int surfind, Point<3> & p, PointGeomInfo & gi) const diff --git a/libsrc/csg/csgeom.hpp b/libsrc/csg/csgeom.hpp index 07844c91..32473ed0 100644 --- a/libsrc/csg/csgeom.hpp +++ b/libsrc/csg/csgeom.hpp @@ -188,7 +188,7 @@ namespace netgen virtual void SaveToMeshFile (ostream & ost) const override; - void ProjectPoint(INDEX surfind, Point<3> & p) const override; + PointGeomInfo ProjectPoint(INDEX surfind, Point<3> & p) const override; bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override; void ProjectPointEdge(INDEX surfind, INDEX surfind2, Point<3> & p) const override; Vec<3> GetNormal(int surfind, const Point<3> & p) const override; diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp index d9542149..ec3930fa 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -95,6 +95,9 @@ namespace netgen mesh.SetLocalH(bounding_box.PMin(), bounding_box.PMax(), mparam.grading); + // only set meshsize for edges longer than this + double mincurvelength = 1e-3 * bounding_box.Diam(); + if(mparam.uselocalh) { double eps = 1e-12 * bounding_box.Diam(); @@ -108,7 +111,7 @@ namespace netgen const auto & edge = edges[i]; auto length = edge->GetLength(); // skip very short edges - if(length < eps) + if(length < mincurvelength) continue; static constexpr int npts = 20; // restrict mesh size based on edge length diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index b9e7bb16..53808739 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -50,7 +50,7 @@ namespace netgen virtual size_t GetNBoundaries() const = 0; virtual Array> GetBoundary(size_t index) const = 0; virtual string GetName() const { return "default"; } - virtual void Project(Point<3>& p) const = 0; + virtual PointGeomInfo Project(Point<3>& p) const = 0; // 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; @@ -124,9 +124,9 @@ namespace netgen virtual void FinalizeMesh(Mesh& mesh) const {} - virtual void ProjectPoint (int surfind, Point<3> & p) const + virtual PointGeomInfo ProjectPoint (int surfind, Point<3> & p) const { - faces[surfind-1]->Project(p); + return faces[surfind-1]->Project(p); } virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p) const diff --git a/libsrc/meshing/meshing2.cpp b/libsrc/meshing/meshing2.cpp index 070a31ef..e275bb22 100644 --- a/libsrc/meshing/meshing2.cpp +++ b/libsrc/meshing/meshing2.cpp @@ -176,7 +176,7 @@ namespace netgen { locpoint = p1 + (h*plainpoint(0)) * ex + (h* plainpoint(1)) * ey; if (!geo.ProjectPointGI(gi.trignum, locpoint, gi)) - geo.ProjectPoint(gi.trignum, locpoint); + gi = geo.ProjectPoint(gi.trignum, locpoint); return 0; } diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index 96434678..6b54ed66 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -1034,7 +1034,7 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a SetCenter(); } - void OCCGeometry :: ProjectPoint(int surfi, Point<3> & p) const + PointGeomInfo OCCGeometry :: ProjectPoint(int surfi, Point<3> & p) const { static int cnt = 0; if (++cnt % 1000 == 0) cout << "Project cnt = " << cnt << endl; @@ -1048,9 +1048,12 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a suval.Coord( u, v); pnt = thesurf->Value( u, v ); - + PointGeomInfo gi; + gi.trignum = surfi; + gi.u = u; + gi.v = v; p = Point<3> (pnt.X(), pnt.Y(), pnt.Z()); - + return gi; } bool OCCGeometry :: ProjectPointGI(int surfind, Point<3>& p, PointGeomInfo& gi) const diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index 09561b29..6dd2aab7 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -280,7 +280,7 @@ namespace netgen void DoArchive(Archive& ar) override; - void ProjectPoint(int surfind, Point<3> & p) const override; + PointGeomInfo ProjectPoint(int surfind, Point<3> & p) const override; void ProjectPointEdge (int surfind, int surfind2, Point<3> & p) const override; bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override; Vec<3> GetNormal(int surfind, const Point<3> & p) const override; diff --git a/libsrc/stlgeom/stlgeom.cpp b/libsrc/stlgeom/stlgeom.cpp index 7d254162..8adcaaaa 100644 --- a/libsrc/stlgeom/stlgeom.cpp +++ b/libsrc/stlgeom/stlgeom.cpp @@ -144,7 +144,7 @@ bool STLGeometry :: ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & g return true; } -void STLGeometry :: ProjectPoint (INDEX surfind, Point<3> & p) const +PointGeomInfo STLGeometry :: ProjectPoint (INDEX surfind, Point<3> & p) const { throw Exception("ProjectPoint without PointGeomInfo not implemented"); } diff --git a/libsrc/stlgeom/stlgeom.hpp b/libsrc/stlgeom/stlgeom.hpp index e29e4824..f542d11b 100644 --- a/libsrc/stlgeom/stlgeom.hpp +++ b/libsrc/stlgeom/stlgeom.hpp @@ -193,7 +193,7 @@ namespace netgen virtual void Save (string filename) const override; bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const override; - void ProjectPoint(INDEX surfind, Point<3> & p) const override; + PointGeomInfo ProjectPoint(INDEX surfind, Point<3> & p) const override; bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override; Vec<3> GetNormal(int surfind, const Point<3> & p) const override; Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const override; From 1b1c4700ad9ad231cdaae526101895e80257c185 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Mon, 4 Nov 2019 11:27:01 +0100 Subject: [PATCH 5/5] geo GetNormal and ProjectPointEdge with geoinfo pointer --- libsrc/csg/csgeom.cpp | 5 +++-- libsrc/csg/csgeom.hpp | 6 ++++-- libsrc/geom2d/geometry2d.cpp | 2 +- libsrc/geom2d/geometry2d.hpp | 2 +- libsrc/meshing/basegeom.hpp | 29 ++++++----------------------- libsrc/meshing/curvedelems.cpp | 6 +++--- libsrc/meshing/improve2.cpp | 6 +++--- libsrc/meshing/improve2gen.cpp | 2 +- libsrc/meshing/meshing2.cpp | 6 +++--- libsrc/meshing/smoothing2.cpp | 10 +++++----- libsrc/occ/occgeom.cpp | 31 +++++++++++++++---------------- libsrc/occ/occgeom.hpp | 6 +++--- libsrc/stlgeom/stlgeom.cpp | 11 ++++------- libsrc/stlgeom/stlgeom.hpp | 3 +-- 14 files changed, 53 insertions(+), 72 deletions(-) diff --git a/libsrc/csg/csgeom.cpp b/libsrc/csg/csgeom.cpp index ad3912d4..fa4c78f7 100644 --- a/libsrc/csg/csgeom.cpp +++ b/libsrc/csg/csgeom.cpp @@ -87,7 +87,7 @@ namespace netgen } void CSGeometry :: ProjectPointEdge(int surfind, INDEX surfind2, - Point<3> & p) const + Point<3> & p, EdgePointGeomInfo* /*unused*/) const { Point<3> hp = p; ProjectToEdge (GetSurface(surfind), @@ -96,7 +96,8 @@ namespace netgen } - Vec<3> CSGeometry :: GetNormal(int surfind, const Point<3> & p) const + Vec<3> CSGeometry :: GetNormal(int surfind, const Point<3> & p, + const PointGeomInfo* /*unused*/) const { Vec<3> hn; GetSurface(surfind)->CalcGradient(p, hn); diff --git a/libsrc/csg/csgeom.hpp b/libsrc/csg/csgeom.hpp index 32473ed0..c419c9a4 100644 --- a/libsrc/csg/csgeom.hpp +++ b/libsrc/csg/csgeom.hpp @@ -190,8 +190,10 @@ namespace netgen PointGeomInfo ProjectPoint(INDEX surfind, Point<3> & p) const override; bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override; - void ProjectPointEdge(INDEX surfind, INDEX surfind2, Point<3> & p) const override; - Vec<3> GetNormal(int surfind, const Point<3> & p) const override; + void ProjectPointEdge(INDEX surfind, INDEX surfind2, Point<3> & p, + EdgePointGeomInfo* gi = nullptr) const override; + Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi = nullptr) const override; + void PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint, int surfi, const PointGeomInfo & gi1, diff --git a/libsrc/geom2d/geometry2d.cpp b/libsrc/geom2d/geometry2d.cpp index 43282cff..11ea2df9 100644 --- a/libsrc/geom2d/geometry2d.cpp +++ b/libsrc/geom2d/geometry2d.cpp @@ -86,7 +86,7 @@ namespace netgen } Vec<3> SplineGeometry2d :: GetNormal(int surfi1, const Point<3> & p, - const PointGeomInfo & gi) const + const PointGeomInfo* gi) const { return Vec<3> (0,0,1); } diff --git a/libsrc/geom2d/geometry2d.hpp b/libsrc/geom2d/geometry2d.hpp index 7dfcbe9f..ff4c459c 100644 --- a/libsrc/geom2d/geometry2d.hpp +++ b/libsrc/geom2d/geometry2d.hpp @@ -182,7 +182,7 @@ namespace netgen Vec<3> GetTangent (const Point<3> & p, int surfi1, int surfi2, const EdgePointGeomInfo & ap1) const override; Vec<3> GetNormal(int surfi1, const Point<3> & p, - const PointGeomInfo & gi) const override; + const PointGeomInfo* gi) const override; const SplineSegExt & GetSpline (const int i) const { diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index 53808739..05dd1a90 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -32,7 +32,7 @@ namespace netgen virtual double CalcStep(double t, double sag) const = 0; virtual bool OrientedLikeGlobal() const = 0; virtual size_t GetHash() const = 0; - virtual void ProjectPoint(Point<3>& p, EdgePointGeomInfo& gi) const = 0; + virtual void ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const = 0; virtual void PointBetween(const Point<3>& p1, const Point<3>& p2, double secpoint, @@ -65,14 +65,8 @@ namespace netgen virtual double GetCurvature(const PointGeomInfo& gi) const = 0; virtual void RestrictH(Mesh& mesh, const MeshingParameters& mparam) const = 0; - virtual Vec<3> GetNormal(const Point<3>& p) const - { - return {0.,0.,1.}; - } - virtual Vec<3> GetNormal(const Point<3>& p, const PointGeomInfo& gi) const - { - return GetNormal(p); - } + virtual Vec<3> GetNormal(const Point<3>& p, const PointGeomInfo* gi = nullptr) const = 0; + virtual void PointBetween(const Point<3>& p1, const Point<3>& p2, double secpoint, @@ -129,13 +123,9 @@ namespace netgen return faces[surfind-1]->Project(p); } - virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p) const - { - throw Exception("In ProjectPointEdge of basegeometry"); - } - virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, EdgePointGeomInfo& gi) const + virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, EdgePointGeomInfo* gi = nullptr) const { - edges[gi.edgenr]->ProjectPoint(p, gi); + edges[gi->edgenr]->ProjectPoint(p, gi); } virtual bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const @@ -147,15 +137,8 @@ namespace netgen return faces[surfind-1]->ProjectPointGI(p, gi); } - virtual Vec<3> GetNormal(int surfind, const Point<3> & p) const - { return faces[surfind-1]->GetNormal(p); } - virtual Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const + virtual Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi = nullptr) const { return faces[surfind-1]->GetNormal(p, gi); } - [[deprecated]] - void GetNormal(int surfind, const Point<3> & p, Vec<3> & n) const - { - n = GetNormal(surfind, p); - } virtual void PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint, diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index 43d8ae35..71fb7d83 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -839,8 +839,8 @@ namespace netgen { Point<3> pm = Center (p1, p2); - Vec<3> n1 = geo.GetNormal (surfnr[e], p1, gi0[e]); - Vec<3> n2 = geo.GetNormal (surfnr[e], p2, gi1[e]); + Vec<3> n1 = geo.GetNormal (surfnr[e], p1, &gi0[e]); + Vec<3> n2 = geo.GetNormal (surfnr[e], p2, &gi1[e]); // p3 = pm + alpha1 n1 + alpha2 n2 @@ -1084,7 +1084,7 @@ namespace netgen v05 /= 1 + (w-1) * 0.5; Point<3> p05 (v05), pp05(v05); geo.ProjectPointEdge(edge_surfnr1[edgenr], edge_surfnr2[edgenr], pp05, - edge_gi0[edgenr]); + &edge_gi0[edgenr]); double d = Dist (pp05, p05); if (d < dold) diff --git a/libsrc/meshing/improve2.cpp b/libsrc/meshing/improve2.cpp index 524ba239..fe40259f 100644 --- a/libsrc/meshing/improve2.cpp +++ b/libsrc/meshing/improve2.cpp @@ -85,11 +85,11 @@ namespace netgen nv1.Normalize(); nv2.Normalize(); - auto nvp3 = geo.GetNormal (surfnr, mesh.Point(pi3), gi3); + auto nvp3 = geo.GetNormal (surfnr, mesh.Point(pi3), &gi3); nvp3.Normalize(); - auto nvp4 = geo.GetNormal (surfnr, mesh.Point(pi4), gi4); + auto nvp4 = geo.GetNormal (surfnr, mesh.Point(pi4), &gi4); nvp4.Normalize(); @@ -648,7 +648,7 @@ namespace netgen { const int faceindex = hel.GetIndex(); const int surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr(); - normals[pi] = geo.GetNormal (surfnr, mesh[pi], hel.GeomInfoPi(k+1)); + normals[pi] = geo.GetNormal (surfnr, mesh[pi], &hel.GeomInfoPi(k+1)); break; } } diff --git a/libsrc/meshing/improve2gen.cpp b/libsrc/meshing/improve2gen.cpp index 1d4ebbc9..ddbd397e 100644 --- a/libsrc/meshing/improve2gen.cpp +++ b/libsrc/meshing/improve2gen.cpp @@ -397,7 +397,7 @@ namespace netgen // calc metric badness double bad1 = 0, bad2 = 0; // SelectSurfaceOfPoint (mesh.Point(pmap.Get(1)), pgi.Get(1)); - auto n = geo.GetNormal(surfnr, mesh.Point(pmap.Get(1)), pgi.Elem(1)); + auto n = geo.GetNormal(surfnr, mesh.Point(pmap.Get(1)), &pgi.Elem(1)); for (int j = 0; j < rule.oldels.Size(); j++) bad1 += mesh[elmap[j]].CalcJacobianBadness (mesh.Points(), n); diff --git a/libsrc/meshing/meshing2.cpp b/libsrc/meshing/meshing2.cpp index e275bb22..bda90593 100644 --- a/libsrc/meshing/meshing2.cpp +++ b/libsrc/meshing/meshing2.cpp @@ -142,8 +142,8 @@ namespace netgen { p1 = ap1; p2 = ap2; - auto n1 = geo.GetNormal(gi1->trignum, p1, *gi1); - auto n2 = geo.GetNormal(gi2->trignum, p2, *gi2); + auto n1 = geo.GetNormal(gi1->trignum, p1, gi1); + auto n2 = geo.GetNormal(gi2->trignum, p2, gi2); ez = 0.5 * (n1+n2); ez.Normalize(); @@ -158,7 +158,7 @@ namespace netgen Point<2> & plainpoint, double h, int & zone) { auto& gi = geominfo.GetPGI(1); - auto n = geo.GetNormal(gi.trignum, locpoint, gi); + auto n = geo.GetNormal(gi.trignum, locpoint, &gi); auto p1p = locpoint - p1; plainpoint(0) = (p1p * ex) / h; plainpoint(1) = (p1p * ey) / h; diff --git a/libsrc/meshing/smoothing2.cpp b/libsrc/meshing/smoothing2.cpp index b4d5ab9a..11f21cb8 100644 --- a/libsrc/meshing/smoothing2.cpp +++ b/libsrc/meshing/smoothing2.cpp @@ -218,7 +218,7 @@ namespace netgen { double badness = 0; - auto n = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1); + auto n = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1); Point<3> pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2; for (int j = 0; j < ld.locelements.Size(); j++) @@ -359,7 +359,7 @@ namespace netgen vgrad = 0; double badness = 0; - auto n = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1); + auto n = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1); pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2; // meshthis -> ProjectPoint (surfi, pp1); @@ -414,7 +414,7 @@ namespace netgen vgrad = 0; double badness = 0; - auto n = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1); + auto n = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1); pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2; for (int j = 0; j < ld.locelements.Size(); j++) @@ -577,7 +577,7 @@ namespace netgen vgrad = 0; badness = 0; - auto n = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1); + // auto n = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1); pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2; @@ -933,7 +933,7 @@ namespace netgen } - ld.normal = geo.GetNormal(ld.surfi, ld.sp1, ld.gi1); + ld.normal = geo.GetNormal(ld.surfi, ld.sp1, &ld.gi1); ld.t1 = ld.normal.GetNormal (); ld.t2 = Cross (ld.normal, ld.t1); diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index 6b54ed66..6a441bf5 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -1072,7 +1072,7 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a } void OCCGeometry :: ProjectPointEdge(int surfind, INDEX surfind2, - Point<3> & p) const + Point<3> & p, EdgePointGeomInfo* gi) const { TopExp_Explorer exp0, exp1; bool done = false; @@ -1154,26 +1154,25 @@ void STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * a return true; } - Vec<3> OCCGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & geominfo) const + Vec<3> OCCGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* geominfo) const { - gp_Pnt pnt; - gp_Vec du, dv; + if(geominfo) + { + gp_Pnt pnt; + gp_Vec du, dv; - Handle(Geom_Surface) occface; - occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind))); + Handle(Geom_Surface) occface; + occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind))); - occface->D1(geominfo.u,geominfo.v,pnt,du,dv); + occface->D1(geominfo->u,geominfo->v,pnt,du,dv); - auto n = Cross (Vec<3>(du.X(), du.Y(), du.Z()), - Vec<3>(dv.X(), dv.Y(), dv.Z())); - n.Normalize(); + auto n = Cross (Vec<3>(du.X(), du.Y(), du.Z()), + Vec<3>(dv.X(), dv.Y(), dv.Z())); + n.Normalize(); - if (fmap(surfind).Orientation() == TopAbs_REVERSED) n *= -1; - return n; - } - - Vec<3> OCCGeometry :: GetNormal(int surfind, const Point<3> & p) const - { + if (fmap(surfind).Orientation() == TopAbs_REVERSED) n *= -1; + return n; + } Standard_Real u,v; gp_Pnt pnt(p(0), p(1), p(2)); diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index 6dd2aab7..94f3fd27 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -281,10 +281,10 @@ namespace netgen void DoArchive(Archive& ar) override; PointGeomInfo ProjectPoint(int surfind, Point<3> & p) const override; - void ProjectPointEdge (int surfind, int surfind2, Point<3> & p) const override; + void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, + EdgePointGeomInfo* gi = nullptr) const override; bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override; - Vec<3> GetNormal(int surfind, const Point<3> & p) const override; - Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const override; + Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const override; bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const override; void PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint, diff --git a/libsrc/stlgeom/stlgeom.cpp b/libsrc/stlgeom/stlgeom.cpp index 8adcaaaa..e65c0090 100644 --- a/libsrc/stlgeom/stlgeom.cpp +++ b/libsrc/stlgeom/stlgeom.cpp @@ -100,14 +100,11 @@ int STLGeometry :: GenerateMesh (shared_ptr & mesh, MeshingParameters & mp return STLMeshingDummy (this, mesh, mparam, stlpar); } -Vec<3> STLGeometry :: GetNormal(INDEX surfind, const Point<3> & p) const +Vec<3> STLGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const { - throw Exception("STLGeometry::GetNormal without PointGeomInfo called"); -} - -Vec<3> STLGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const -{ - return GetChart(GetChartNr(gi.trignum)).GetNormal(); + if(!gi) + throw Exception("STLGeometry::GetNormal without PointGeomInfo called"); + return GetChart(GetChartNr(gi->trignum)).GetNormal(); } bool STLGeometry :: CalcPointGeomInfo(int /*surfind*/, PointGeomInfo& gi, const Point<3> & p3) const diff --git a/libsrc/stlgeom/stlgeom.hpp b/libsrc/stlgeom/stlgeom.hpp index f542d11b..85212d55 100644 --- a/libsrc/stlgeom/stlgeom.hpp +++ b/libsrc/stlgeom/stlgeom.hpp @@ -195,8 +195,7 @@ namespace netgen bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const override; PointGeomInfo ProjectPoint(INDEX surfind, Point<3> & p) const override; bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override; - Vec<3> GetNormal(int surfind, const Point<3> & p) const override; - Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const override; + Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi = nullptr) const override; void PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint, int surfi, const PointGeomInfo & gi1,