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 {