diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp index d7489a81..7f160bac 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -28,6 +28,157 @@ namespace netgen void NetgenGeometry :: FindEdges(Mesh& mesh, const MeshingParameters& mparam) const { + static Timer t1("MeshEdges"); RegionTimer regt(t1); + + // create face descriptors and set bc names + mesh.SetNBCNames(faces.Size()); + for(auto i : Range(faces.Size())) + { + mesh.SetBCName(i, faces[i]->GetName()); + // todo find attached solids + FaceDescriptor fd(i+1, 1, 0, i+1); + fd.SetBCName(mesh.GetBCNamePtr(i)); + mesh.AddFaceDescriptor(fd); + } + + std::map vert2meshpt; + for(auto i : Range(vertices)) + { + const auto& vert = *vertices[i]; + MeshPoint mp(vert.GetPoint()); + vert2meshpt[vert.GetHash()] = mesh.AddPoint(mp); + } + + size_t segnr = 0; + for(auto facenr : Range(faces.Size())) + { + const auto& face = *faces[facenr]; + for(auto facebndnr : Range(face.GetNBoundaries())) + { + auto boundary = face.GetBoundary(facebndnr); + for(auto enr : Range(boundary)) + { + const auto& oriented_edge = *boundary[enr]; + auto edgenr = GetEdgeIndex(oriented_edge); + const auto& edge = edges[edgenr]; + PointIndex startp, endp; + // throws if points are not found + startp = vert2meshpt.at(edge->GetStartVertex().GetHash()); + endp = vert2meshpt.at(edge->GetEndVertex().GetHash()); + + // ignore collapsed edges + if(startp == endp && edge->GetLength() < 1e-10 * bounding_box.Diam()) + continue; + + + Array mps; + Array params; + // -------------------- DivideEdge ----------------- + static constexpr int divide_edge_sections = 1000; + double hvalue[divide_edge_sections+1]; + hvalue[0] = 0; + + Point<3> oldpnt; + auto pnt = edge->GetPoint(0.); + + // calc local h for edge + for(auto i : Range(divide_edge_sections)) + { + oldpnt = pnt; + pnt = edge->GetPoint(double(i+1)/divide_edge_sections); + hvalue[i+1] = hvalue[i] + 1./mesh.GetH(pnt) * (pnt-oldpnt).Length(); + } + int nsubedges = max2(1, int(floor(hvalue[divide_edge_sections]+0.5))); + mps.SetSize(nsubedges-1); + params.SetSize(nsubedges+1); + + int i = 1; + int i1 = 0; + do + { + if (hvalue[i1]/hvalue[divide_edge_sections]*nsubedges >= i) + { + params[i] = (i1/double(divide_edge_sections)); + pnt = edge->GetPoint(params[i]); + mps[i-1] = MeshPoint(pnt); + i++; + } + i1++; + if (i1 > divide_edge_sections) + { + nsubedges = i; + mps.SetSize(nsubedges-1); + params.SetSize(nsubedges+1); + cout << "divide edge: local h too small" << endl; + } + + } while(i < nsubedges); + + params[0] = 0.; + params[nsubedges] = 1.; + + if(params[nsubedges] <= params[nsubedges-1]) + { + cout << "CORRECTED" << endl; + mps.SetSize (nsubedges-2); + params.SetSize (nsubedges); + params[nsubedges] = 1.; + } + // ----------- Add Points to mesh and create segments ----- + + Array pnums(mps.Size() + 2); + pnums[0] = startp; + pnums[mps.Size()+1] = endp; + + double eps = bounding_box.Diam() * 1e-8; + + for(auto i : Range(mps)) + { + bool exists = false; + for(auto pi : Range(mesh.Points())) + { + if((mesh[pi] - mps[i]).Length() < eps) + { + exists = true; + pnums[i+1] = pi; + break; + } + } + if(!exists) + pnums[i+1] = mesh.AddPoint(mps[i]); + } + + for(auto i : Range(pnums.Size()-1)) + { + segnr++; + Segment seg; + seg[0] = pnums[i]; + seg[1] = pnums[i+1]; + seg.edgenr = segnr; + seg.epgeominfo[0].dist = params[i]; + seg.epgeominfo[1].dist = params[i+1]; + seg.epgeominfo[0].edgenr = edgenr; + seg.epgeominfo[1].edgenr = edgenr; + seg.si = facenr+1; + seg.surfnr1 = facenr+1; + + // TODO: implement functionality to transfer edge parameter t to face parameters u,v + for(auto j : Range(2)) + face.CalcEdgePointGI(*edge, params[i+j], + seg.epgeominfo[j]); + + if(!oriented_edge.OrientedLikeGlobal()) + { + swap (seg[0], seg[1]); + swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist); + swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u); + swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v); + } + mesh.AddSegment(seg); + } + } + } + } } void NetgenGeometry :: MeshSurface(Mesh& mesh, diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index 4efe2cbb..3835937c 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -12,10 +12,24 @@ struct Tcl_Interp; namespace netgen { + class GeometryVertex + { + public: + virtual ~GeometryVertex() {} + virtual Point<3> GetPoint() const = 0; + virtual size_t GetHash() const = 0; + }; + class GeometryEdge { public: virtual ~GeometryEdge() {} + virtual const GeometryVertex& GetStartVertex() const = 0; + virtual const GeometryVertex& GetEndVertex() const = 0; + virtual double GetLength() const = 0; + virtual Point<3> GetPoint(double t) const = 0; + virtual bool OrientedLikeGlobal() const = 0; + virtual size_t GetHash() const = 0; }; class GeometryFace @@ -23,10 +37,14 @@ namespace netgen public: virtual ~GeometryFace() {} virtual size_t GetNBoundaries() const = 0; - virtual Array GetBoundary(size_t index) const = 0; + virtual Array> GetBoundary(size_t index) const = 0; + virtual string GetName() const { return "default"; } // 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 void CalcEdgePointGI(const GeometryEdge& edge, + double t, + EdgePointGeomInfo& egi) const = 0; virtual Box<3> GetBoundingBox() const = 0; }; @@ -34,6 +52,8 @@ namespace netgen { unique_ptr ref; protected: + Array> vertices; + Array> edges; Array> faces; Box<3> bounding_box; public: @@ -112,6 +132,14 @@ namespace netgen int surfi2, const EdgePointGeomInfo & egi) const { throw Exception("Call GetTangent of " + Demangle(typeid(*this).name())); } + + virtual size_t GetEdgeIndex(const GeometryEdge& edge) const + { + for(auto i : Range(edges)) + if(edge.GetHash() == edges[i]->GetHash()) + return i; + throw Exception("Couldn't find edge index"); + } virtual void Save (string filename) const; virtual void SaveToMeshFile (ostream & /* ost */) const { ; } };