From 16b88e8e67a0e52d4de3e17b43343d0a09a3c790 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Sun, 28 Nov 2021 15:14:41 +0000 Subject: [PATCH] OCC Mesher Cleanup --- cmake/SuperBuild.cmake | 6 +- libsrc/gprim/adtree.hpp | 2 + libsrc/meshing/basegeom.cpp | 699 ++++++++++++++------ libsrc/meshing/basegeom.hpp | 93 ++- libsrc/meshing/meshtype.cpp | 1 + libsrc/meshing/meshtype.hpp | 15 +- libsrc/meshing/python_mesh.cpp | 13 +- libsrc/occ/CMakeLists.txt | 12 +- libsrc/occ/occ_edge.cpp | 84 +++ libsrc/occ/occ_edge.hpp | 44 ++ libsrc/occ/occ_face.cpp | 254 ++++++++ libsrc/occ/occ_face.hpp | 49 ++ libsrc/occ/occ_solid.hpp | 27 + libsrc/occ/occ_utils.cpp | 24 + libsrc/occ/occ_utils.hpp | 267 ++++++++ libsrc/occ/occ_vertex.cpp | 25 + libsrc/occ/occ_vertex.hpp | 28 + libsrc/occ/occgenmesh.cpp | 1044 ++++++------------------------ libsrc/occ/occgeom.cpp | 694 ++++++++++---------- libsrc/occ/occgeom.hpp | 254 +------- libsrc/occ/occmeshsurf.cpp | 6 +- libsrc/occ/occmeshsurf.hpp | 12 +- libsrc/occ/occpkg.cpp | 3 + libsrc/occ/python_occ.cpp | 61 +- libsrc/occ/python_occ.hpp | 68 -- libsrc/occ/python_occ_basic.cpp | 52 +- libsrc/occ/python_occ_shapes.cpp | 284 +++----- libsrc/occ/vsocc.cpp | 25 +- nglib/nglib.cpp | 6 +- python/occ.py | 4 + tests/dockerfile | 3 + tests/dockerfile_mpi | 23 +- 32 files changed, 2139 insertions(+), 2043 deletions(-) create mode 100644 libsrc/occ/occ_edge.cpp create mode 100644 libsrc/occ/occ_edge.hpp create mode 100644 libsrc/occ/occ_face.cpp create mode 100644 libsrc/occ/occ_face.hpp create mode 100644 libsrc/occ/occ_solid.hpp create mode 100644 libsrc/occ/occ_utils.cpp create mode 100644 libsrc/occ/occ_utils.hpp create mode 100644 libsrc/occ/occ_vertex.cpp create mode 100644 libsrc/occ/occ_vertex.hpp delete mode 100644 libsrc/occ/python_occ.hpp diff --git a/cmake/SuperBuild.cmake b/cmake/SuperBuild.cmake index c5a1b01b..1301bde1 100644 --- a/cmake/SuperBuild.cmake +++ b/cmake/SuperBuild.cmake @@ -101,8 +101,8 @@ if(BUILD_OCC) ExternalProject_Add(project_occ DEPENDS project_freetype - URL https://github.com/Open-Cascade-SAS/OCCT/archive/refs/tags/V7_5_0.zip - URL_MD5 a24e6d3cf2d24bf9347d2d4aee9dd80a + URL https://github.com/Open-Cascade-SAS/OCCT/archive/refs/tags/V7_6_0.zip + URL_MD5 37519251c99cb3469ccfa82a9241d528 DOWNLOAD_DIR ${CMAKE_CURRENT_SOURCE_DIR}/external_dependencies ${SUBPROJECT_ARGS} CMAKE_ARGS @@ -116,7 +116,7 @@ if(BUILD_OCC) -DBUILD_MODULE_DataExchange:BOOL=ON -DBUILD_MODULE_ApplicationFramework:BOOL=OFF -DBUILD_MODULE_Draw:BOOL=OFF - -DUSE_FREETYPE=OFF + -DUSE_FREETYPE=ON ${SUBPROJECT_CMAKE_ARGS} UPDATE_COMMAND "" ) diff --git a/libsrc/gprim/adtree.hpp b/libsrc/gprim/adtree.hpp index f05f7170..a1768493 100644 --- a/libsrc/gprim/adtree.hpp +++ b/libsrc/gprim/adtree.hpp @@ -817,6 +817,8 @@ public: : BoxTree(box.PMin(), box.PMax()) { } + double GetTolerance() { return tol; } + size_t GetNLeaves() { return n_leaves; diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp index 62c0b981..870c1dd3 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -4,13 +4,67 @@ namespace netgen { + struct PointTree + { + BoxTree<3> tree; - DLL_HEADER GeometryRegisterArray geometryregister; + PointTree( Box<3> bb ) : tree(bb) {} + + void Insert(Point<3> p, PointIndex n) + { + tree.Insert(p, p, n); + } + + PointIndex Find(Point<3> p) const + { + ArrayMem points; + tree.GetIntersecting(p, p, points); + if(points.Size()==0) + throw Exception("cannot find mapped point"); + return points[0]; + } + + double GetTolerance() { return tree.GetTolerance(); } + }; + + DLL_HEADER GeometryRegisterArray geometryregister; //DLL_HEADER NgArray geometryregister; GeometryRegister :: ~GeometryRegister() { ; } + bool GeometryShape :: IsMappedShape( const GeometryShape & other_, const Transformation<3> & trafo, double tol ) const + { + throw Exception("GeometryShape::IsMappedShape not implemented for class " + Demangle(typeid(this).name())); + } + + bool GeometryVertex :: IsMappedShape( const GeometryShape & other_, const Transformation<3> & trafo, double tol ) const + { + const auto other_ptr = dynamic_cast(&other_); + if(!other_ptr) + return false; + + return Dist(trafo(GetPoint()), other_ptr->GetPoint()) < tol; + } + + bool GeometryEdge :: IsMappedShape( const GeometryShape & other_, const Transformation<3> & trafo, double tol ) const + { + const auto other_ptr = dynamic_cast(&other_); + if(!other_ptr) + return false; + auto & e = *other_ptr; + + if(tol < Dist(GetCenter(), e.GetCenter())) + return false; + + auto &v0 = GetStartVertex(); + auto &v1 = GetEndVertex(); + auto &w0 = e.GetStartVertex(); + auto &w1 = e.GetEndVertex(); + return( (v0.IsMappedShape(w0, trafo, tol) && v1.IsMappedShape(w1, trafo, tol)) || + (v0.IsMappedShape(w1, trafo, tol) && v1.IsMappedShape(w0, trafo, tol)) ); + } + void GeometryFace :: RestrictHTrig(Mesh& mesh, const PointGeomInfo& gi0, const PointGeomInfo& gi1, @@ -103,6 +157,64 @@ namespace netgen } }; + void NetgenGeometry :: ProcessIdentifications() + { + auto mirror_identifications = [&] ( auto & shapes ) + { + for(auto i : Range(shapes)) + { + auto &s = shapes[i]; + s->nr = i; + for(auto & ident : s->identifications) + if(s.get() == ident.from) + ident.to->identifications.Append(ident); + } + }; + + mirror_identifications(vertices); + mirror_identifications(edges); + mirror_identifications(faces); + + // todo: propagate identifications faces -> edges -> vertices + + auto find_primary = [&] (auto & shapes) + { + for(auto &s : shapes) + { + s->primary = s.get(); + s->primary_to_me = Transformation<3>{ Vec<3> {0,0,0} }; // init with identity + } + + bool changed = true; + + while(changed) { + changed = false; + for(auto &s : shapes) + { + auto current = s->primary; + for(auto & ident : current->identifications) + { + bool need_inverse = ident.from == s.get(); + auto other = need_inverse ? ident.to : ident.from; + if(other->nr < current->nr) + { + auto trafo = ident.trafo; + if(need_inverse) + trafo = trafo.CalcInverse(); + s->primary = other; + s->primary_to_me.Combine(trafo, s->primary_to_me); + changed = true; + } + } + } + } + }; + + find_primary(vertices); + find_primary(edges); + find_primary(faces); + } + void NetgenGeometry :: Analyse(Mesh& mesh, const MeshingParameters& mparam) const { @@ -241,165 +353,257 @@ namespace netgen mesh.LoadLocalMeshSize(mparam.meshsizefilename); } + void DivideEdge(GeometryEdge * edge, const MeshingParameters & mparam, const Mesh & mesh, Array> & points, Array & params) + { + static Timer tdivedgesections("Divide edge sections"); + static Timer tdivide("Divide Edges"); + RegionTimer rt(tdivide); + // -------------------- DivideEdge ----------------- + static constexpr size_t divide_edge_sections = 1000; + double hvalue[divide_edge_sections+1]; + hvalue[0] = 0; + + Point<3> old_pt = edge->GetPoint(0.); + // calc local h for edge + tdivedgesections.Start(); + 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(); + points.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] = (double(i1)/divide_edge_sections); + points[i-1] = MeshPoint(edge->GetPoint(params[i])); + i++; + } + i1++; + if (i1 > divide_edge_sections) + { + nsubedges = i; + points.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; + points.SetSize (nsubedges-2); + params.SetSize (nsubedges); + params[nsubedges-1] = 1.; + } + } + void NetgenGeometry :: FindEdges(Mesh& mesh, const MeshingParameters& mparam) const { 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()); - 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); - } + PointTree tree( bounding_box ); + + auto & identifications = mesh.GetIdentifications(); std::map vert2meshpt; - for(auto i : Range(vertices)) + for(auto & vert : vertices) { - const auto& vert = *vertices[i]; - MeshPoint mp(vert.GetPoint()); - vert2meshpt[vert.GetHash()] = mesh.AddPoint(mp); + auto pi = mesh.AddPoint(vert->GetPoint()); + tree.Insert(mesh[pi], pi); + vert2meshpt[vert->GetHash()] = pi; + mesh[pi].Singularity(vert->properties.hpref); + + if(vert->properties.name) + { + Element0d el(pi, pi); + el.name = vert->properties.GetName(); + mesh.SetCD3Name(pi, el.name); + mesh.pointelements.Append (el); + } } + for(auto & vert : vertices) + for(auto & ident : vert->identifications) + identifications.Add(vert2meshpt[ident.from->GetHash()], + vert2meshpt[ident.to->GetHash()], + ident.name, + ident.type); + size_t segnr = 0; - for(auto facenr : Range(faces.Size())) - { - const auto& face = *faces[facenr]; - for(auto facebndnr : Range(face.GetNBoundaries())) + auto nedges = edges.Size(); + Array> all_pnums(nedges); + Array> all_params(nedges); + + for (auto edgenr : Range(edges)) + { + auto edge = edges[edgenr].get(); + 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; + + // ----------- Add Points to mesh and create segments ----- + auto & pnums = all_pnums[edgenr]; + auto & params = all_params[edgenr]; + Array> edge_points; + Array edge_params; + + if(edge->primary == edge) + { + DivideEdge(edge, mparam, mesh, edge_points, edge_params); + } + else + { + auto nr_primary = edge->primary->nr; + auto & pnums_primary = all_pnums[nr_primary]; + auto & params_primary = all_params[nr_primary]; + auto trafo = edge->primary_to_me; + + auto np = pnums_primary.Size(); + edge_points.SetSize(np-2); + edge_params.SetSize(np-2); + for(auto i : Range(np-2)) + { + edge_points[i] = trafo(mesh[pnums_primary[i+1]]); + EdgePointGeomInfo gi; + edge->ProjectPoint(edge_points[i], &gi); + edge_params[i] = gi.dist; + } + + // reverse entries if we have decreasing parameters + if(edge_params.Size()>2 && edge_params[0] > edge_params.Last()) + for(auto i : Range((np-2)/2)) + { + swap(edge_points[i], edge_points[np-3-i]); + swap(edge_params[i], edge_params[np-3-i]); + } + } + + pnums.SetSize(edge_points.Size() + 2); + pnums[0] = startp; + pnums.Last() = endp; + + params.SetSize(edge_points.Size()+2); + params[0] = 0.; + params.Last() = 1.; + + for(auto i : Range(edge_points)) + { + auto pi = mesh.AddPoint(edge_points[i]); + tree.Insert(mesh[pi], pi); + pnums[i+1] = pi; + params[i+1] = edge_params[i]; + } + + for(auto i : Range(pnums.Size()-1)) + { + segnr++; + Segment seg; + seg[0] = pnums[i]; + seg[1] = pnums[i+1]; + seg.edgenr = edgenr+1; + seg.si = edgenr+1; + seg.epgeominfo[0].dist = params[i]; + seg.epgeominfo[1].dist = params[i+1]; + seg.epgeominfo[0].edgenr = edgenr; + seg.epgeominfo[1].edgenr = edgenr; + seg.singedge_left = edge->properties.hpref; + seg.singedge_right = edge->properties.hpref; + mesh.AddSegment(seg); + } + mesh.SetCD2Name(edgenr+1, edge->properties.GetName()); + } + + for (auto & edge : edges) + { + // identify points on edge + for(auto & ident : edge->identifications) + if(ident.from == edge.get()) { - auto boundary = face.GetBoundary(facebndnr); - for(auto enr : Range(boundary)) + auto & pnums = all_pnums[edge->nr]; + // start and end vertex are already identified + for(auto pi : pnums.Range(1, pnums.Size()-1)) + { + auto pi_other = tree.Find(ident.trafo(mesh[pi])); + identifications.Add(pi, pi_other, ident.name, ident.type); + } + } + } + mesh.CalcSurfacesOfNode(); + multithread.task = savetask; + } + + bool NetgenGeometry :: MeshFace(Mesh& mesh, const MeshingParameters& mparam, + int k, FlatArray glob2loc) const + { + multithread.percent = 100. * k/faces.Size(); + const auto& face = *faces[k]; + auto bb = face.GetBoundingBox(); + bb.Increase(bb.Diam()/10); + Meshing2 meshing(*this, mparam, bb); + glob2loc = 0; + int cntp = 0; + + auto segments = face.GetBoundary(mesh); + for(auto& seg : segments) + { + for(auto j : Range(2)) + { + auto pi = seg[j]; + if(glob2loc[pi] == 0) { - 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]; - 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 size_t divide_edge_sections = 1000; - tdivide.Start(); - double hvalue[divide_edge_sections+1]; - hvalue[0] = 0; - - Point<3> old_pt = edge->GetPoint(0.); - // calc local h for edge - tdivedgesections.Start(); - 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); - params.SetSize(nsubedges+1); - - int i = 1; - int i1 = 0; - do - { - if (hvalue[i1]/hvalue[divide_edge_sections]*nsubedges >= i) - { - params[i] = (double(i1)/divide_edge_sections); - mps[i-1] = MeshPoint(edge->GetPoint(params[i])); - 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] = 1.; - } - tdivide.Stop(); - // ----------- 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); - } + meshing.AddPoint(mesh[pi], pi); + cntp++; + glob2loc[pi] = cntp; } } } - mesh.CalcSurfacesOfNode(); - multithread.task = savetask; + for(auto & seg : segments) + { + PointGeomInfo gi0, gi1; + gi0.trignum = gi1.trignum = k+1; + gi0.u = seg.epgeominfo[0].u; + gi0.v = seg.epgeominfo[0].v; + gi1.u = seg.epgeominfo[1].u; + gi1.v = seg.epgeominfo[1].v; + meshing.AddBoundaryElement(glob2loc[seg[0]], + glob2loc[seg[1]], + gi0, gi1); + } + + // TODO Set max area 2* area of face + + auto noldsurfels = mesh.GetNSE(); + + + static Timer t("GenerateMesh"); RegionTimer reg(t); + MESHING2_RESULT res = meshing.GenerateMesh(mesh, mparam, mparam.maxh, k+1); + + for(auto i : Range(noldsurfels, mesh.GetNSE())) + { + mesh.SurfaceElements()[i].SetIndex(k+1); + } + return res != MESHING2_OK; } void NetgenGeometry :: MeshSurface(Mesh& mesh, @@ -409,65 +613,160 @@ namespace netgen const char* savetask = multithread.task; multithread.task = "Mesh Surface"; + size_t n_failed_faces = 0; 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); - Meshing2 meshing(*this, mparam, bb); - glob2loc = 0; - int cntp = 0; + { + auto & face = *faces[k]; + mesh.SetBCName(k, face.properties.GetName()); + // todo find attached solids + FaceDescriptor fd(k+1, face.domin+1, face.domout+1, k+1); + fd.SetBCName(mesh.GetBCNamePtr(k)); + mesh.AddFaceDescriptor(fd); + if(face.primary == &face) + { + if(MeshFace(mesh, mparam, k, glob2loc)) + n_failed_faces++; + } + } - for(auto& seg : mesh.LineSegments()) - { - if(seg.si == k+1) - { - for(auto j : Range(2)) - { - auto pi = seg[j]; - if(glob2loc[pi] == 0) - { - meshing.AddPoint(mesh[pi], pi); - cntp++; - glob2loc[pi] = cntp; - } - } - } - } - for(auto & seg : mesh.LineSegments()) - { - if(seg.si == k+1) - { - PointGeomInfo gi0, gi1; - gi0.trignum = gi1.trignum = k+1; - gi0.u = seg.epgeominfo[0].u; - gi0.v = seg.epgeominfo[0].v; - gi1.u = seg.epgeominfo[1].u; - gi1.v = seg.epgeominfo[1].v; - meshing.AddBoundaryElement(glob2loc[seg[0]], - glob2loc[seg[1]], - gi0, gi1); - } - } + if(n_failed_faces) + { + cout << "WARNING! NOT ALL FACES HAVE BEEN MESHED" << endl; + cout << "SURFACE MESHING ERROR OCCURRED IN " << n_failed_faces << " FACES:" << endl; + return; + } - // TODO Set max area 2* area of face + if (mparam.perfstepsend >= MESHCONST_OPTSURFACE) + { + mesh.CalcSurfacesOfNode(); + OptimizeSurface(mesh, mparam); + } - auto noldsurfels = mesh.GetNSE(); + bool have_identifications = false; + for(auto & face : faces) + if(face->primary != face.get()) + { + have_identifications = true; + MapSurfaceMesh(mesh, *face); + } + // identify points on faces + if(have_identifications) + { + mesh.CalcSurfacesOfNode(); + BitArray is_identified_face(faces.Size()); + is_identified_face = false; + for(auto & face : faces) + for(auto & ident : face->identifications) + { + is_identified_face.SetBit(ident.from->nr); + is_identified_face.SetBit(ident.to->nr); + } - static Timer t("GenerateMesh"); RegionTimer reg(t); - MESHING2_RESULT res = meshing.GenerateMesh(mesh, mparam, mparam.maxh, k+1); + PointTree tree( bounding_box ); + Array pi_to_face(mesh.GetNP()); + pi_to_face = -1; + Array si_of_face; + Array> pi_of_face(faces.Size()); + for(auto & face : faces) + if(is_identified_face[face->nr]) + { + mesh.GetSurfaceElementsOfFace(face->nr+1, si_of_face); + for(auto si : si_of_face) + for(auto pi : mesh[si].PNums()) + { + if(mesh[pi].Type() == SURFACEPOINT && pi_to_face[pi]==-1) + { + pi_to_face[pi] = face->nr; + tree.Insert(mesh[pi], pi); + pi_of_face[face->nr].Append(pi); + } + } + } - for(auto i : Range(noldsurfels, mesh.GetNSE())) - { - mesh.SurfaceElements()[i].SetIndex(k+1); - } - } + auto & mesh_ident = mesh.GetIdentifications(); + for(auto & face : faces) + for(auto & ident : face->identifications) + { + if(ident.from == face.get()) + for(auto pi : pi_of_face[face->nr]) + { + auto pi_other = tree.Find(ident.trafo(mesh[pi])); + mesh_ident.Add(pi, pi_other, ident.name, ident.type); + } + } + } + + mesh.CalcSurfacesOfNode(); multithread.task = savetask; } + void NetgenGeometry :: MapSurfaceMesh( Mesh & mesh, const GeometryFace & dst ) const + { + static Timer timer("MapSurfaceMesh"); + RegionTimer rt(timer); + + const auto & src = dynamic_cast(*dst.primary); + auto trafo = dst.primary_to_me; + + PrintMessage(2, "Map face ", src.nr+1, " -> ", dst.nr+1); + + // point map from src to dst + Array pmap(mesh.Points().Size()); + pmap = PointIndex::INVALID; + + // first map points on edges (mapped points alread in mesh, use search tree) + Array is_point_in_tree(mesh.Points().Size()); + is_point_in_tree = false; + PointTree tree( bounding_box ); + + for (Segment & seg : src.GetBoundary(mesh)) + for(auto i : Range(2)) + { + auto pi = seg[i]; + if(!is_point_in_tree[pi]) + { + tree.Insert(trafo(mesh[pi]), pi); + is_point_in_tree[pi] = true; + } + } + + for (Segment & seg : dst.GetBoundary(mesh)) + for(auto i : Range(2)) + { + auto pi = seg[i]; + if(pmap[pi].IsValid()) + continue; + + pmap[tree.Find(mesh[pi])] = pi; + } + + // now insert mapped surface elements + for(auto sei : mesh.SurfaceElements().Range()) + { + auto sel = mesh[sei]; + if(sel.GetIndex() != src.nr+1) + continue; + + auto sel_new = sel; + sel_new.SetIndex(dst.nr+1); + for(auto i : Range(sel.PNums())) + { + auto pi = sel[i]; + if(!pmap[pi].IsValid()) + { + pmap[pi] = mesh.AddPoint(trafo(mesh[pi]), 1, SURFACEPOINT); + } + sel_new[i] = pmap[pi]; + } + sel_new.Invert(); + for(auto i : Range(sel.PNums())) + dst.CalcPointGeomInfo(mesh[sel_new[i]], sel_new.GeomInfo()[i]); + mesh.AddSurfaceElement(sel_new); + } + } + void NetgenGeometry :: OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const { const auto savetask = multithread.task; @@ -504,6 +803,13 @@ namespace netgen mesh.Compress(); multithread.task = savetask; } + + void NetgenGeometry :: FinalizeMesh(Mesh& mesh) const + { + for (int i = 0; i < mesh.GetNDomains(); i++) + if (auto name = solids[i]->properties.name) + mesh.SetMaterial (i+1, *name); + } shared_ptr GeometryRegisterArray :: LoadFromMeshFile (istream & ist) const { @@ -567,18 +873,17 @@ namespace netgen if (mparam.perfstepsstart <= MESHCONST_MESHSURFACE) { MeshSurface(*mesh, mparam); - mesh->CalcSurfacesOfNode(); } - if (multithread.terminate || mparam.perfstepsend <= MESHCONST_MESHSURFACE) - return 0; - - if (mparam.perfstepsstart <= MESHCONST_OPTSURFACE) - OptimizeSurface(*mesh, mparam); - if (multithread.terminate || mparam.perfstepsend <= MESHCONST_OPTSURFACE) return 0; + if(dimension == 2) + { + FinalizeMesh(*mesh); + mesh->SetDimension(2); + return 0; + } if(mparam.perfstepsstart <= MESHCONST_MESHVOLUME) { diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index 8d367966..c813be83 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -12,26 +12,72 @@ struct Tcl_Interp; namespace netgen { - class GeometryVertex + struct ShapeProperties { - public: - virtual ~GeometryVertex() {} - virtual Point<3> GetPoint() const = 0; - virtual size_t GetHash() const = 0; + optional name; + optional> col; + double maxh = 1e99; + double hpref = 0; // number of hp refinement levels (will be multiplied by factor later) + void Merge(const ShapeProperties & prop2) + { + if (prop2.name) name = prop2.name; + if (prop2.col) col = prop2.col; + maxh = min2(maxh, prop2.maxh); + hpref = max2(hpref, prop2.hpref); + } + + string GetName() const { return name ? *name : "default"; } + Vec<4> GetColor() { return col ? *col : Vec<4>{0., 1., 0., 1.}; } + + void DoArchive(Archive& ar) + { + ar & name & col & maxh & hpref; + } }; - class GeometryEdge + class GeometryShape; + + struct ShapeIdentification + { + GeometryShape * from; + GeometryShape * to; + Transformation<3> trafo; + Identifications::ID_TYPE type; + string name = ""; + }; + + class DLL_HEADER GeometryShape + { + public: + int nr = -1; + ShapeProperties properties; + Array identifications; + GeometryShape * primary; + Transformation<3> primary_to_me; + + virtual ~GeometryShape() {} + virtual size_t GetHash() const = 0; + virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const; + }; + + + class DLL_HEADER GeometryVertex : public GeometryShape + { + public: + virtual Point<3> GetPoint() const = 0; + virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const override; + }; + + class DLL_HEADER GeometryEdge : public GeometryShape { public: - virtual ~GeometryEdge() {} virtual const GeometryVertex& GetStartVertex() const = 0; virtual const GeometryVertex& GetEndVertex() const = 0; virtual double GetLength() const = 0; + virtual Point<3> GetCenter() 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; virtual void ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const = 0; virtual void PointBetween(const Point<3>& p1, const Point<3>& p2, @@ -46,15 +92,17 @@ namespace netgen ProjectPoint(newp, &newgi); } virtual Vec<3> GetTangent(double t) const = 0; + virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const override; }; - class GeometryFace + class DLL_HEADER GeometryFace : public GeometryShape { public: - virtual ~GeometryFace() {} + int domin=-1, domout=-1; + + virtual Point<3> GetCenter() const = 0; virtual size_t GetNBoundaries() const = 0; - virtual Array> GetBoundary(size_t index) const = 0; - virtual string GetName() const { return "default"; } + virtual Array GetBoundary(const Mesh& mesh) 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. @@ -102,6 +150,9 @@ namespace netgen int depth = 0, double h = 0.) const; }; + class DLL_HEADER GeometrySolid : public GeometryShape + { }; + class DLL_HEADER NetgenGeometry { unique_ptr ref; @@ -109,15 +160,25 @@ namespace netgen Array> vertices; Array> edges; Array> faces; + Array> solids; Array, double>> restricted_h; Box<3> bounding_box; + int dimension = 3; public: + NetgenGeometry() { ref = make_unique(*this); } virtual ~NetgenGeometry () { ; } + size_t GetNVertices() const { return vertices.Size(); } + size_t GetNEdges() const { return edges.Size(); } + size_t GetNFaces() const { return faces.Size(); } + + const GeometryFace & GetFace(int i) const { return *faces[i]; } + + virtual int GenerateMesh (shared_ptr & mesh, MeshingParameters & mparam); void RestrictH(const Point<3>& pnt, double maxh) @@ -134,13 +195,17 @@ namespace netgen { throw NgException("DoArchive not implemented for " + Demangle(typeid(*this).name())); } virtual Mesh::GEOM_TYPE GetGeomType() const { return Mesh::NO_GEOM; } + virtual void ProcessIdentifications(); virtual void Analyse(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 bool MeshFace(Mesh& mesh, const MeshingParameters& mparam, + int nr, FlatArray glob2loc) const; + virtual void MapSurfaceMesh( Mesh & mesh, const GeometryFace & dst ) const; virtual void OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const; - virtual void FinalizeMesh(Mesh& mesh) const {} + virtual void FinalizeMesh(Mesh& mesh) const; virtual PointGeomInfo ProjectPoint (int surfind, Point<3> & p) const { diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index c9d4b4fa..8596c753 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -2687,6 +2687,7 @@ namespace netgen identifiedpoints_nr.Set (tripl, 1); if (identnr > maxidentnr) maxidentnr = identnr; + names.SetSize(maxidentnr); if (identnr+1 > idpoints_table.Size()) idpoints_table.ChangeSize (identnr+1); diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index f3c5e536..1c03401a 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -1497,6 +1497,7 @@ namespace netgen /// number of identifications (or, actually used identifications ?) int maxidentnr; + Array names; public: /// @@ -1511,7 +1512,12 @@ namespace netgen identification nr identnr */ DLL_HEADER void Add (PointIndex pi1, PointIndex pi2, int identnr); - + void Add (PointIndex pi1, PointIndex pi2, string name, ID_TYPE type) + { + auto nr = GetNr(name); + Add(pi1, pi2, nr); + SetType(nr, type); + } int Get (PointIndex pi1, PointIndex pi2) const; int GetSymmetric (PointIndex pi1, PointIndex pi2) const; @@ -1560,6 +1566,13 @@ namespace netgen /// int GetMaxNr () const { return maxidentnr; } + int GetNr(string name) + { + if(!names.Contains(name)) + names.Append(name); + return names.Pos(name)+1; + } + /// remove secondorder void SetMaxPointNr (int maxpnum); diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 40d63e78..1d22044b 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -117,6 +117,13 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) m.def("_PushStatus", [](string s) { PushStatus(MyStr(s)); }); m.def("_SetThreadPercentage", [](double percent) { SetThreadPercent(percent); }); + py::enum_(m,"IdentificationType") + .value("UNDEFINED", Identifications::UNDEFINED) + .value("PERIODIC", Identifications::PERIODIC) + .value("CLOSESURFACES", Identifications::CLOSESURFACES) + .value("CLOSEEDGES", Identifications::CLOSEEDGES) + ; + py::class_ (m, "MPI_Comm") #ifdef NG_MPI4PY .def(py::init([] (mpi4py_comm comm) @@ -944,19 +951,19 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def ("GetCD3Name", &Mesh::GetCD3Name) .def ("SetCD3Name", &Mesh::SetCD3Name) - .def ("AddPointIdentification", [](Mesh & self, py::object pindex1, py::object pindex2, int identnr, int type) + .def ("AddPointIdentification", [](Mesh & self, py::object pindex1, py::object pindex2, int identnr, Identifications::ID_TYPE type) { if(py::extract(pindex1).check() && py::extract(pindex2).check()) { self.GetIdentifications().Add (py::extract(pindex1)(), py::extract(pindex2)(), identnr); - self.GetIdentifications().SetType(identnr, Identifications::ID_TYPE(type)); // type = 2 ... periodic + self.GetIdentifications().SetType(identnr, type); // type = 2 ... periodic } }, //py::default_call_policies(), py::arg("pid1"), py::arg("pid2"), py::arg("identnr"), - py::arg("type")) + py::arg("type")=Identifications::PERIODIC) .def("IdentifyPeriodicBoundaries", &Mesh::IdentifyPeriodicBoundaries, py::arg("face1"), py::arg("face2"), py::arg("mapping"), py::arg("point_tolerance") = -1.) .def("GetNrIdentifications", [](Mesh& self) diff --git a/libsrc/occ/CMakeLists.txt b/libsrc/occ/CMakeLists.txt index 3219f07c..68eaf33f 100644 --- a/libsrc/occ/CMakeLists.txt +++ b/libsrc/occ/CMakeLists.txt @@ -1,9 +1,12 @@ +if(USE_OCC) + add_definitions(-DNGINTERFACE_EXPORTS) add_library(occ ${NG_LIB_TYPE} Partition_Inter2d.cxx Partition_Inter3d.cxx Partition_Loop.cxx Partition_Loop2d.cxx Partition_Loop3d.cxx Partition_Spliter.cxx occgenmesh.cpp occgeom.cpp occmeshsurf.cpp python_occ.cpp python_occ_basic.cpp python_occ_shapes.cpp + occ_face.cpp occ_edge.cpp occ_vertex.cpp occ_utils.cpp ) if(USE_GUI) add_library(occvis ${NG_LIB_TYPE} vsocc.cpp) @@ -14,11 +17,11 @@ target_link_libraries(occ PUBLIC ngcore PRIVATE "$ +#include +#include + +#include "occ_edge.hpp" +#include "occgeom.hpp" + +namespace netgen +{ + OCCEdge::OCCEdge(TopoDS_Shape edge_) + : tedge(edge_.TShape()), + edge(TopoDS::Edge(edge_)) + { + curve = BRep_Tool::Curve(edge, s0, s1); + BRepGProp::LinearProperties(edge, props); + + auto verts = GetVertices(edge); + if(verts.size() != 2) + throw Exception("OCC edge does not have 2 vertices"); + + start = OCCVertex(verts[0]); + end = OCCVertex(verts[1]); + + // swap start/end if necessary + double d00 = Dist(GetPoint(0), start.GetPoint()); + double d01 = Dist(GetPoint(0), end.GetPoint()); + if(d01 < d00) + swap(start, end); + } + + const GeometryVertex& OCCEdge::GetStartVertex() const + { + return start; + } + + const GeometryVertex& OCCEdge::GetEndVertex() const + { + return end; + } + + double OCCEdge::GetLength() const + { + return props.Mass(); + } + + Point<3> OCCEdge::GetCenter() const + { + return occ2ng( props.CentreOfMass() ); + } + + Point<3> OCCEdge::GetPoint(double t) const + { + return occ2ng( curve->Value(s0+t*(s1-s0)) ); + } + + double OCCEdge::CalcStep(double t, double sag) const + { + throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__)); + } + + size_t OCCEdge::GetHash() const + { + return reinterpret_cast(tedge.get()); + } + + void OCCEdge::ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const + { + auto pnt = ng2occ(p); + GeomAPI_ProjectPointOnCurve proj(pnt, curve); + pnt = proj.NearestPoint(); + if(gi) + gi->dist = proj.LowerDistanceParameter(); + p = occ2ng(pnt); + } + + Vec<3> OCCEdge::GetTangent(double t) const + { + t = s0 + t*(s1-s0); + gp_Pnt p; + gp_Vec v; + curve->D1(t, p, v); + return occ2ng(v); + } +} diff --git a/libsrc/occ/occ_edge.hpp b/libsrc/occ/occ_edge.hpp new file mode 100644 index 00000000..871fa252 --- /dev/null +++ b/libsrc/occ/occ_edge.hpp @@ -0,0 +1,44 @@ +#ifndef FILE_OCC_EDGE_INCLUDED +#define FILE_OCC_EDGE_INCLUDED + +#include +#include +#include +#include +#include + +#include "occ_vertex.hpp" +#include "meshing.hpp" + +namespace netgen +{ + class OCCEdge : public GeometryEdge + { + T_Shape tedge; + TopoDS_Edge edge; + Handle(Geom_Curve) curve; + double s0, s1; + GProp_GProps props; + + OCCVertex start; + OCCVertex end; + + public: + OCCEdge(TopoDS_Shape edge_); + + auto Shape() const { return edge; } + T_Shape TShape() const { return tedge; } + + const GeometryVertex& GetStartVertex() const override; + const GeometryVertex& GetEndVertex() const override; + double GetLength() const override; + Point<3> GetCenter() const override; + Point<3> GetPoint(double t) const override; + double CalcStep(double t, double sag) const override; + size_t GetHash() const override; + void ProjectPoint(Point<3>& p, EdgePointGeomInfo* gi) const override; + Vec<3> GetTangent(double t) const override; + }; +} + +#endif // FILE_OCCEDGE_INCLUDED diff --git a/libsrc/occ/occ_face.cpp b/libsrc/occ/occ_face.cpp new file mode 100644 index 00000000..6988808d --- /dev/null +++ b/libsrc/occ/occ_face.cpp @@ -0,0 +1,254 @@ +#include +#include +#include + +#include "occ_edge.hpp" +#include "occ_face.hpp" +#include "occgeom.hpp" + +namespace netgen +{ + OCCFace::OCCFace(TopoDS_Shape dshape) + : tface(dshape.TShape()), + face(TopoDS::Face(dshape)) + { + BRepGProp::LinearProperties(face, props); + bbox = ::netgen::GetBoundingBox(face); + + surface = BRep_Tool::Surface(face); + shape_analysis = new ShapeAnalysis_Surface( surface ); + tolerance = BRep_Tool::Tolerance( face ); + } + + size_t OCCFace::GetNBoundaries() const + { + return 0; + } + + size_t OCCFace::GetHash() const + { + return reinterpret_cast(tface.get()); + } + + Point<3> OCCFace::GetCenter() const + { + return occ2ng( props.CentreOfMass() ); + } + + Array OCCFace::GetBoundary(const Mesh& mesh) const + { + auto & geom = dynamic_cast(*mesh.GetGeometry()); + + auto n_edges = geom.edge_map.size(); + constexpr int UNUSED = 0; + constexpr int FORWARD = 1; + constexpr int REVERSED = 2; + constexpr int BOTH = 3; + + Array edge_orientation(n_edges); + edge_orientation = UNUSED; + + Array curve_on_face[BOTH]; + curve_on_face[FORWARD].SetSize(n_edges); + curve_on_face[REVERSED].SetSize(n_edges); + + Array edge_on_face[BOTH]; + edge_on_face[FORWARD].SetSize(n_edges); + edge_on_face[REVERSED].SetSize(n_edges); + + for(auto edge_ : GetEdges(face)) + { + auto edge = TopoDS::Edge(edge_); + if(geom.edge_map.count(edge.TShape())==0) + continue; + auto edgenr = geom.edge_map[edge.TShape()]; + auto & orientation = edge_orientation[edgenr]; + double s0, s1; + auto cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1); + if(edge.Orientation() == TopAbs_FORWARD) + { + curve_on_face[FORWARD][edgenr] = cof; + orientation += FORWARD; + edge_on_face[FORWARD][edgenr] = edge; + } + if(edge.Orientation() == TopAbs_REVERSED) + { + curve_on_face[REVERSED][edgenr] = cof; + orientation += REVERSED; + edge_on_face[REVERSED][edgenr] = edge; + } + + if(orientation > BOTH) + throw Exception("have edge more than twice in face " + ToString(nr) + " " + properties.GetName() + ", orientation: " + ToString(orientation)); + } + + Array boundary; + for (auto seg : mesh.LineSegments()) + { + auto edgenr = seg.epgeominfo[0].edgenr; + auto orientation = edge_orientation[edgenr]; + + if(orientation == UNUSED) + continue; + + for(const auto ORIENTATION : {FORWARD, REVERSED}) + { + if((orientation & ORIENTATION) == 0) + continue; + + // auto cof = curve_on_face[ORIENTATION][edgenr]; + auto edge = edge_on_face[ORIENTATION][edgenr]; + OCCEdge gedge(edge); + double s0, s1; + auto cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1); + + for(auto i : Range(2)) + { + Point<3> p = mesh[seg[i]]; + gedge.ProjectPoint(p, &seg.epgeominfo[i]); + } + double s[2] = { seg.epgeominfo[0].dist, seg.epgeominfo[1].dist }; + + // fixes normal-vector roundoff problem when endpoint is cone-tip + double delta = s[1]-s[0]; + s[0] += 1e-10*delta; + s[1] -= 1e-10*delta; + + for(auto i : Range(2)) + { + auto uv = cof->Value(s[i]); + seg.epgeominfo[i].u = uv.X(); + seg.epgeominfo[i].v = uv.Y(); + } + + if(ORIENTATION == REVERSED) + { + 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); + } + + boundary.Append(seg); + } + } + return boundary; + } + + PointGeomInfo OCCFace::Project(Point<3>& p) const + { + auto suval = shape_analysis->ValueOfUV(ng2occ(p), tolerance); + double u,v; + suval.Coord(u, v); + p = occ2ng(surface->Value( u, v )); + + PointGeomInfo gi; + gi.trignum = nr+1; + gi.u = u; + gi.v = v; + return gi; + } + + bool OCCFace::ProjectPointGI(Point<3>& p_, PointGeomInfo& gi) const + { + double u = gi.u; + double v = gi.v; + auto p = ng2occ(p_); + auto x = surface->Value (u,v); + + if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true; + + gp_Vec du, dv; + surface->D1(u,v,x,du,dv); + + int count = 0; + gp_Pnt xold; + gp_Vec n; + double det, lambda, mu; + + do { + count++; + + n = du^dv; + + det = Det3 (n.X(), du.X(), dv.X(), + n.Y(), du.Y(), dv.Y(), + n.Z(), du.Z(), dv.Z()); + + if (det < 1e-15) return false; + + lambda = Det3 (n.X(), p.X()-x.X(), dv.X(), + n.Y(), p.Y()-x.Y(), dv.Y(), + n.Z(), p.Z()-x.Z(), dv.Z())/det; + + mu = Det3 (n.X(), du.X(), p.X()-x.X(), + n.Y(), du.Y(), p.Y()-x.Y(), + n.Z(), du.Z(), p.Z()-x.Z())/det; + + u += lambda; + v += mu; + + xold = x; + surface->D1(u,v,x,du,dv); + + } while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50); + + // (*testout) << "FastProject count: " << count << endl; + + if (count == 50) return false; + + p_ = occ2ng(x); + + return true; + } + + Point<3> OCCFace::GetPoint(const PointGeomInfo& gi) const + { + return occ2ng(surface->Value( gi.u, gi.v )); + } + + void OCCFace::CalcEdgePointGI(const GeometryEdge& edge, + double t, + EdgePointGeomInfo& egi) const + { + throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__)); + } + + Box<3> OCCFace::GetBoundingBox() const + { + return bbox; + } + + + double OCCFace::GetCurvature(const PointGeomInfo& gi) const + { + throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__)); + } + + void OCCFace::RestrictH(Mesh& mesh, const MeshingParameters& mparam) const + { + throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__)); + } + + Vec<3> OCCFace::GetNormal(const Point<3>& p, const PointGeomInfo* gi) const + { + PointGeomInfo gi_; + if(gi==nullptr) + { + auto p_ = p; + gi_ = Project(p_); + gi = &gi_; + } + + gp_Pnt pnt; + gp_Vec du, dv; + surface->D1(gi->u,gi->v,pnt,du,dv); + auto n = Cross (occ2ng(du), occ2ng(dv)); + n.Normalize(); + if (face.Orientation() == TopAbs_REVERSED) + n *= -1; + return n; + } + + +} diff --git a/libsrc/occ/occ_face.hpp b/libsrc/occ/occ_face.hpp new file mode 100644 index 00000000..63cc7715 --- /dev/null +++ b/libsrc/occ/occ_face.hpp @@ -0,0 +1,49 @@ +#ifndef FILE_OCC_FACE_INCLUDED +#define FILE_OCC_FACE_INCLUDED + +#include +#include +#include +#include + +#include "occ_vertex.hpp" +#include "meshing.hpp" + +namespace netgen +{ + class OCCFace : public GeometryFace + { + T_Shape tface; + TopoDS_Face face; + GProp_GProps props; + Box<3> bbox; + + Handle( Geom_Surface ) surface; + Handle( ShapeAnalysis_Surface ) shape_analysis; + double tolerance; + + public: + OCCFace(TopoDS_Shape dshape); + + T_Shape TShape() { return tface; } + + size_t GetHash() const override; + Point<3> GetCenter() const override; + virtual size_t GetNBoundaries() const override; + virtual Array GetBoundary(const Mesh& mesh) const override; + virtual PointGeomInfo Project(Point<3>& p) const override; + virtual bool ProjectPointGI(Point<3>& p, PointGeomInfo& gi) const override; + virtual Point<3> GetPoint(const PointGeomInfo& gi) const override; + virtual void CalcEdgePointGI(const GeometryEdge& edge, + double t, + EdgePointGeomInfo& egi) const override; + virtual Box<3> GetBoundingBox() const override; + + virtual double GetCurvature(const PointGeomInfo& gi) const override; + + virtual void RestrictH(Mesh& mesh, const MeshingParameters& mparam) const override; + virtual Vec<3> GetNormal(const Point<3>& p, const PointGeomInfo* gi = nullptr) const override; + }; +} + +#endif // FILE_OCC_FACE_INCLUDED diff --git a/libsrc/occ/occ_solid.hpp b/libsrc/occ/occ_solid.hpp new file mode 100644 index 00000000..869df215 --- /dev/null +++ b/libsrc/occ/occ_solid.hpp @@ -0,0 +1,27 @@ +#ifndef FILE_OCC_SOLID_INCLUDED +#define FILE_OCC_SOLID_INCLUDED + +#include +#include + +#include "meshing.hpp" + +namespace netgen +{ + class OCCSolid : public GeometrySolid + { + T_Shape tsolid; + TopoDS_Solid solid; + + public: + OCCSolid(TopoDS_Shape dshape) + : tsolid(dshape.TShape()), + solid(TopoDS::Solid(dshape)) + { } + + T_Shape TShape() { return tsolid; } + size_t GetHash() const override { return reinterpret_cast(tsolid.get()); } + }; +} + +#endif // FILE_OCC_SOLID_INCLUDED diff --git a/libsrc/occ/occ_utils.cpp b/libsrc/occ/occ_utils.cpp new file mode 100644 index 00000000..8b747137 --- /dev/null +++ b/libsrc/occ/occ_utils.cpp @@ -0,0 +1,24 @@ +#include +#include +#include + +#include "occ_utils.hpp" + +namespace netgen +{ + Point<3> occ2ng (Handle(TopoDS_TShape) shape) + { + return occ2ng( Handle(BRep_TVertex)::DownCast(shape)->Pnt() ); + } + + Box<3> GetBoundingBox( const TopoDS_Shape & shape ) + { + Bnd_Box bb; +#if OCC_VERSION_HEX < 0x070000 + BRepBndLib::Add (shape, bb); +#else + BRepBndLib::Add (shape, bb, true); +#endif + return {occ2ng(bb.CornerMin()), occ2ng(bb.CornerMax())}; + } +} diff --git a/libsrc/occ/occ_utils.hpp b/libsrc/occ/occ_utils.hpp new file mode 100644 index 00000000..52930206 --- /dev/null +++ b/libsrc/occ/occ_utils.hpp @@ -0,0 +1,267 @@ +#ifndef FILE_OCC_UTILS_INCLUDED +#define FILE_OCC_UTILS_INCLUDED + +#include +#include +#include +#include +#include +#include +#include + +#include "meshing.hpp" + +#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4 +#define OCC_HAVE_DUMP_JSON +#endif + +namespace netgen +{ + typedef Handle(TopoDS_TShape) T_Shape; + + inline Point<3> occ2ng (const gp_Pnt & p) + { + return Point<3> (p.X(), p.Y(), p.Z()); + } + + inline Point<2> occ2ng (const gp_Pnt2d & p) + { + return Point<2> (p.X(), p.Y()); + } + + inline Vec<3> occ2ng (const gp_Vec & v) + { + return Vec<3> (v.X(), v.Y(), v.Z()); + } + + DLL_HEADER Point<3> occ2ng (T_Shape shape); + + inline Point<3> occ2ng (const TopoDS_Shape & s) + { + return occ2ng(s.TShape()); + } + + inline Point<3> occ2ng (const TopoDS_Vertex & v) + { + return occ2ng (BRep_Tool::Pnt (v)); + } + + inline gp_Pnt ng2occ (const Point<3> & p) + { + return gp_Pnt(p(0), p(1), p(2)); + } + + DLL_HEADER Box<3> GetBoundingBox( const TopoDS_Shape & shape ); + + class OCCIdentification + { + public: + T_Shape from; + T_Shape to; + Transformation<3> trafo; + string name; + Identifications::ID_TYPE type; + bool opposite_direction; + }; + + + class MyExplorer + { + class Iterator + { + TopExp_Explorer exp; + public: + Iterator (TopoDS_Shape ashape, TopAbs_ShapeEnum atoFind, TopAbs_ShapeEnum atoAvoid) + : exp(ashape, atoFind, atoAvoid) { } + auto operator*() { return exp.Current(); } + Iterator & operator++() { exp.Next(); return *this; } + bool operator!= (nullptr_t nu) { return exp.More(); } + }; + + public: + TopoDS_Shape shape; + TopAbs_ShapeEnum toFind; + TopAbs_ShapeEnum toAvoid; + MyExplorer (TopoDS_Shape ashape, TopAbs_ShapeEnum atoFind, TopAbs_ShapeEnum atoAvoid = TopAbs_SHAPE) + : shape(ashape), toFind(atoFind), toAvoid(atoAvoid) { ; } + Iterator begin() { return Iterator(shape, toFind, toAvoid); } + auto end() { return nullptr; } + }; + + inline auto Explore (TopoDS_Shape shape, TopAbs_ShapeEnum toFind, TopAbs_ShapeEnum toAvoid = TopAbs_SHAPE) + { + return MyExplorer (shape, toFind, toAvoid); + } + + + class IndexMapIterator + { + class Iterator + { + const TopTools_IndexedMapOfShape & indmap; + int i; + public: + Iterator (const TopTools_IndexedMapOfShape & aindmap, int ai) + : indmap(aindmap), i(ai) { ; } + auto operator*() { return tuple(i, indmap(i)); } + Iterator & operator++() { i++; return *this; } + bool operator!= (const Iterator & i2) { return i != i2.i; } + }; + + public: + const TopTools_IndexedMapOfShape & indmap; + IndexMapIterator (const TopTools_IndexedMapOfShape & aindmap) : indmap(aindmap) { } + Iterator begin() { return Iterator(indmap, 1); } + Iterator end() { return Iterator(indmap, indmap.Extent()+1); } + }; + + inline auto Enumerate (const TopTools_IndexedMapOfShape & indmap) + { + return IndexMapIterator(indmap); + } + + struct ShapeLess + { + bool operator() (const TopoDS_Shape& s1, const TopoDS_Shape& s2) const + { + return s1.TShape() < s2.TShape(); + } + }; + + class ListOfShapes : public std::vector + { + public: + DLL_HEADER TopoDS_Shape Max(gp_Vec dir); + DLL_HEADER TopoDS_Shape Nearest(gp_Pnt pnt); + DLL_HEADER ListOfShapes SubShapes(TopAbs_ShapeEnum type) const; + + ListOfShapes Solids() const + { + return SubShapes(TopAbs_SOLID); + } + ListOfShapes Faces() const + { + return SubShapes(TopAbs_FACE); + } + ListOfShapes Edges() const + { + return SubShapes(TopAbs_EDGE); + } + ListOfShapes Vertices() const + { + return SubShapes(TopAbs_VERTEX); + } + + ListOfShapes operator*(const ListOfShapes& other) const + { + ListOfShapes common; + for(const auto& shape : (*this)) + for(const auto& shape_o : other) + if(shape.IsSame(shape_o)) + common.push_back(shape); + return common; + } + }; + + inline ListOfShapes GetSolids(const TopoDS_Shape & shape) + { + ListOfShapes sub; + for (TopExp_Explorer e(shape, TopAbs_SOLID); e.More(); e.Next()) + sub.push_back(e.Current()); + return sub; + } + + inline ListOfShapes GetFaces(const TopoDS_Shape & shape) + { + ListOfShapes sub; + for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next()) + sub.push_back(e.Current()); + return sub; + } + + inline ListOfShapes GetEdges(const TopoDS_Shape & shape) + { + ListOfShapes sub; + for (TopExp_Explorer e(shape, TopAbs_EDGE); e.More(); e.Next()) + sub.push_back(e.Current()); + return sub; + } + + inline ListOfShapes GetVertices(const TopoDS_Shape & shape) + { + ListOfShapes sub; + for (TopExp_Explorer e(shape, TopAbs_VERTEX); e.More(); e.Next()) + sub.push_back(e.Current()); + return sub; + } + + class DirectionalInterval + { + public: + gp_Vec dir; + double minval = -1e99; + double maxval = 1e99; + bool openmin = false, openmax = false; + + DirectionalInterval (gp_Vec adir) : dir(adir) { ; } + DirectionalInterval (const DirectionalInterval & i2) + : dir(i2.dir), minval(i2.minval), maxval(i2.maxval) { ; } + + DirectionalInterval operator< (double val) const + { + DirectionalInterval i2 = *this; + i2.maxval = val; + return i2; + } + + DirectionalInterval operator> (double val) const + { + DirectionalInterval i2 = *this; + i2.minval = val; + return i2; + } + + + DirectionalInterval Intersect (const DirectionalInterval & i2) + { + DirectionalInterval res = *this; + res.minval = max(res.minval, i2.minval); + res.maxval = min(res.maxval, i2.maxval); + return res; + } + + bool Contains (gp_Pnt p, double eps = 1e-8) + { + // cout << "Contains point " << p.X() << "," << p.Y() << "," << p.Z() << " ? " << endl; + double val = dir.X()*p.X() + dir.Y()*p.Y() + dir.Z() * p.Z(); + // cout << "minval = " << minval << ", val = " << val << " maxval = " << maxval << endl; + if (openmin) { + if (val < minval+eps) return false; + } else { + if (val < minval-eps) return false; + } + if (openmax) { + if (val > maxval-eps) return false; + } else { + if (val > maxval+eps) return false; + } + return true; + } + }; + + + inline gp_Pnt Center (TopoDS_Shape shape) + { + GProp_GProps props; + switch (shape.ShapeType()) + { + case TopAbs_FACE: + BRepGProp::SurfaceProperties (shape, props); break; + default: + BRepGProp::LinearProperties(shape, props); + } + return props.CentreOfMass(); + } + +} +#endif // FILE_OCC_UTILS_INCLUDED diff --git a/libsrc/occ/occ_vertex.cpp b/libsrc/occ/occ_vertex.cpp new file mode 100644 index 00000000..0f114651 --- /dev/null +++ b/libsrc/occ/occ_vertex.cpp @@ -0,0 +1,25 @@ +#include +#include + +#include "occ_vertex.hpp" + +namespace netgen +{ + + OCCVertex::OCCVertex( TopoDS_Shape s ) + : vertex(TopoDS::Vertex(s)), + tvertex(s.TShape()) + { + p = occ2ng(vertex); + } + + Point<3> OCCVertex::GetPoint() const + { + return p; + } + + size_t OCCVertex::GetHash() const + { + return reinterpret_cast(tvertex.get()); + } +} diff --git a/libsrc/occ/occ_vertex.hpp b/libsrc/occ/occ_vertex.hpp new file mode 100644 index 00000000..c1d6a488 --- /dev/null +++ b/libsrc/occ/occ_vertex.hpp @@ -0,0 +1,28 @@ +#ifndef FILE_OCC_VERTEX_INCLUDED +#define FILE_OCC_VERTEX_INCLUDED + +#include +#include + +#include "meshing.hpp" +#include "occ_utils.hpp" + +namespace netgen +{ + class OCCVertex : public GeometryVertex + { + TopoDS_Vertex vertex; + T_Shape tvertex; + Point<3> p; + + public: + OCCVertex( ) = default; + OCCVertex( TopoDS_Shape s ); + ~OCCVertex() {} + Point<3> GetPoint() const override; + size_t GetHash() const override; + T_Shape TShape() { return tvertex; } + }; +} + +#endif // FILE_OCC_VERTEX_INCLUDED diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp index f4b8c160..6f87e389 100644 --- a/libsrc/occ/occgenmesh.cpp +++ b/libsrc/occ/occgenmesh.cpp @@ -1,14 +1,27 @@ #ifdef OCCGEOMETRY #include -#include #include +#include "occgeom.hpp" +#include "occmeshsurf.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include namespace netgen { -#include "occmeshsurf.hpp" #define TCL_OK 0 #define TCL_ERROR 1 @@ -218,890 +231,239 @@ namespace netgen } } - - - void DivideEdge (TopoDS_Edge & edge, NgArray & ps, - Array & params, Mesh & mesh, - const MeshingParameters & mparam) + bool OCCMeshFace (const OCCGeometry & geom, Mesh & mesh, FlatArray glob2loc, + const MeshingParameters & mparam, int nr, int projecttype, bool delete_on_failure) { - double s0, s1; - int nsubedges = 1; - gp_Pnt pnt, oldpnt; - double svalue[DIVIDEEDGESECTIONS]; - - GProp_GProps system; - BRepGProp::LinearProperties(edge, system); - double L = system.Mass(); - - Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); - - double hvalue[DIVIDEEDGESECTIONS+1]; - hvalue[0] = 0; - pnt = c->Value(s0); - - int tmpVal = (int)(DIVIDEEDGESECTIONS); - - for (int i = 1; i <= tmpVal; i++) + auto k = nr+1; + if(1==0 && !geom.fvispar[k-1].IsDrawable()) { - oldpnt = pnt; - pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0)); - hvalue[i] = hvalue[i-1] + - 1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))* - pnt.Distance(oldpnt); - - //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) - // << " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl; + (*testout) << "ignoring face " << k << endl; + cout << "ignoring face " << k << endl; + return true; } - // nsubedges = int(ceil(hvalue[DIVIDEEDGESECTIONS])); - nsubedges = max (1, int(floor(hvalue[DIVIDEEDGESECTIONS]+0.5))); + // if(master_faces[k]!=k) + // continue; - ps.SetSize(nsubedges-1); - params.SetSize(nsubedges+1); + (*testout) << "mesh face " << k << endl; + multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL); + geom.facemeshstatus[k-1] = -1; - int i = 1; - int i1 = 0; - do - { - if (hvalue[i1]/hvalue[DIVIDEEDGESECTIONS]*nsubedges >= i) - { - params[i] = s0+(i1/double(DIVIDEEDGESECTIONS))*(s1-s0); - pnt = c->Value(params[i]); - ps[i-1] = MeshPoint (Point3d(pnt.X(), pnt.Y(), pnt.Z())); - i++; - } - i1++; - if (i1 > DIVIDEEDGESECTIONS) - { - nsubedges = i; - ps.SetSize(nsubedges-1); - params.SetSize(nsubedges+1); - cout << "divide edge: local h too small" << endl; - } - } while (i < nsubedges); + FaceDescriptor & fd = mesh.GetFaceDescriptor(k); + auto face = TopoDS::Face(geom.fmap(k)); + auto fshape = face.TShape(); - params[0] = s0; - params[nsubedges] = s1; + int oldnf = mesh.GetNSE(); - if (params[nsubedges] <= params[nsubedges-1]) - { - cout << "CORRECTED" << endl; - ps.SetSize (nsubedges-2); - params.SetSize (nsubedges); - params[nsubedges] = s1; - } - } + Box<3> bb = geom.GetBoundingBox(); - - - - void OCCFindEdges (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam) - { - static Timer t("OCCFindEdges"); RegionTimer r(t); - static Timer tsearch("OCCFindEdges - search point"); - const char * savetask = multithread.task; - multithread.task = "Edge meshing"; - - (*testout) << "edge meshing" << endl; - - int nvertices = geom.vmap.Extent(); - int nedges = geom.emap.Extent(); - - Array> alledgepnums(nedges); - Array> alledgeparams(nedges); + // int projecttype = PLANESPACE; + // int projecttype = PARAMETERSPACE; - (*testout) << "nvertices = " << nvertices << endl; - (*testout) << "nedges = " << nedges << endl; + static Timer tinit("init"); + tinit.Start(); + Meshing2OCCSurfaces meshing(geom, face, bb, projecttype, mparam); + tinit.Stop(); - double eps = 1e-6 * geom.GetBoundingBox().Diam(); - tsearch.Start(); - for (auto [i,vshape] : Enumerate(geom.vmap)) - { - TopoDS_Vertex vertex = TopoDS::Vertex(vshape); - gp_Pnt pnt = BRep_Tool::Pnt (vertex); + static Timer tprint("print"); + tprint.Start(); + if (meshing.GetProjectionType() == PLANESPACE) + PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (plane space projection)"); + else + PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (parameter space projection)"); + tprint.Stop(); - mesh.AddPoint (occ2ng(pnt)); - - double hpref = OCCGeometry::global_shape_properties[vertex.TShape()].hpref; - mesh.Points().Last().Singularity(hpref); - - double maxh = OCCGeometry::global_shape_properties[vertex.TShape()].maxh; - mesh.RestrictLocalH (occ2ng(pnt), maxh); - } - tsearch.Stop(); - (*testout) << "different vertices = " << mesh.GetNP() << endl; - - // int first_ep = mesh.GetNP()+1; - // PointIndex first_ep = mesh.Points().End(); - PointIndex first_ep = *mesh.Points().Range().end(); - auto vertexrange = mesh.Points().Range(); - - NgArray face2solid[2]; - for (int i = 0; i < 2; i++) - { - face2solid[i].SetSize (geom.fmap.Extent()); - face2solid[i] = 0; - } - - /* - int solidnr = 0; - for (TopExp_Explorer exp0(geom.shape, TopAbs_SOLID); exp0.More(); exp0.Next()) - { - solidnr++; - for (TopExp_Explorer exp1(exp0.Current(), TopAbs_FACE); exp1.More(); exp1.Next()) - { - TopoDS_Face face = TopoDS::Face(exp1.Current()); - int facenr = geom.fmap.FindIndex(face); - if(facenr < 1) continue; - - if (face2solid[0][facenr-1] == 0) - face2solid[0][facenr-1] = solidnr; - else - face2solid[1][facenr-1] = solidnr; - } - } - */ - int solidnr = 0; - for (auto solid : Explore(geom.shape, TopAbs_SOLID)) - { - solidnr++; - for (auto face : Explore (solid, TopAbs_FACE)) - if (geom.fmap.Contains(face)) - { - int facenr = geom.fmap.FindIndex(face); - if (face2solid[0][facenr-1] == 0) - face2solid[0][facenr-1] = solidnr; - else - face2solid[1][facenr-1] = solidnr; - } - } + // Meshing2OCCSurfaces meshing(f2, bb); + // meshing.SetStartTime (starttime); + //(*testout) << "Face " << k << endl << endl; + auto segments = geom.GetFace(k-1).GetBoundary(mesh); - /* - int total = 0; - for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++) - for (TopExp_Explorer exp2(geom.fmap(i3), TopAbs_WIRE); exp2.More(); exp2.Next()) - for (TopExp_Explorer exp3(exp2.Current(), TopAbs_EDGE); exp3.More(); exp3.Next()) - total++; - */ - int total = 0; - for (auto [i3, face] : Enumerate(geom.fmap)) - for (auto wire : Explore(face, TopAbs_WIRE)) - for (auto edge : Explore(wire, TopAbs_EDGE)) - total++; - - - - int facenr = 0; - // int edgenr = mesh.GetNSeg(); - - (*testout) << "faces = " << geom.fmap.Extent() << endl; - int curr = 0; - - /* - for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++) + if (meshing.GetProjectionType() == PLANESPACE) { - TopoDS_Face face = TopoDS::Face(geom.fmap(i3)); - */ - for (auto [i3,faceshape] : Enumerate(geom.fmap)) - { - TopoDS_Face face = TopoDS::Face(faceshape); - facenr = geom.fmap.FindIndex (face); // sollte doch immer == i3 sein ??? JS - if (facenr != i3) - cout << "info: facenr != i3, no problem, but please report to developers" << endl; - - int solidnr0 = face2solid[0][i3-1]; - int solidnr1 = face2solid[1][i3-1]; + static Timer t("MeshSurface: Find edges and points - Physical"); RegionTimer r(t); + int cntp = 0; + glob2loc = 0; - /* auskommentiert am 3.3.05 von robert - for (exp2.Init (geom.somap(solidnr0), TopAbs_FACE); exp2.More(); exp2.Next()) - { - TopoDS_Face face2 = TopoDS::Face(exp2.Current()); - if (geom.fmap.FindIndex(face2) == facenr) - { - // if (face.Orientation() != face2.Orientation()) swap (solidnr0, solidnr1); - } - } - */ - - mesh.AddFaceDescriptor (FaceDescriptor(facenr, solidnr0, solidnr1, 0)); - - Vec<4> col = OCCGeometry::global_shape_properties[face.TShape()].col.value_or(Vec<4>(0,1,0,1)); - mesh.GetFaceDescriptor(facenr).SetSurfColour(col); - - if(auto & opt_name = geom.fprops[facenr-1]->name) - mesh.GetFaceDescriptor(facenr).SetBCName(*opt_name); - else - mesh.GetFaceDescriptor(facenr).SetBCName("bc_"+ToString(facenr)); - mesh.GetFaceDescriptor(facenr).SetBCProperty(facenr); - // ACHTUNG! STIMMT NICHT ALLGEMEIN (RG) - // kA was RG damit meinte - - - Handle(Geom_Surface) occface = BRep_Tool::Surface(face); + for (Segment & seg : segments) + // if (seg.si == k) + for (int j = 0; j < 2; j++) + { + PointIndex pi = seg[j]; + if (glob2loc[pi] == 0) + { + meshing.AddPoint (mesh.Point(pi), pi); + cntp++; + glob2loc[pi] = cntp; + } + } /* - for (TopExp_Explorer exp2 (face, TopAbs_WIRE); exp2.More(); exp2.Next()) + for (int i = 1; i <= mesh.GetNSeg(); i++) { - TopoDS_Shape wire = exp2.Current(); + Segment & seg = mesh.LineSegment(i); */ - for (auto wire : MyExplorer (face, TopAbs_WIRE)) - { - // for (TopExp_Explorer exp3 (wire, TopAbs_EDGE); exp3.More(); exp3.Next()) - for (auto edgeshape : MyExplorer (wire, TopAbs_EDGE)) - { - TopoDS_Edge edge = TopoDS::Edge(edgeshape); - curr++; - (*testout) << "edge nr " << curr << endl; - multithread.percent = 100 * curr / double (total); - if (multithread.terminate) return; - - // TopoDS_Edge edge = TopoDS::Edge (exp3.Current()); - if (BRep_Tool::Degenerated(edge)) - { - //(*testout) << "ignoring degenerated edge" << endl; - continue; - } - - if (!geom.emap.Contains(edge)) continue; - - if (geom.vmap.FindIndex(TopExp::FirstVertex (edge)) == - geom.vmap.FindIndex(TopExp::LastVertex (edge))) - { - GProp_GProps system; - BRepGProp::LinearProperties(edge, system); - - if (system.Mass() < eps) - { - cout << "ignoring edge " << geom.emap.FindIndex (edge) - << ". closed edge with length < " << eps << endl; - continue; - } - } - - double s0, s1; - Handle(Geom2d_Curve) cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1); - - int geomedgenr = geom.emap.FindIndex(edge); - Array pnums; - Array params; - - - // check for identifications - bool copy_identified = false; - if (auto it = geom.identifications.find(edge.TShape()); it != geom.identifications.end()) - for (auto & ids : it->second) - { - cout << "edge has identification with trafo " << ids.name << ", inv = " << ids.inverse << endl; - int otherind = geom.emap.FindIndex(ids.other); - Array othersegs; - for (auto & seg : mesh.LineSegments()) - if (seg.edgenr == otherind) - othersegs.Append (seg); - - if (othersegs.Size()) - { - cout << "other has already segs" << endl; - copy_identified = true; - - Array pnums_other; - pnums_other.Append (othersegs[0][0]); - for (auto & seg : othersegs) - pnums_other.Append (seg[1]); - - auto inv = ids.trafo.CalcInverse(); - // for (auto & pi : pnums) - for (auto oi : Range(pnums_other)) - { - PointIndex piother = pnums_other[pnums_other.Size()-oi-1]; - Point<3> pother = mesh[piother]; - Point<3> p = inv(pother); - - bool found = false; - PointIndex pi; - for (PointIndex piv : vertexrange) - if (Dist2 (mesh[piv], p) < eps*eps) - { - pi = piv; - found = true; - } - - if (!found) - pi = mesh.AddPoint (p); - - // params.Add ( find parameter p ); - double s0, s1; - Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, s0, s1); - - GeomAPI_ProjectPointOnCurve proj(ng2occ(p), curve); - params.Append (proj.LowerDistanceParameter()); - pnums.Append (pi); - mesh.GetIdentifications().Add (pi, piother, geomedgenr); - - } - mesh.GetIdentifications().SetType(geomedgenr,Identifications::PERIODIC); - - copy_identified = true; - break; - } - } - - - if (!copy_identified) - { - - - if (alledgepnums[geomedgenr-1].Size()) - { - pnums = alledgepnums[geomedgenr-1]; - params = alledgeparams[geomedgenr-1]; - } - else - { - NgArray mp; - DivideEdge (edge, mp, params, mesh, mparam); - - pnums.SetSize(mp.Size()+2); - if (!merge_solids) - { - pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1; - pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1; - } - else - { - Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge))); - Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge))); - - pnums[0] = PointIndex::INVALID; - pnums.Last() = PointIndex::INVALID; - for (PointIndex pi : vertexrange) - { - if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi; - if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi; - } - } - - for (size_t i = 1; i <= mp.Size(); i++) - { - bool exists = false; - tsearch.Start(); - - /* - // for (PointIndex j = first_ep; j < mesh.Points().End(); j++) - for (PointIndex j = first_ep; j < *mesh.Points().Range().end(); j++) - if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps) - { - exists = true; - pnums[i] = j; - break; - } - */ - tsearch.Stop(); - - if (!exists) - pnums[i] = mesh.AddPoint (mp[i-1]); - } - } - - - alledgepnums[geomedgenr-1] = pnums; - alledgeparams[geomedgenr-1] = params; - } - - auto name = geom.eprops[geom.emap.FindIndex(edge)-1]->name.value_or(""); - mesh.SetCD2Name(geomedgenr, name); - - (*testout) << "NP = " << mesh.GetNP() << endl; - //(*testout) << pnums[pnums.Size()-1] << endl; - - double hpref = OCCGeometry::global_shape_properties[edge.TShape()].hpref; - - // for (size_t i = 1; i <= mp.Size()+1; i++) - for (size_t i = 1; i < pnums.Size(); i++) - { - // edgenr++; - Segment seg; - - seg[0] = pnums[i-1]; - seg[1] = pnums[i]; - // seg.edgenr = edgenr; - seg.edgenr = geomedgenr; - seg.si = facenr; - seg.epgeominfo[0].dist = params[i-1]; - seg.epgeominfo[1].dist = params[i]; - seg.epgeominfo[0].edgenr = geomedgenr; - seg.epgeominfo[1].edgenr = geomedgenr; - - double s0 = params[i-1]; - double s1 = params[i]; - double delta = s1-s0; - s0 += 1e-10*delta; // fixes normal-vector roundoff problem when endpoint is cone-tip - s1 -= 1e-10*delta; - gp_Pnt2d p2d1 = cof->Value(s0); - gp_Pnt2d p2d2 = cof->Value(s1); - seg.epgeominfo[0].u = p2d1.X(); - seg.epgeominfo[0].v = p2d1.Y(); - seg.epgeominfo[1].u = p2d2.X(); - seg.epgeominfo[1].v = p2d2.Y(); - - seg.singedge_left = hpref; - seg.singedge_right = hpref; - /* - if (occface->IsUPeriodic()) - { - cout << "U Periodic" << endl; - if (fabs(seg.epgeominfo[1].u-seg.epgeominfo[0].u) > - fabs(seg.epgeominfo[1].u- - (seg.epgeominfo[0].u-occface->UPeriod()))) - seg.epgeominfo[0].u = p2d.X()+occface->UPeriod(); - - if (fabs(seg.epgeominfo[1].u-seg.epgeominfo[0].u) > - fabs(seg.epgeominfo[1].u- - (seg.epgeominfo[0].u+occface->UPeriod()))) - seg.epgeominfo[0].u = p2d.X()-occface->UPeriod(); - } - - if (occface->IsVPeriodic()) - { - cout << "V Periodic" << endl; - if (fabs(seg.epgeominfo[1].v-seg.epgeominfo[0].v) > - fabs(seg.epgeominfo[1].v- - (seg.epgeominfo[0].v-occface->VPeriod()))) - seg.epgeominfo[0].v = p2d.Y()+occface->VPeriod(); - - if (fabs(seg.epgeominfo[1].v-seg.epgeominfo[0].v) > - fabs(seg.epgeominfo[1].v- - (seg.epgeominfo[0].v+occface->VPeriod()))) - seg.epgeominfo[0].v = p2d.Y()-occface->VPeriod(); - } - */ - - if (edge.Orientation() == TopAbs_REVERSED) - { - 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); - - //edgesegments[geomedgenr-1]->Append(mesh.GetNSeg()); - - } - } - } - } - - // for(i=1; i<=mesh.GetNSeg(); i++) - // (*testout) << "edge " << mesh.LineSegment(i).edgenr << " face " << mesh.LineSegment(i).si - // << " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl; - // exit(10); - - mesh.CalcSurfacesOfNode(); - multithread.task = savetask; - } - - - - - void OCCMeshSurface (const OCCGeometry & geom, Mesh & mesh, - const MeshingParameters & mparam) - { - static Timer t("OCCMeshSurface"); RegionTimer r(t); - - // int i, j, k; - // int changed; - - const char * savetask = multithread.task; - multithread.task = "Surface meshing"; - - geom.facemeshstatus = 0; - - int noldp = mesh.GetNP(); - - double starttime = GetTime(); - - NgArray glob2loc(noldp); - - //int projecttype = PARAMETERSPACE; - - int projecttype = PARAMETERSPACE; - - int notrys = 1; - - int surfmesherror = 0; - - for (int k = 1; k <= mesh.GetNFD(); k++) - { - if(1==0 && !geom.fvispar[k-1].IsDrawable()) - { - (*testout) << "ignoring face " << k << endl; - cout << "ignoring face " << k << endl; - continue; - } - - (*testout) << "mesh face " << k << endl; - multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL); - geom.facemeshstatus[k-1] = -1; - - FaceDescriptor & fd = mesh.GetFaceDescriptor(k); - - int oldnf = mesh.GetNSE(); - - Box<3> bb = geom.GetBoundingBox(); - - // int projecttype = PLANESPACE; - // int projecttype = PARAMETERSPACE; - - static Timer tinit("init"); - tinit.Start(); - Meshing2OCCSurfaces meshing(geom, TopoDS::Face(geom.fmap(k)), bb, projecttype, mparam); - tinit.Stop(); - - - static Timer tprint("print"); - tprint.Start(); - if (meshing.GetProjectionType() == PLANESPACE) - PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (plane space projection)"); - else - PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (parameter space projection)"); - tprint.Stop(); - if (surfmesherror) - cout << "Surface meshing error occurred before (in " << surfmesherror << " faces)" << endl; - - // Meshing2OCCSurfaces meshing(f2, bb); - meshing.SetStartTime (starttime); - //(*testout) << "Face " << k << endl << endl; - - - if (meshing.GetProjectionType() == PLANESPACE) - { - static Timer t("MeshSurface: Find edges and points - Physical"); RegionTimer r(t); - int cntp = 0; - glob2loc = 0; - - for (Segment & seg : mesh.LineSegments()) - if (seg.si == k) - for (int j = 0; j < 2; j++) - { - PointIndex pi = seg[j]; - if (glob2loc[pi] == 0) - { - meshing.AddPoint (mesh.Point(pi), pi); - cntp++; - glob2loc[pi] = cntp; - } - } - - /* - for (int i = 1; i <= mesh.GetNSeg(); i++) - { - Segment & seg = mesh.LineSegment(i); - */ - for (Segment & seg : mesh.LineSegments()) - if (seg.si == k) - { - PointGeomInfo gi0, gi1; - gi0.trignum = gi1.trignum = k; - gi0.u = seg.epgeominfo[0].u; - gi0.v = seg.epgeominfo[0].v; - gi1.u = seg.epgeominfo[1].u; - gi1.v = seg.epgeominfo[1].v; - - meshing.AddBoundaryElement (glob2loc[seg[0]], glob2loc[seg[1]], gi0, gi1); - } - } - else - { - static Timer t("MeshSurface: Find edges and points - Parameter"); RegionTimer r(t); - - int cntp = 0; - - /* - for (int i = 1; i <= mesh.GetNSeg(); i++) - if (mesh.LineSegment(i).si == k) - cntp+=2; - */ - for (Segment & seg : mesh.LineSegments()) - if (seg.si == k) - cntp += 2; - - NgArray gis; - - gis.SetAllocSize (cntp); - gis.SetSize (0); - - for (int i = 1; i <= mesh.GetNSeg(); i++) - { - Segment & seg = mesh.LineSegment(i); - if (seg.si == k) - { - PointGeomInfo gi0, gi1; - gi0.trignum = gi1.trignum = k; - gi0.u = seg.epgeominfo[0].u; - gi0.v = seg.epgeominfo[0].v; - gi1.u = seg.epgeominfo[1].u; - gi1.v = seg.epgeominfo[1].v; - - int locpnum[2] = {0, 0}; - - for (int j = 0; j < 2; j++) - { - PointGeomInfo gi = (j == 0) ? gi0 : gi1; - - /* - int l; - for (l = 0; l < gis.Size() && locpnum[j] == 0; l++) - { - double dist = sqr (gis[l].u-gi.u)+sqr(gis[l].v-gi.v); - - if (dist < 1e-10) - locpnum[j] = l+1; - } - - if (locpnum[j] == 0) - { - PointIndex pi = seg[j]; - meshing.AddPoint (mesh.Point(pi), pi); - - gis.SetSize (gis.Size()+1); - gis[l] = gi; - locpnum[j] = l+1; - } - */ - for (int l = 0; l < gis.Size(); l++) - { - double dist = sqr (gis[l].u-gi.u)+sqr(gis[l].v-gi.v); - if (dist < 1e-10) - { - locpnum[j] = l+1; - break; - } - } - - if (locpnum[j] == 0) - { - PointIndex pi = seg[j]; - meshing.AddPoint (mesh.Point(pi), pi); - - gis.Append (gi); - locpnum[j] = gis.Size(); - } - } - - meshing.AddBoundaryElement (locpnum[0], locpnum[1], gi0, gi1); - } - } - } - - - - - // Philippose - 15/01/2009 - double maxh = min2(geom.face_maxh[k-1], OCCGeometry::global_shape_properties[TopoDS::Face(geom.fmap(k)).TShape()].maxh); - //double maxh = mparam.maxh; - // int noldpoints = mesh->GetNP(); - int noldsurfel = mesh.GetNSE(); - - static Timer tsurfprop("surfprop"); - tsurfprop.Start(); - GProp_GProps sprops; - BRepGProp::SurfaceProperties(TopoDS::Face(geom.fmap(k)),sprops); - tsurfprop.Stop(); - meshing.SetMaxArea(2.*sprops.Mass()); - - MESHING2_RESULT res; - - // TODO: check overlap not correctly working here - MeshingParameters mparam_without_overlap = mparam; - mparam_without_overlap.checkoverlap = false; - - try { - static Timer t("GenerateMesh"); RegionTimer reg(t); - res = meshing.GenerateMesh (mesh, mparam_without_overlap, maxh, k); - } - - catch (SingularMatrixException) - { - // (*myerr) << "Singular Matrix" << endl; - res = MESHING2_GIVEUP; - } - - catch (UVBoundsException) - { - // (*myerr) << "UV bounds exceeded" << endl; - res = MESHING2_GIVEUP; - } - - projecttype = PARAMETERSPACE; - static Timer t1("rest of loop"); RegionTimer reg1(t1); - - if (res != MESHING2_OK) - { - if (notrys == 1) - { - for (SurfaceElementIndex sei = noldsurfel; sei < mesh.GetNSE(); sei++) - mesh.Delete(sei); - - mesh.Compress(); - - (*testout) << "retry Surface " << k << endl; - - k--; - // projecttype*=-1; - projecttype = PLANESPACE; - notrys++; - continue; - } - else - { - geom.facemeshstatus[k-1] = -1; - PrintError ("Problem in Surface mesh generation"); - surfmesherror++; - // throw NgException ("Problem in Surface mesh generation"); - } - } - else - { - geom.facemeshstatus[k-1] = 1; - } - - notrys = 1; - - for (SurfaceElementIndex sei = oldnf; sei < mesh.GetNSE(); sei++) - mesh[sei].SetIndex (k); - - auto n_illegal_trigs = mesh.FindIllegalTrigs(); - PrintMessage (3, n_illegal_trigs, " illegal triangles"); - } - - // ofstream problemfile("occmesh.rep"); - - // problemfile << "SURFACEMESHING" << endl << endl; - - if (surfmesherror) - { - cout << "WARNING! NOT ALL FACES HAVE BEEN MESHED" << endl; - cout << "SURFACE MESHING ERROR OCCURRED IN " << surfmesherror << " FACES:" << endl; - for (int i = 1; i <= geom.fmap.Extent(); i++) - if (geom.facemeshstatus[i-1] == -1) + // for (Segment & seg : mesh.LineSegments()) + for (Segment & seg : segments) + //if (seg.si == k) { - cout << "Face " << i << endl; - // problemfile << "problem with face " << i << endl; - // problemfile << "vertices: " << endl; - TopExp_Explorer exp0,exp1,exp2; - for ( exp0.Init(TopoDS::Face (geom.fmap(i)), TopAbs_WIRE); exp0.More(); exp0.Next() ) - { - TopoDS_Wire wire = TopoDS::Wire(exp0.Current()); - for ( exp1.Init(wire,TopAbs_EDGE); exp1.More(); exp1.Next() ) - { - TopoDS_Edge edge = TopoDS::Edge(exp1.Current()); - for ( exp2.Init(edge,TopAbs_VERTEX); exp2.More(); exp2.Next() ) - { - TopoDS_Vertex vertex = TopoDS::Vertex(exp2.Current()); - gp_Pnt point = BRep_Tool::Pnt(vertex); - // problemfile << point.X() << " " << point.Y() << " " << point.Z() << endl; - } - } - } - // problemfile << endl; + PointGeomInfo gi0, gi1; + gi0.trignum = gi1.trignum = k; + gi0.u = seg.epgeominfo[0].u; + gi0.v = seg.epgeominfo[0].v; + gi1.u = seg.epgeominfo[1].u; + gi1.v = seg.epgeominfo[1].v; + + //if(orientation & 1) + meshing.AddBoundaryElement (glob2loc[seg[0]], glob2loc[seg[1]], gi0, gi1); } - cout << endl << endl; - cout << "for more information open IGES/STEP Topology Explorer" << endl; - // problemfile.close(); - throw NgException ("Problem in Surface mesh generation"); } else { - // problemfile << "OK" << endl << endl; - // problemfile.close(); - } - - for (int i = 0; i < mesh.GetNFD(); i++) - mesh.SetBCName (i, mesh.GetFaceDescriptor(i+1).GetBCName()); - multithread.task = savetask; - } + static Timer t("MeshSurface: Find edges and points - Parameter"); RegionTimer r(t); + + int cntp = 0; - void OCCOptimizeSurface(OCCGeometry & geom, Mesh & mesh, - const MeshingParameters & mparam) - { - const char * savetask = multithread.task; - multithread.task = "Optimizing surface"; + /* + for (int i = 1; i <= mesh.GetNSeg(); i++) + if (mesh.LineSegment(i).si == k) + cntp+=2; + */ + cntp = 2*segments.Size(); + //for (Segment & seg : mesh.LineSegments()) + //if (seg.si == k) + //cntp += 2; - static Timer timer_opt2d("Optimization 2D"); - timer_opt2d.Start(); + NgArray gis; - for (int k = 1; k <= mesh.GetNFD(); k++) - { - // if (k != 42) continue; - // if (k != 36) continue; + gis.SetAllocSize (cntp); + gis.SetSize (0); - // (*testout) << "optimize face " << k << endl; - multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL); - - FaceDescriptor & fd = mesh.GetFaceDescriptor(k); - - PrintMessage (1, "Optimize Surface ", k); - for (int i = 1; i <= mparam.optsteps2d; i++) + //for (int i = 1; i <= mesh.GetNSeg(); i++) + for(auto & seg : segments) { - // (*testout) << "optstep " << i << endl; - if (multithread.terminate) return; + //Segment & seg = mesh.LineSegment(i); + //if (seg.si == k) + { + PointGeomInfo gi0, gi1; + gi0.trignum = gi1.trignum = k; + gi0.u = seg.epgeominfo[0].u; + gi0.v = seg.epgeominfo[0].v; + gi1.u = seg.epgeominfo[1].u; + gi1.v = seg.epgeominfo[1].v; - { - MeshOptimize2d meshopt(mesh); - meshopt.SetFaceIndex (k); - meshopt.SetImproveEdges (0); - meshopt.SetMetricWeight (mparam.elsizeweight); - meshopt.SetWriteStatus (0); - meshopt.EdgeSwapping (i > mparam.optsteps2d/2); - } + int locpnum[2] = {0, 0}; - if (multithread.terminate) return; - { - MeshOptimize2d meshopt(mesh); - meshopt.SetFaceIndex (k); - meshopt.SetImproveEdges (0); - meshopt.SetMetricWeight (mparam.elsizeweight); - meshopt.SetWriteStatus (0); - meshopt.ImproveMesh (mparam); - } + for (int j = 0; j < 2; j++) + { + PointGeomInfo gi = (j == 0) ? gi0 : gi1; - { - MeshOptimize2d meshopt(mesh); - meshopt.SetFaceIndex (k); - meshopt.SetImproveEdges (0); - meshopt.SetMetricWeight (mparam.elsizeweight); - meshopt.SetWriteStatus (0); - meshopt.CombineImprove (); - } + /* + int l; + for (l = 0; l < gis.Size() && locpnum[j] == 0; l++) + { + double dist = sqr (gis[l].u-gi.u)+sqr(gis[l].v-gi.v); - if (multithread.terminate) return; - { - MeshOptimize2d meshopt(mesh); - meshopt.SetFaceIndex (k); - meshopt.SetImproveEdges (0); - meshopt.SetMetricWeight (mparam.elsizeweight); - meshopt.SetWriteStatus (0); - meshopt.ImproveMesh (mparam); - } + if (dist < 1e-10) + locpnum[j] = l+1; + } + + if (locpnum[j] == 0) + { + PointIndex pi = seg[j]; + meshing.AddPoint (mesh.Point(pi), pi); + + gis.SetSize (gis.Size()+1); + gis[l] = gi; + locpnum[j] = l+1; + } + */ + for (int l = 0; l < gis.Size(); l++) + { + double dist = sqr (gis[l].u-gi.u)+sqr(gis[l].v-gi.v); + if (dist < 1e-10) + { + locpnum[j] = l+1; + break; + } + } + + if (locpnum[j] == 0) + { + PointIndex pi = seg[j]; + meshing.AddPoint (mesh.Point(pi), pi); + + gis.Append (gi); + locpnum[j] = gis.Size(); + } + } + + meshing.AddBoundaryElement (locpnum[0], locpnum[1], gi0, gi1); + } } - } - mesh.CalcSurfacesOfNode(); - mesh.Compress(); - timer_opt2d.Stop(); + // Philippose - 15/01/2009 + double maxh = min2(geom.face_maxh[k-1], OCCGeometry::global_shape_properties[TopoDS::Face(geom.fmap(k)).TShape()].maxh); + //double maxh = mparam.maxh; + // int noldpoints = mesh->GetNP(); + int noldsurfel = mesh.GetNSE(); - multithread.task = savetask; + static Timer tsurfprop("surfprop"); + tsurfprop.Start(); + GProp_GProps sprops; + BRepGProp::SurfaceProperties(TopoDS::Face(geom.fmap(k)),sprops); + tsurfprop.Stop(); + meshing.SetMaxArea(2.*sprops.Mass()); + + MESHING2_RESULT res; + + // TODO: check overlap not correctly working here + MeshingParameters mparam_without_overlap = mparam; + mparam_without_overlap.checkoverlap = false; + + try { + static Timer t("GenerateMesh"); RegionTimer reg(t); + res = meshing.GenerateMesh (mesh, mparam_without_overlap, maxh, k); + } + + catch (SingularMatrixException) + { + // (*myerr) << "Singular Matrix" << endl; + res = MESHING2_GIVEUP; + } + + catch (UVBoundsException) + { + // (*myerr) << "UV bounds exceeded" << endl; + res = MESHING2_GIVEUP; + } + + projecttype = PARAMETERSPACE; + static Timer t1("rest of loop"); RegionTimer reg1(t1); + + bool meshing_failed = res != MESHING2_OK; + if(meshing_failed && delete_on_failure) + { + for (SurfaceElementIndex sei = noldsurfel; sei < mesh.GetNSE(); sei++) + mesh.Delete(sei); + + mesh.Compress(); + } + + for (SurfaceElementIndex sei = oldnf; sei < mesh.GetNSE(); sei++) + mesh[sei].SetIndex (k); + + auto n_illegal_trigs = mesh.FindIllegalTrigs(); + PrintMessage (3, n_illegal_trigs, " illegal triangles"); + return meshing_failed; } - void OCCSetLocalMeshSize(const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam, const OCCParameters& occparam) { diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index d5a85e97..2ed2c002 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -1,52 +1,69 @@ #ifdef OCCGEOMETRY -#include -#include -#include #include -#include "ShapeAnalysis_ShapeTolerance.hxx" -#include "ShapeAnalysis_ShapeContents.hxx" -#include "ShapeAnalysis_CheckSmallFace.hxx" -#include "ShapeAnalysis_DataMapOfShapeListOfReal.hxx" -#include "ShapeAnalysis_Surface.hxx" +#include -#include "BRepCheck_Analyzer.hxx" -#include "BRepLib.hxx" -#include "ShapeBuild_ReShape.hxx" -#include "ShapeFix.hxx" -#include "ShapeFix_FixSmallFace.hxx" -#include "Partition_Spliter.hxx" -#include "BRepAlgoAPI_Fuse.hxx" -#include "Interface_InterfaceModel.hxx" +#include +#include -#include "XSControl_WorkSession.hxx" -#include "XSControl_TransferReader.hxx" -#include "StepRepr_RepresentationItem.hxx" -#include "StepBasic_ProductDefinitionRelationship.hxx" -#include "Transfer_TransientProcess.hxx" -#include "TransferBRep.hxx" -#ifndef _Standard_Version_HeaderFile -#include -#endif +#include "occ_vertex.hpp" +#include "occ_edge.hpp" +#include "occ_face.hpp" +#include "occ_solid.hpp" +#include "occgeom.hpp" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include -#include -#include -#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include +#include +#include +#include +#include +#include +#include +#include +#include #include -#include -#include -#include +#include +#include #if OCC_VERSION_HEX < 0x070000 // pass #elif OCC_VERSION_HEX < 0x070200 - #include "StlTransfer.hxx" - #include "TopoDS_Iterator.hxx" + #include + #include #else - #include "TopoDS_Iterator.hxx" + #include #endif namespace netgen @@ -57,6 +74,67 @@ namespace netgen std::map OCCGeometry::global_shape_properties; std::map> OCCGeometry::identifications; + TopoDS_Shape ListOfShapes::Max(gp_Vec dir) + { + double maxval = -1e99; + TopoDS_Shape maxshape; + for (auto shape : *this) + { + GProp_GProps props; + gp_Pnt center; + + switch (shape.ShapeType()) + { + case TopAbs_VERTEX: + center = BRep_Tool::Pnt (TopoDS::Vertex(shape)); break; + case TopAbs_FACE: + BRepGProp::SurfaceProperties (shape, props); + center = props.CentreOfMass(); + break; + default: + BRepGProp::LinearProperties(shape, props); + center = props.CentreOfMass(); + } + + double val = center.X()*dir.X() + center.Y()*dir.Y() + center.Z() * dir.Z(); + if (val > maxval) + { + maxval = val; + maxshape = shape; + } + } + return maxshape; + } + TopoDS_Shape ListOfShapes::Nearest(gp_Pnt pnt) + { + double mindist = 1e99; + TopoDS_Shape nearestshape; + auto vertex = BRepBuilderAPI_MakeVertex (pnt).Vertex(); + + for (auto shape : *this) + { + double dist = BRepExtrema_DistShapeShape(shape, vertex).Value(); + if (dist < mindist) + { + nearestshape = shape; + mindist = dist; + } + } + return nearestshape; + } + + ListOfShapes ListOfShapes::SubShapes(TopAbs_ShapeEnum type) const + { + std::set unique_shapes; + for(const auto& shape : *this) + for(TopExp_Explorer e(shape, type); e.More(); e.Next()) + unique_shapes.insert(e.Current()); + ListOfShapes sub; + for(const auto& shape : unique_shapes) + sub.push_back(shape); + return sub; + } + OCCGeometry::OCCGeometry(const TopoDS_Shape& _shape, int aoccdim, bool copy) { if(copy) @@ -64,14 +142,14 @@ namespace netgen auto filename = GetTempFilename(); step_utils::WriteSTEP(_shape, filename.c_str()); LoadOCCInto(this, filename.c_str()); - occdim = aoccdim; + dimension = aoccdim; std::remove(filename.c_str()); } else { shape = _shape; changed = 1; - occdim = aoccdim; + dimension = aoccdim; BuildFMap(); CalcBoundingBox(); PrintContents (this); @@ -117,23 +195,23 @@ namespace netgen OCCSetLocalMeshSize(*this, mesh, mparam, occparam); } - void OCCGeometry :: FindEdges(Mesh& mesh, - const MeshingParameters& mparam) const + bool OCCGeometry :: MeshFace(Mesh& mesh, + const MeshingParameters& mparam, int nr, FlatArray glob2loc) const { - OCCFindEdges(*this, mesh, mparam); - } + bool failed = OCCMeshFace(*this, mesh, glob2loc, mparam, nr, PLANESPACE, true); + if(failed) + failed = OCCMeshFace(*this, mesh, glob2loc, mparam, nr, PARAMETERSPACE, false); - void OCCGeometry :: MeshSurface(Mesh& mesh, - const MeshingParameters& mparam) const - { - OCCMeshSurface(*this, mesh, mparam); - } - - void OCCGeometry :: FinalizeMesh(Mesh& mesh) const - { - for (int i = 0; i < mesh.GetNDomains(); i++) - if (auto name = sprops[i]->name) - mesh.SetMaterial (i+1, *name); + if(failed) + { + facemeshstatus[nr] = -1; + PrintError ("Problem in Surface mesh generation"); + } + else + { + facemeshstatus[nr] = 1; + } + return failed; } void OCCGeometry :: PrintNrShapes () @@ -1029,31 +1107,97 @@ namespace netgen fsingular = esingular = vsingular = false; - - sprops.SetSize(somap.Extent()); - for (TopExp_Explorer e(shape, TopAbs_SOLID); e.More(); e.Next()) + // Add shapes + for(auto v : GetVertices(shape)) { - auto s = e.Current(); - sprops[somap.FindIndex(s)-1] = &global_shape_properties[s.TShape()]; + auto tshape = v.TShape(); + if(vertex_map.count(tshape)!=0) + continue; + vertex_map[tshape] = vertices.Size(); + auto occ_vertex = make_unique(TopoDS::Vertex(v)); + if(global_shape_properties.count(tshape)>0) + occ_vertex->properties = global_shape_properties[tshape]; + vertices.Append(std::move(occ_vertex)); } - fprops.SetSize(fmap.Extent()); - for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next()) + for(auto e : GetEdges(shape)) { - auto s = e.Current(); - fprops[fmap.FindIndex(s)-1] = &global_shape_properties[s.TShape()]; + auto tshape = e.TShape(); + auto edge = TopoDS::Edge(e); + if(edge_map.count(tshape)!=0) + continue; + if(BRep_Tool::Degenerated(edge)) + continue; + edge_map[tshape] = edges.Size(); + auto occ_edge = make_unique(edge); + occ_edge->properties = global_shape_properties[tshape]; + edges.Append(std::move(occ_edge)); } - eprops.SetSize(emap.Extent()); - /* - for (TopExp_Explorer e(shape, TopAbs_EDGE); e.More(); e.Next()) + for(auto f : GetFaces(shape)) { - auto s = e.Current(); - eprops[emap.FindIndex(s)-1] = &global_shape_properties[s.TShape()]; + auto tshape = f.TShape(); + if(face_map.count(tshape)==0) + { + + auto k = faces.Size(); + face_map[tshape] = k; + auto occ_face = make_unique(f); + if(global_shape_properties.count(tshape)>0) + occ_face->properties = global_shape_properties[tshape]; + faces.Append(std::move(occ_face)); + } } - */ - for (auto [nr,s] : Enumerate(emap)) - eprops[nr-1] = &global_shape_properties[s.TShape()]; + + for(auto s : GetSolids(shape)) + { + auto tshape = s.TShape(); + int k; + if(solid_map.count(tshape)==0) + { + k = solids.Size(); + solid_map[tshape] = k; + auto occ_solid = make_unique(s); + if(global_shape_properties.count(tshape)>0) + occ_solid->properties = global_shape_properties[tshape]; + solids.Append(std::move(occ_solid)); + } + + for(auto f : GetFaces(s)) + { + auto face_nr = face_map[f.TShape()]; + auto & face = faces[face_nr]; + if(face->domin==-1) + face->domin = k; + else + face->domout = k; + } + } + + // Add identifications + auto add_identifications = [&](auto & shapes, auto & shape_map) + { + for(auto &[tshape, nr] : shape_map) + if(identifications.count(tshape)) + for(auto & ident : identifications[tshape]) + { + ShapeIdentification si{ + shapes[shape_map[ident.from]].get(), + shapes[shape_map[ident.to]].get(), + ident.trafo, + ident.type, + ident.name + }; + shapes[nr]->identifications.Append(si); + } + }; + add_identifications( vertices, vertex_map ); + add_identifications( edges, edge_map ); + add_identifications( faces, face_map ); + + ProcessIdentifications(); + + bounding_box = ::netgen::GetBoundingBox( shape ); } @@ -1156,286 +1300,11 @@ namespace netgen void OCCGeometry :: CalcBoundingBox () { - Bnd_Box bb; -#if OCC_VERSION_HEX < 0x070000 - BRepBndLib::Add (shape, bb); -#else - BRepBndLib::Add ((const TopoDS_Shape) shape, bb,(Standard_Boolean)true); -#endif - - double x1,y1,z1,x2,y2,z2; - bb.Get (x1,y1,z1,x2,y2,z2); - Point<3> p1 = Point<3> (x1,y1,z1); - Point<3> p2 = Point<3> (x2,y2,z2); - - (*testout) << "Bounding Box = [" << p1 << " - " << p2 << "]" << endl; - boundingbox = Box<3> (p1,p2); + boundingbox = ::netgen::GetBoundingBox(shape); + (*testout) << "Bounding Box = [" << boundingbox.PMin() << " - " << boundingbox.PMax() << "]" << endl; SetCenter(); } - PointGeomInfo OCCGeometry :: ProjectPoint(int surfi, Point<3> & p) const - { - static int cnt = 0; - if (++cnt % 1000 == 0) cout << "Project cnt = " << cnt << endl; - - gp_Pnt pnt(p(0), p(1), p(2)); - - double u,v; - Handle( Geom_Surface ) thesurf = BRep_Tool::Surface(TopoDS::Face(fmap(surfi))); - Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( thesurf ); - gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfi)) ) ); - 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 - { - double u = gi.u; - double v = gi.v; - - Point<3> hp = p; - if (FastProject (surfind, hp, u, v)) - { - p = hp; - return 1; - } - ProjectPoint (surfind, p); - return CalcPointGeomInfo (surfind, gi, p); - } - - void OCCGeometry :: ProjectPointEdge(int surfind, INDEX surfind2, - Point<3> & p, EdgePointGeomInfo* gi) const - { - TopExp_Explorer exp0, exp1; - bool done = false; - Handle(Geom_Curve) c; - - for (exp0.Init(fmap(surfind), TopAbs_EDGE); !done && exp0.More(); exp0.Next()) - for (exp1.Init(fmap(surfind2), TopAbs_EDGE); !done && exp1.More(); exp1.Next()) - { - if (TopoDS::Edge(exp0.Current()).IsSame(TopoDS::Edge(exp1.Current()))) - { - done = true; - double s0, s1; - c = BRep_Tool::Curve(TopoDS::Edge(exp0.Current()), s0, s1); - } - } - - gp_Pnt pnt(p(0), p(1), p(2)); - GeomAPI_ProjectPointOnCurve proj(pnt, c); - pnt = proj.NearestPoint(); - p(0) = pnt.X(); - p(1) = pnt.Y(); - p(2) = pnt.Z(); - - } - - bool OCCGeometry :: FastProject (int surfi, Point<3> & ap, double& u, double& v) const - { - gp_Pnt p(ap(0), ap(1), ap(2)); - - Handle(Geom_Surface) surface = BRep_Tool::Surface(TopoDS::Face(fmap(surfi))); - - gp_Pnt x = surface->Value (u,v); - - if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true; - - gp_Vec du, dv; - - surface->D1(u,v,x,du,dv); - - int count = 0; - - gp_Pnt xold; - gp_Vec n; - double det, lambda, mu; - - do { - count++; - - n = du^dv; - - det = Det3 (n.X(), du.X(), dv.X(), - n.Y(), du.Y(), dv.Y(), - n.Z(), du.Z(), dv.Z()); - - if (det < 1e-15) return false; - - lambda = Det3 (n.X(), p.X()-x.X(), dv.X(), - n.Y(), p.Y()-x.Y(), dv.Y(), - n.Z(), p.Z()-x.Z(), dv.Z())/det; - - mu = Det3 (n.X(), du.X(), p.X()-x.X(), - n.Y(), du.Y(), p.Y()-x.Y(), - n.Z(), du.Z(), p.Z()-x.Z())/det; - - u += lambda; - v += mu; - - xold = x; - surface->D1(u,v,x,du,dv); - - } while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50); - - // (*testout) << "FastProject count: " << count << endl; - - if (count == 50) return false; - - ap = Point<3> (x.X(), x.Y(), x.Z()); - - return true; - } - - Vec<3> OCCGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* geominfo) const - { - if(geominfo) - { - gp_Pnt pnt; - gp_Vec du, dv; - - Handle(Geom_Surface) occface; - occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind))); - - 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(); - - if (fmap(surfind).Orientation() == TopAbs_REVERSED) n *= -1; - return n; - } - Standard_Real u,v; - - gp_Pnt pnt(p(0), p(1), p(2)); - - Handle(Geom_Surface) occface; - occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind))); - - /* - GeomAPI_ProjectPointOnSurf proj(pnt, occface); - - if (proj.NbPoints() < 1) - { - cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!" - << endl; - cout << p << endl; - return; - } - - proj.LowerDistanceParameters (u, v); - */ - - Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface ); - gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfind)) ) ); - suval.Coord( u, v); - pnt = occface->Value( u, v ); - - gp_Vec du, dv; - occface->D1(u,v,pnt,du,dv); - - /* - if (!occface->IsCNu (1) || !occface->IsCNv (1)) - (*testout) << "SurfOpt: Differentiation FAIL" << endl; - */ - - auto n = Cross (Vec3d(du.X(), du.Y(), du.Z()), - Vec3d(dv.X(), dv.Y(), dv.Z())); - n.Normalize(); - - if (fmap(surfind).Orientation() == TopAbs_REVERSED) n *= -1; - return n; - } - - bool OCCGeometry :: CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p) const - { - Standard_Real u,v; - - gp_Pnt pnt(p(0), p(1), p(2)); - - Handle(Geom_Surface) occface; - occface = BRep_Tool::Surface(TopoDS::Face(fmap(surfind))); - - /* - GeomAPI_ProjectPointOnSurf proj(pnt, occface); - - if (proj.NbPoints() < 1) - { - cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!" - << endl; - cout << p << endl; - return 0; - } - - proj.LowerDistanceParameters (u, v); - */ - - Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface ); - gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfind)) ) ); - suval.Coord( u, v); - //pnt = occface->Value( u, v ); - - - gi.u = u; - gi.v = v; - return true; - } - - void OCCGeometry :: PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint, - int surfi, - const PointGeomInfo & gi1, - const PointGeomInfo & gi2, - Point<3> & newp, PointGeomInfo & newgi) const - { - Point<3> hnewp; - hnewp = p1+secpoint*(p2-p1); - - if (surfi > 0) - { - double u = gi1.u+secpoint*(gi2.u-gi1.u); - double v = gi1.v+secpoint*(gi2.v-gi1.v); - - auto savept = hnewp; - if (!FastProject(surfi, hnewp, u, v) || Dist(hnewp, savept) > Dist(p1,p2)) - { - // cout << "Fast projection to surface fails! Using OCC projection" << endl; - hnewp = savept; - ProjectPoint(surfi, hnewp); - } - newgi.trignum = 1; - newgi.u = u; - newgi.v = v; - } - newp = hnewp; - } - - - void OCCGeometry :: PointBetweenEdge(const Point<3> & p1, - const Point<3> & p2, double secpoint, - int surfi1, int surfi2, - const EdgePointGeomInfo & ap1, - const EdgePointGeomInfo & ap2, - Point<3> & newp, EdgePointGeomInfo & newgi) const - { - double s0, s1; - - newp = p1+secpoint*(p2-p1); - if(ap1.edgenr > emap.Size() || ap1.edgenr == 0) - return; - gp_Pnt pnt(newp(0), newp(1), newp(2)); - GeomAPI_ProjectPointOnCurve proj(pnt, BRep_Tool::Curve(TopoDS::Edge(emap(ap1.edgenr)), s0, s1)); - pnt = proj.NearestPoint(); - newp = Point<3> (pnt.X(), pnt.Y(), pnt.Z()); - newgi = ap1; - }; - // void OCCGeometry :: WriteOCC_STL(char * filename) // { @@ -1699,26 +1568,18 @@ namespace netgen auto occ_hash = key.HashCode(1<<31UL); return std::hash()(occ_hash); }; - std::unordered_map shape_map(10, my_hash); - Array shape_list; - std::map tshape_map; Array tshape_list; - ar & occdim; + ar & dimension; for (auto typ : { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE }) for (TopExp_Explorer e(shape, typ); e.More(); e.Next()) { auto ds = e.Current(); auto ts = ds.TShape(); - if(shape_map.count(ds)==0) - { - shape_map[ds] = shape_list.Size(); - shape_list.Append(ds); - } if(tshape_map.count(ts)==0) { - tshape_map[ts] = shape_list.Size(); + tshape_map[ts] = tshape_list.Size(); tshape_list.Append(ts); } } @@ -1741,12 +1602,18 @@ namespace netgen for(auto i : Range(n_idents)) { auto & id = idents[i]; - int shape_id; + int id_from, id_to; if(ar.Output()) - shape_id = shape_map[id.other]; - ar & shape_id & id.trafo & id.inverse & id.name; + { + id_from = tshape_map[id.from]; + id_to = tshape_map[id.to]; + } + ar & id_from & id_to & id.trafo & id.name; if(ar.Input()) - id.other = shape_list[shape_id]; + { + id.from = tshape_list[id_from]; + id.to = tshape_list[id_to]; + } } } } @@ -2053,6 +1920,97 @@ namespace netgen return false; } + Point<3> GetCenter(const TopoDS_Shape & shape) + { + GProp_GProps props; + BRepGProp::LinearProperties(shape, props); + return occ2ng( props.CentreOfMass() ); + } + + void OCCGeometry :: IdentifyEdges(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type) + { + auto cme = GetCenter(me); + auto cyou = GetCenter(you); + Transformation<3> trafo{cyou-cme}; + identifications[me.TShape()].push_back( {me.TShape(), you.TShape(), Transformation<3>(cyou - cme), name, type} ); + + auto vme = GetVertices(me); + auto vyou = GetVertices(you); + Point<3> pme0 = trafo(occ2ng(vme[0])); + Point<3> pme1 = trafo(occ2ng(vme[1])); + Point<3> pyou = occ2ng(vyou[0]); + + bool do_swap = Dist(pme1, pyou) < Dist(pme0, pyou); + + for(auto i : Range(2)) + identifications[vme[i].TShape()].push_back( {vme[i].TShape(), vyou[do_swap ? 1-i : i].TShape(), trafo, name, type} ); + } + + bool IsMappedShape(const Transformation<3> & trafo, const TopoDS_Shape & me, const TopoDS_Shape & you) + { + if(me.ShapeType() != you.ShapeType()) return false; + + Bnd_Box bbox; + BRepBndLib::Add(me, bbox); + BRepBndLib::Add(you, bbox); + BoxTree<3> tree( occ2ng(bbox.CornerMin()), occ2ng(bbox.CornerMax()) ); + + Point<3> c_me = GetCenter(me); + Point<3> c_you = GetCenter(you); + if(tree.GetTolerance() < Dist(trafo(c_me), c_you)) + return false; + + std::map vmap; + + auto verts_me = GetVertices(me); + for (auto i : Range(verts_me.size())) + { + auto s = verts_me[i].TShape(); + if(vmap.count(s)>0) + throw Exception("vertex mapped twice!"); + auto p = trafo(occ2ng(s)); + tree.Insert( p, i ); + vmap[s] = nullptr; + } + + bool all_verts_mapped = true; + for (auto vert : GetVertices(you)) + { + auto s = vert.TShape(); + auto p = occ2ng(s); + bool vert_mapped = false; + tree.GetFirstIntersecting( p, p, [&](size_t i ) { + vmap[verts_me[i].TShape()] = s; + vert_mapped = true; + return true; + }); + if(!vert_mapped) + { + all_verts_mapped = false; + break; + } + } + return all_verts_mapped; + } + + void OCCGeometry :: IdentifyFaces(const TopoDS_Shape & solid, const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type) + { + auto cme = GetCenter(me); + auto cyou = GetCenter(you); + Transformation<3> trafo(cyou-cme); + + identifications[me.TShape()].push_back + (OCCIdentification { me.TShape(), you.TShape(), trafo, name, type }); + + auto edges_me = GetEdges(me); + auto edges_you = GetEdges(you); + + for (auto e_me : edges_me) + for (auto e_you : edges_you) + if(IsMappedShape(trafo, e_me, e_you)) + IdentifyEdges(e_me, e_you, name, type); + } + void OCCParameters :: Print(ostream & ost) const { ost << "OCC Parameters:" << endl @@ -2063,6 +2021,18 @@ namespace netgen DLL_HEADER extern OCCParameters occparam; OCCParameters occparam; + + + + + + + + + + + + // int OCCGeometry :: GenerateMesh (shared_ptr & mesh, MeshingParameters & mparam) // { // return OCCGenerateMesh (*this, mesh, mparam, occparam); @@ -2208,8 +2178,7 @@ namespace netgen for(auto & ident : identifications) { Array items; - items.Append(STEPConstruct::FindEntity(finder, ident.other)); - items.Append(MakeInt(static_cast(ident.inverse))); + // items.Append(STEPConstruct::FindEntity(finder, ident.other)); // TODO! auto & m = ident.trafo.GetMatrix(); for(auto i : Range(9)) items.Append(MakeReal(m(i))); @@ -2239,8 +2208,7 @@ namespace netgen auto id_item = Handle(StepRepr_CompoundRepresentationItem)::DownCast(idents->ItemElementValue(i)); OCCIdentification ident; ident.name = id_item->Name()->ToCString(); - ident.other = TransferBRep::ShapeResult(transProc->Find(id_item->ItemElementValue(1))); - ident.inverse = static_cast(ReadInt(id_item->ItemElementValue(2))); + // ident.other = TransferBRep::ShapeResult(transProc->Find(id_item->ItemElementValue(1))); /TODO! auto & m = ident.trafo.GetMatrix(); for(auto i : Range(9)) diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index 381a3df1..27eae5b7 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -10,99 +10,27 @@ #ifdef OCCGEOMETRY #include +#include "occ_utils.hpp" +#include "occmeshsurf.hpp" -#include -#include "BRep_Tool.hxx" -#include "Geom_Curve.hxx" -#include "Geom2d_Curve.hxx" -#include "Geom_Surface.hxx" -#include "GeomAPI_ProjectPointOnSurf.hxx" -#include "GeomAPI_ProjectPointOnCurve.hxx" -#include "BRepTools.hxx" -#include "TopExp.hxx" -#include "BRepBuilderAPI_MakeVertex.hxx" -#include "BRepBuilderAPI_MakeShell.hxx" -#include "BRepBuilderAPI_MakeSolid.hxx" -#include "BRepOffsetAPI_Sewing.hxx" -#include "BRepLProp_SLProps.hxx" -#include "BRepAdaptor_Surface.hxx" -#include "Poly_Triangulation.hxx" -#include "Poly_Array1OfTriangle.hxx" -#include "TColgp_Array1OfPnt2d.hxx" -#include "Poly_Triangle.hxx" -#include "GProp_GProps.hxx" -#include "BRepGProp.hxx" -#include "gp_Pnt.hxx" -#include "TopoDS.hxx" -#include "TopoDS_Solid.hxx" -#include "TopExp_Explorer.hxx" -#include "TopTools_ListIteratorOfListOfShape.hxx" -#include "TopoDS_Wire.hxx" -#include "BRepTools_WireExplorer.hxx" -#include "TopTools_IndexedMapOfShape.hxx" -#include "BRepLProp_CLProps.hxx" -#include "BRepAdaptor_Curve.hxx" -#include "TopoDS_Shape.hxx" -#include "TopoDS_Face.hxx" -#include "IGESToBRep_Reader.hxx" -#include "Interface_Static.hxx" -#include "GeomAPI_ExtremaCurveCurve.hxx" -#include "Standard_ErrorHandler.hxx" -#include "Standard_Failure.hxx" -#include "ShapeUpgrade_ShellSewing.hxx" -#include "ShapeFix_Shape.hxx" -#include "ShapeFix_Wireframe.hxx" -#include "BRepMesh_IncrementalMesh.hxx" -#include "BRepBndLib.hxx" -#include "Bnd_Box.hxx" -#include "ShapeAnalysis.hxx" -#include "ShapeBuild_ReShape.hxx" -#include "BOPAlgo_Builder.hxx" - -// Philippose - 29/01/2009 -// OpenCascade XDE Support -// Include support for OpenCascade XDE Features -#include "TDocStd_Document.hxx" -#include "Quantity_Color.hxx" -#include "XCAFApp_Application.hxx" -#include "XCAFDoc_ShapeTool.hxx" -#include "XCAFDoc_Color.hxx" -#include "XCAFDoc_ColorTool.hxx" -#include "XCAFDoc_ColorType.hxx" -#include "XCAFDoc_LayerTool.hxx" -#include "XCAFDoc_DimTolTool.hxx" -#include "XCAFDoc_MaterialTool.hxx" -#include "XCAFDoc_DocumentTool.hxx" -#include "TDF_Label.hxx" -#include "TDF_LabelSequence.hxx" -#include "STEPCAFControl_Reader.hxx" -#include "STEPCAFControl_Writer.hxx" -#include "IGESCAFControl_Reader.hxx" -#include "IGESCAFControl_Writer.hxx" - -#include "IGESControl_Reader.hxx" -#include "STEPControl_Reader.hxx" -#include "IGESControl_Writer.hxx" -#include "STEPControl_Writer.hxx" - -#include -#include -#include +#include +#include #include - -#include "StlAPI_Writer.hxx" -#include "STEPControl_StepModelType.hxx" +#include +#include +#include +#include +#include +#include +#include +#include #if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4 #define OCC_HAVE_HISTORY #endif - - - namespace netgen { -#include "occmeshsurf.hpp" // extern DLL_HEADER MeshingParameters mparam; @@ -117,28 +45,6 @@ namespace netgen #define OCCGEOMETRYVISUALIZATIONHALFCHANGE 2 // Redraw - inline Point<3> occ2ng (const gp_Pnt & p) - { - return Point<3> (p.X(), p.Y(), p.Z()); - } - - inline Point<2> occ2ng (const gp_Pnt2d & p) - { - return Point<2> (p.X(), p.Y()); - } - - inline Vec<3> occ2ng (const gp_Vec & v) - { - return Vec<3> (v.X(), v.Y(), v.Z()); - } - - inline gp_Pnt ng2occ (const Point<3> & p) - { - return gp_Pnt(p(0), p(1), p(2)); - } - - - class EntityVisualizationCode { int code; @@ -219,108 +125,20 @@ namespace netgen }; - class ShapeProperties - { - public: - optional name; - optional> col; - double maxh = 1e99; - double hpref = 0; // number of hp refinement levels (will be multiplied by factor later) - void Merge(const ShapeProperties & prop2) - { - if (prop2.name) name = prop2.name; - if (prop2.col) col = prop2.col; - maxh = min2(maxh, prop2.maxh); - } - - void DoArchive(Archive& ar) - { - ar & name & col & maxh & hpref; - } - }; - - - class OCCIdentification - { - public: - TopoDS_Shape other; - Transformation<3> trafo; - bool inverse; - string name; - }; - - - class MyExplorer - { - class Iterator - { - TopExp_Explorer exp; - public: - Iterator (TopoDS_Shape ashape, TopAbs_ShapeEnum atoFind, TopAbs_ShapeEnum atoAvoid) - : exp(ashape, atoFind, atoAvoid) { } - auto operator*() { return exp.Current(); } - Iterator & operator++() { exp.Next(); return *this; } - bool operator!= (nullptr_t nu) { return exp.More(); } - }; - - public: - TopoDS_Shape shape; - TopAbs_ShapeEnum toFind; - TopAbs_ShapeEnum toAvoid; - MyExplorer (TopoDS_Shape ashape, TopAbs_ShapeEnum atoFind, TopAbs_ShapeEnum atoAvoid = TopAbs_SHAPE) - : shape(ashape), toFind(atoFind), toAvoid(atoAvoid) { ; } - Iterator begin() { return Iterator(shape, toFind, toAvoid); } - auto end() { return nullptr; } - }; - - inline auto Explore (TopoDS_Shape shape, TopAbs_ShapeEnum toFind, TopAbs_ShapeEnum toAvoid = TopAbs_SHAPE) - { - return MyExplorer (shape, toFind, toAvoid); - } - - - class IndexMapIterator - { - class Iterator - { - const TopTools_IndexedMapOfShape & indmap; - int i; - public: - Iterator (const TopTools_IndexedMapOfShape & aindmap, int ai) - : indmap(aindmap), i(ai) { ; } - auto operator*() { return tuple(i, indmap(i)); } - Iterator & operator++() { i++; return *this; } - bool operator!= (const Iterator & i2) { return i != i2.i; } - }; - - public: - const TopTools_IndexedMapOfShape & indmap; - IndexMapIterator (const TopTools_IndexedMapOfShape & aindmap) : indmap(aindmap) { } - Iterator begin() { return Iterator(indmap, 1); } - Iterator end() { return Iterator(indmap, indmap.Extent()+1); } - }; - - inline auto Enumerate (const TopTools_IndexedMapOfShape & indmap) - { - return IndexMapIterator(indmap); - } - - class DLL_HEADER OCCGeometry : public NetgenGeometry { Point<3> center; OCCParameters occparam; public: - static std::map global_shape_properties; - static std::map> identifications; + static std::map global_shape_properties; + static std::map> identifications; TopoDS_Shape shape; - TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap; + TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap; // legacy maps NgArray fsingular, esingular, vsingular; Box<3> boundingbox; - // should we use 1-based arrays (JS->MH) ? - Array fprops, eprops, sprops; // pointers to the gobal property map + std::map edge_map, vertex_map, face_map, solid_map; mutable int changed; mutable NgArray facemeshstatus; @@ -350,8 +168,6 @@ namespace netgen bool makesolids; bool splitpartitions; - int occdim = 3; // meshing is always done 3D, changed to 2D later of occdim=2 - OCCGeometry() { somap.Clear(); @@ -372,36 +188,15 @@ namespace netgen void Analyse(Mesh& mesh, const MeshingParameters& mparam) const override; - void FindEdges(Mesh& mesh, - const MeshingParameters& mparam) const override; - void MeshSurface(Mesh& mesh, - const MeshingParameters& mparam) const override; + bool MeshFace(Mesh& mesh, const MeshingParameters& mparam, + int nr, FlatArray glob2loc) const override; + // void OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const override {} - void FinalizeMesh(Mesh& mesh) const override; - void Save (string filename) const override; void SaveToMeshFile (ostream & /* ost */) const override; void DoArchive(Archive& ar) override; - PointGeomInfo ProjectPoint(int surfind, 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 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, - int surfi1, int surfi2, - const EdgePointGeomInfo & ap1, - const EdgePointGeomInfo & ap2, - Point<3> & newp, EdgePointGeomInfo & newgi) const override; - void PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint, - int surfi, - const PointGeomInfo & gi1, - const PointGeomInfo & gi2, - Point<3> & newp, PointGeomInfo & newgi) const override; - void BuildFMap(); auto GetShape() const { return shape; } @@ -546,9 +341,11 @@ namespace netgen bool ErrorInSurfaceMeshing (); // void WriteOCC_STL(char * filename); + static void IdentifyEdges(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type); + static void IdentifyFaces(const TopoDS_Shape & solid,const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type); private: - bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const; + //bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const; }; @@ -564,11 +361,8 @@ namespace netgen DLL_HEADER extern void OCCSetLocalMeshSize(const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam, const OCCParameters& occparam); - DLL_HEADER extern void OCCMeshSurface (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam); - - DLL_HEADER extern void OCCOptimizeSurface (OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam); - - DLL_HEADER extern void OCCFindEdges (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam); + DLL_HEADER extern bool OCCMeshFace (const OCCGeometry & geom, Mesh & mesh, FlatArray glob2loc, + const MeshingParameters & mparam, int nr, int projecttype, bool delete_on_failure); namespace step_utils diff --git a/libsrc/occ/occmeshsurf.cpp b/libsrc/occ/occmeshsurf.cpp index 8f815b5c..8cadda3e 100644 --- a/libsrc/occ/occmeshsurf.cpp +++ b/libsrc/occ/occmeshsurf.cpp @@ -2,16 +2,16 @@ #include -#include #include +#include "occgeom.hpp" + #include #include +#include "occmeshsurf.hpp" namespace netgen { -#include "occmeshsurf.hpp" - bool glob_testout(false); diff --git a/libsrc/occ/occmeshsurf.hpp b/libsrc/occ/occmeshsurf.hpp index 2488a18f..b4c4b06e 100644 --- a/libsrc/occ/occmeshsurf.hpp +++ b/libsrc/occ/occmeshsurf.hpp @@ -6,9 +6,15 @@ #include "occgeom.hpp" #include "mydefs.hpp" +#include +#include +#include + #define PARAMETERSPACE -1 #define PLANESPACE 1 +namespace netgen +{ class OCCGeometry; class SingularMatrixException @@ -144,8 +150,8 @@ protected: class OCCGeometry; -#endif - - +} // namespace netgen + +#endif #endif diff --git a/libsrc/occ/occpkg.cpp b/libsrc/occ/occpkg.cpp index 21be4007..21cf0822 100644 --- a/libsrc/occ/occpkg.cpp +++ b/libsrc/occ/occpkg.cpp @@ -15,6 +15,9 @@ #include "vsocc.hpp" +#include +#include + // __declspec(dllimport) void AutoColourBcProps(Mesh & mesh, const char *bccolourfile); // __declspec(dllimport) void GetFaceColours(Mesh & mesh, NgArray & face_colours); // __declspec(dllimport) bool ColourMatch(Vec3d col1, Vec3d col2, double eps = 2.5e-05); diff --git a/libsrc/occ/python_occ.cpp b/libsrc/occ/python_occ.cpp index 4c1e1330..957d3bfa 100644 --- a/libsrc/occ/python_occ.cpp +++ b/libsrc/occ/python_occ.cpp @@ -1,58 +1,25 @@ #ifdef NG_PYTHON #ifdef OCCGEOMETRY -#include <../general/ngpython.hpp> -#include -#include "../meshing/python_mesh.hpp" #include +#include +#include +#include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -// #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include -#include +#include "occgeom.hpp" +#include +#include #include - -#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4 -#define OCC_HAVE_DUMP_JSON -#endif +#include +#include +#include +#include +#include +#include +#include using namespace netgen; @@ -297,8 +264,6 @@ DLL_HEADER void ExportNgOCC(py::module &m) auto result = geo->GenerateMesh(mesh, mp); if(result != 0) throw Exception("Meshing failed!"); - if (geo->occdim==2) - mesh->SetDimension(2); SetGlobalMesh(mesh); ng_geometry = geo; return mesh; diff --git a/libsrc/occ/python_occ.hpp b/libsrc/occ/python_occ.hpp deleted file mode 100644 index 605875a0..00000000 --- a/libsrc/occ/python_occ.hpp +++ /dev/null @@ -1,68 +0,0 @@ -class DirectionalInterval -{ -public: - gp_Vec dir; - double minval = -1e99; - double maxval = 1e99; - bool openmin = false, openmax = false; - - DirectionalInterval (gp_Vec adir) : dir(adir) { ; } - DirectionalInterval (const DirectionalInterval & i2) - : dir(i2.dir), minval(i2.minval), maxval(i2.maxval) { ; } - - DirectionalInterval operator< (double val) const - { - DirectionalInterval i2 = *this; - i2.maxval = val; - return i2; - } - - DirectionalInterval operator> (double val) const - { - DirectionalInterval i2 = *this; - i2.minval = val; - return i2; - } - - - DirectionalInterval Intersect (const DirectionalInterval & i2) - { - DirectionalInterval res = *this; - res.minval = max(res.minval, i2.minval); - res.maxval = min(res.maxval, i2.maxval); - return res; - } - - bool Contains (gp_Pnt p, double eps = 1e-8) - { - // cout << "Contains point " << p.X() << "," << p.Y() << "," << p.Z() << " ? " << endl; - double val = dir.X()*p.X() + dir.Y()*p.Y() + dir.Z() * p.Z(); - // cout << "minval = " << minval << ", val = " << val << " maxval = " << maxval << endl; - if (openmin) { - if (val < minval+eps) return false; - } else { - if (val < minval-eps) return false; - } - if (openmax) { - if (val > maxval-eps) return false; - } else { - if (val > maxval+eps) return false; - } - return true; - } -}; - - -inline gp_Pnt Center (TopoDS_Shape shape) -{ - GProp_GProps props; - switch (shape.ShapeType()) - { - case TopAbs_FACE: - BRepGProp::SurfaceProperties (shape, props); break; - default: - BRepGProp::LinearProperties(shape, props); - } - return props.CentreOfMass(); -} - diff --git a/libsrc/occ/python_occ_basic.cpp b/libsrc/occ/python_occ_basic.cpp index 2d955770..f658998c 100644 --- a/libsrc/occ/python_occ_basic.cpp +++ b/libsrc/occ/python_occ_basic.cpp @@ -1,59 +1,21 @@ #ifdef NG_PYTHON #ifdef OCCGEOMETRY -#include <../general/ngpython.hpp> +#include #include -#include "../meshing/python_mesh.hpp" - +#include #include -#include +#include "occgeom.hpp" + +#include #include #include #include +#include #include -#include -#include -#include -#include -#include -#include -#include -#include -// #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include -#include - -#include - - -#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4 -#define OCC_HAVE_DUMP_JSON -#endif - - - +using namespace netgen; DLL_HEADER void ExportNgOCCBasic(py::module &m) { diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index 99546609..a05907de 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -3,179 +3,68 @@ #include -#include <../general/ngpython.hpp> +#include #include -#include "../meshing/python_mesh.hpp" - +#include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include "occgeom.hpp" + +#include +#include #include +#include #include -// #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include #include +#include +#include #include #include -#include -#include #include -#include -#include -#include - +#include #include -#include +#include #include - +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include #include #include -#include -#include -#include -#include +#include #include - -#include +#include +#include +#include +#include #include -#include #include - - -#include - -#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4 -#define OCC_HAVE_DUMP_JSON -#endif +#include +#include +#include +#include +#include +#include using namespace netgen; -struct ShapeLess -{ - bool operator() (const TopoDS_Shape& s1, const TopoDS_Shape& s2) const - { - return s1.TShape() < s2.TShape(); - } -}; - -class ListOfShapes : public std::vector -{ -public: - TopoDS_Shape Max(gp_Vec dir) - { - double maxval = -1e99; - TopoDS_Shape maxshape; - for (auto shape : *this) - { - GProp_GProps props; - gp_Pnt center; - - switch (shape.ShapeType()) - { - case TopAbs_VERTEX: - center = BRep_Tool::Pnt (TopoDS::Vertex(shape)); break; - case TopAbs_FACE: - BRepGProp::SurfaceProperties (shape, props); - center = props.CentreOfMass(); - break; - default: - BRepGProp::LinearProperties(shape, props); - center = props.CentreOfMass(); - } - - double val = center.X()*dir.X() + center.Y()*dir.Y() + center.Z() * dir.Z(); - if (val > maxval) - { - maxval = val; - maxshape = shape; - } - } - return maxshape; - } - - TopoDS_Shape Nearest(gp_Pnt pnt) - { - double mindist = 1e99; - TopoDS_Shape nearestshape; - auto vertex = BRepBuilderAPI_MakeVertex (pnt).Vertex(); - - for (auto shape : *this) - { - double dist = BRepExtrema_DistShapeShape(shape, vertex).Value(); - if (dist < mindist) - { - nearestshape = shape; - mindist = dist; - } - } - return nearestshape; - } - - ListOfShapes SubShapes(TopAbs_ShapeEnum type) const - { - std::set unique_shapes; - for(const auto& shape : *this) - for(TopExp_Explorer e(shape, type); e.More(); e.Next()) - unique_shapes.insert(e.Current()); - ListOfShapes sub; - for(const auto& shape : unique_shapes) - sub.push_back(shape); - return sub; - } - ListOfShapes Solids() const - { - return SubShapes(TopAbs_SOLID); - } - ListOfShapes Faces() const - { - return SubShapes(TopAbs_FACE); - } - ListOfShapes Edges() const - { - return SubShapes(TopAbs_EDGE); - } - ListOfShapes Vertices() const - { - return SubShapes(TopAbs_VERTEX); - } - - ListOfShapes operator*(const ListOfShapes& other) const - { - ListOfShapes common; - for(const auto& shape : (*this)) - for(const auto& shape_o : other) - if(shape.IsSame(shape_o)) - common.push_back(shape); - return common; - } -}; - - void ExtractEdgeData( const TopoDS_Edge & edge, int index, std::vector * p, Box<3> & box ) { if (BRep_Tool::Degenerated(edge)) return; @@ -717,36 +606,14 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) return sub; }, py::arg("type"), "returns list of sub-shapes of type 'type'") - .def_property_readonly("solids", [] (const TopoDS_Shape & shape) - { - ListOfShapes solids; - for(TopExp_Explorer e(shape, TopAbs_SOLID); e.More(); e.Next()) - solids.push_back(e.Current()); - return solids; - }, "returns all sub-shapes of type 'SOLID'") - .def_property_readonly("faces", [] (const TopoDS_Shape & shape) - { - ListOfShapes sub; - for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next()) - sub.push_back(e.Current()); - return sub; - }, "returns all sub-shapes of type 'FACE'") - - .def_property_readonly("edges", [] (const TopoDS_Shape & shape) - { - ListOfShapes sub; - for (TopExp_Explorer e(shape, TopAbs_EDGE); e.More(); e.Next()) - sub.push_back(e.Current()); - return sub; - }, "returns all sub-shapes of type 'EDGE'") - - .def_property_readonly("vertices", [] (const TopoDS_Shape & shape) - { - ListOfShapes sub; - for (TopExp_Explorer e(shape, TopAbs_VERTEX); e.More(); e.Next()) - sub.push_back(e.Current()); - return sub; - }, "returns all sub-shapes of type 'VERTEX'") + .def_property_readonly("solids", GetSolids, + "returns all sub-shapes of type 'SOLID'") + .def_property_readonly("faces", GetFaces, + "returns all sub-shapes of type 'FACE'") + .def_property_readonly("edges", GetEdges, + "returns all sub-shapes of type 'EDGE'") + .def_property_readonly("vertices", GetVertices, + "returns all sub-shapes of type 'VERTEX'") .def("Properties", [] (const TopoDS_Shape & shape) { @@ -1127,30 +994,27 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) BRepMesh_IncrementalMesh (shape, deflection, true); }) - .def("Identify", [](const TopoDS_Shape & me, const TopoDS_Shape & you, string name) { + .def("Identify", [](const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE idtype) { // only edges supported, by now - auto me_edge = TopoDS::Edge(me); - auto you_edge = TopoDS::Edge(you); + auto type = me.ShapeType(); + auto tyou = you.ShapeType(); + if(type != tyou) + throw NgException ("Identify: cannot identify different shape types"); - GProp_GProps props; - BRepGProp::LinearProperties(me, props); - gp_Pnt cme = props.CentreOfMass(); - BRepGProp::LinearProperties(you, props); - gp_Pnt cyou = props.CentreOfMass(); - - double s0, s1; - auto curve_me = BRep_Tool::Curve(me_edge, s0, s1); - auto vme = occ2ng(curve_me->Value(s1))-occ2ng(curve_me->Value(s0)); - auto curve_you = BRep_Tool::Curve(you_edge, s0, s1); - auto vyou = occ2ng(curve_you->Value(s1))-occ2ng(curve_you->Value(s0)); - - bool inv = vme*vyou < 0; - OCCGeometry::identifications[me.TShape()].push_back - (OCCIdentification { you, Transformation<3>(occ2ng(cyou) - occ2ng(cme)), inv, name }); - OCCGeometry::identifications[you.TShape()].push_back - (OCCIdentification { me, Transformation<3>(occ2ng(cme) - occ2ng(cyou)), inv, name }); - }, py::arg("other"), py::arg("name"), "Identify shapes for periodic meshing") + switch(type) + { + case TopAbs_VERTEX: + case TopAbs_EDGE: + OCCGeometry::IdentifyEdges(me, you, name, idtype); + break; + default: + throw NgException ("Identify: unsupported shape type"); + break; + } + }, py::arg("other"), py::arg("name"), py::arg("type")=Identifications::PERIODIC, "Identify shapes for periodic meshing") + .def("Identify", OCCGeometry::IdentifyFaces, "Identify faces", + py::arg("from"), py::arg("to"), py::arg("name"), py::arg("type")=Identifications::PERIODIC) .def("Triangulation", [](const TopoDS_Shape & shape) { @@ -1665,6 +1529,18 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) OCCGeometry::global_shape_properties[shape.TShape()].maxh = maxh; } }, "set maxh for all elements of list") + .def_property("hpref", [](ListOfShapes& shapes) + { + throw Exception("Cannot get property of ListOfShapes, get the property from individual shapes!"); + }, + [](ListOfShapes& shapes, double hpref) + { + for(auto& shape : shapes) + { + auto& val = OCCGeometry::global_shape_properties[shape.TShape()].hpref; + val = max2(hpref, val); + } + }, "set hpref for all elements of list") ; diff --git a/libsrc/occ/vsocc.cpp b/libsrc/occ/vsocc.cpp index 3f11f47d..595d0de2 100644 --- a/libsrc/occ/vsocc.cpp +++ b/libsrc/occ/vsocc.cpp @@ -8,19 +8,18 @@ #include -#include "TopoDS_Shape.hxx" -#include "TopoDS_Vertex.hxx" -#include "TopExp_Explorer.hxx" -#include "BRep_Tool.hxx" -#include "TopoDS.hxx" -#include "gp_Pnt.hxx" -#include "Geom_Curve.hxx" -#include "Poly_Triangulation.hxx" -#include "Poly_Array1OfTriangle.hxx" -#include "TColgp_Array1OfPnt2d.hxx" -#include "Poly_Triangle.hxx" -#include "Poly_Polygon3D.hxx" -#include "Poly_PolygonOnTriangulation.hxx" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include diff --git a/nglib/nglib.cpp b/nglib/nglib.cpp index 9eb1e8c3..56e19fe8 100644 --- a/nglib/nglib.cpp +++ b/nglib/nglib.cpp @@ -880,7 +880,7 @@ namespace nglib mp->Transfer_Parameters(); - OCCFindEdges(*occgeom, *me, mparam); + occgeom->FindEdges(*me, mparam); if((me->GetNP()) && (me->GetNFD())) { @@ -926,8 +926,8 @@ namespace nglib perfstepsend = MESHCONST_OPTSURFACE; } - OCCMeshSurface(*occgeom, *me, mparam); - OCCOptimizeSurface(*occgeom, *me, mparam); + occgeom->MeshSurface(*me, mparam); + occgeom->OptimizeSurface(*me, mparam); me->CalcSurfacesOfNode(); diff --git a/python/occ.py b/python/occ.py index c22f3ae0..29f95eb9 100644 --- a/python/occ.py +++ b/python/occ.py @@ -9,6 +9,10 @@ For more detailed documentation consult the OCCT docs, a good starting point is https://dev.opencascade.org/doc/refman/html/class_b_rep_builder_a_p_i___make_shape.html """ +from .config import USE_OCC +if not USE_OCC: + raise ImportError("Netgen was not built with Opencascade support") + from .libngpy._NgOCC import * from .meshing import meshsize diff --git a/tests/dockerfile b/tests/dockerfile index c79b1012..a584f4f8 100644 --- a/tests/dockerfile +++ b/tests/dockerfile @@ -11,12 +11,15 @@ RUN apt-get update && apt-get -y install \ libcgns-dev \ libglu1-mesa-dev \ libhdf5-dev \ + libocct-ocaf-dev \ + libocct-visualization-dev \ libocct-data-exchange-dev \ libocct-draw-dev \ libpython3-dev \ libtbb-dev \ libxi-dev \ libxmu-dev \ + occt-misc \ python3 \ python3-distutils \ python3-numpy \ diff --git a/tests/dockerfile_mpi b/tests/dockerfile_mpi index a25a96b9..fe33fcd8 100644 --- a/tests/dockerfile_mpi +++ b/tests/dockerfile_mpi @@ -1,6 +1,27 @@ FROM ubuntu:20.04 ENV DEBIAN_FRONTEND=noninteractive MAINTAINER Matthias Hochsteger -RUN apt-get update && apt-get -y install python3 libpython3-dev python3-pip libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-numpy python3-tk python3-mpi4py clang-tidy python3-distutils clang libopenmpi-dev openmpi-bin gfortran +RUN apt-get update && apt-get -y install \ + ccache \ + clang \ + clang-tidy \ + cmake \ + g++ \ + gfortran \ + git \ + libglu1-mesa-dev \ + libopenmpi-dev \ + libpython3-dev \ + libxmu-dev \ + openmpi-bin \ + python3 \ + python3-distutils \ + python3-mpi4py \ + python3-numpy \ + python3-pip \ + python3-tk \ + tcl-dev \ + tk-dev + RUN python3 -m pip install pytest-mpi pytest-check pytest ADD . /root/src/netgen