diff --git a/CMakeLists.txt b/CMakeLists.txt index a2cc9eab..c7125b6f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -254,6 +254,13 @@ if (USE_GUI) add_definitions(-DTCL -DOPENGL -DUSE_TOGL_2 -DUSE_TCL_STUBS -DUSE_TK_STUBS) include_directories(${TCL_INCLUDE_PATH}) include_directories(${TK_INCLUDE_PATH}) + if(NOT EXISTS ${TK_INCLUDE_PATH}/tkWin.h AND EXISTS ${TK_INCLUDE_PATH}/../win/tkWin.h) + include_directories(${TK_INCLUDE_PATH}/../win) + endif() + if(NOT EXISTS ${TK_INCLUDE_PATH}/x11/Xlib.h AND EXISTS ${TK_INCLUDE_PATH}/../xlib/X11/Xlib.h) + include_directories(${TK_INCLUDE_PATH}/../xlib) + endif() + set(LIBTOGL togl) if(WIN32) diff --git a/cmake/SuperBuild.cmake b/cmake/SuperBuild.cmake index c5a1b01b..77ffb231 100644 --- a/cmake/SuperBuild.cmake +++ b/cmake/SuperBuild.cmake @@ -50,7 +50,6 @@ set_vars(SUBPROJECT_CMAKE_ARGS CMAKE_C_COMPILER) set_vars(SUBPROJECT_CMAKE_ARGS CMAKE_CXX_COMPILER) set_vars(SUBPROJECT_CMAKE_ARGS CMAKE_BUILD_TYPE) -set(SUBPROJECT_CMAKE_ARGS "${SUBPROJECT_CMAKE_ARGS};-DCMAKE_OSX_ARCHITECTURES=${CMAKE_OSX_ARCHITECTURES}" CACHE INTERNAL "") set(SUBPROJECT_CMAKE_ARGS "${SUBPROJECT_CMAKE_ARGS};-DCMAKE_POSITION_INDEPENDENT_CODE=ON" CACHE INTERNAL "") if(USE_CCACHE) @@ -101,8 +100,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 +115,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/cmake/external_projects/tcltk.cmake b/cmake/external_projects/tcltk.cmake index 1d9e180a..637c1899 100644 --- a/cmake/external_projects/tcltk.cmake +++ b/cmake/external_projects/tcltk.cmake @@ -35,7 +35,7 @@ ExternalProject_Add(project_tk ) set(TCL_INCLUDE_PATH ${TCL_DIR}/generic) -set(TK_INCLUDE_PATH ${TK_DIR}/generic;${TK_DIR}/xlib;${TK_DIR}/win) +set(TK_INCLUDE_PATH ${TK_DIR}/generic) list(APPEND NETGEN_DEPENDENCIES project_tcl project_tk) if(APPLE OR WIN32) diff --git a/external_dependencies/pybind11 b/external_dependencies/pybind11 index c16da993..e7e2c79f 160000 --- a/external_dependencies/pybind11 +++ b/external_dependencies/pybind11 @@ -1 +1 @@ -Subproject commit c16da993094988141101ac4d96a9b2c92f9ac714 +Subproject commit e7e2c79f3f520f78ffc39fcb34f7919003102733 diff --git a/libsrc/core/python_ngcore_export.cpp b/libsrc/core/python_ngcore_export.cpp index 9cb04d11..b53c6f47 100644 --- a/libsrc/core/python_ngcore_export.cpp +++ b/libsrc/core/python_ngcore_export.cpp @@ -37,12 +37,14 @@ PYBIND11_MODULE(pyngcore, m) // NOLINT .def("__len__", &BitArray::Size) .def("__getitem__", [] (BitArray & self, int i) { + if (i < 0) i+=self.Size(); if (i < 0 || i >= self.Size()) throw py::index_error(); return self.Test(i); }, py::arg("pos"), "Returns bit from given position") .def("__setitem__", [] (BitArray & self, int i, bool b) { + if (i < 0) i+=self.Size(); if (i < 0 || i >= self.Size()) throw py::index_error(); if (b) self.SetBit(i); else self.Clear(i); diff --git a/libsrc/geom2d/csg2d.cpp b/libsrc/geom2d/csg2d.cpp index dfee7fba..357d6582 100644 --- a/libsrc/geom2d/csg2d.cpp +++ b/libsrc/geom2d/csg2d.cpp @@ -2185,6 +2185,7 @@ shared_ptr CSG2d :: GenerateSplineGeometry() if(!is_solid_degenerated) { geo->SetMaterial(dom, s.name); + geo->SetDomainMaxh(dom, s.maxh); if(s.layer != 1) geo->SetDomainLayer(dom, s.layer); } diff --git a/libsrc/geom2d/csg2d.hpp b/libsrc/geom2d/csg2d.hpp index fbe705a9..dcc06132 100644 --- a/libsrc/geom2d/csg2d.hpp +++ b/libsrc/geom2d/csg2d.hpp @@ -632,6 +632,7 @@ struct Solid2d int layer = 1; string name = MAT_DEFAULT; + double maxh = MAXH_DEFAULT; Solid2d() = default; Solid2d(string name_) : name(name_) {} @@ -697,6 +698,7 @@ struct Solid2d Solid2d & Maxh(double maxh) { + this->maxh = maxh; for(auto & p : polys) for(auto v : p.Vertices(ALL)) v->info.maxh = maxh; diff --git a/libsrc/geom2d/genmesh2d.cpp b/libsrc/geom2d/genmesh2d.cpp index 2ba08b9c..b2ed2dc3 100644 --- a/libsrc/geom2d/genmesh2d.cpp +++ b/libsrc/geom2d/genmesh2d.cpp @@ -429,7 +429,7 @@ namespace netgen PrintMessage (1, "Generate Mesh from spline geometry"); Box<2> bbox = geometry.GetBoundingBox (); - + bbox.Increase (1e-2*bbox.Diam()); t_h.Start(); if (bbox.Diam() < mp.maxh) mp.maxh = bbox.Diam(); 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..c2834d42 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -1,16 +1,112 @@ +#include + #include #include "meshing.hpp" #include 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(trafo(GetCenter()), e.GetCenter())) + return false; + + auto v0 = trafo(GetStartVertex().GetPoint()); + auto v1 = trafo(GetEndVertex().GetPoint()); + auto w0 = e.GetStartVertex().GetPoint(); + auto w1 = e.GetEndVertex().GetPoint(); + + // have two closed edges, use midpoints to compare + if(Dist(v0,v1) < tol && Dist(w0,w1) < tol) + { + v1 = trafo(GetPoint(0.5)); + w1 = other_ptr->GetPoint(0.5); + } + + return( (Dist(v0, w0) < tol && Dist(v1, w1) < tol) || + (Dist(v0, w1) < tol && Dist(v1, w0) < tol) ); + } + + bool GeometryFace :: IsMappedShape( const GeometryShape & other_, const Transformation<3> & trafo, double tol ) const + { + const auto other_ptr = dynamic_cast(&other_); + if(!other_ptr) + return false; + auto & f = *other_ptr; + + if(tol < Dist(GetCenter(), f.GetCenter())) + return false; + + // simple check: check if there is a bijective mapping of mapped edges + auto & other_edges = f.edges; + if(edges.Size() != other_edges.Size()) + return false; + + auto nedges = edges.Size(); + Array is_mapped(nedges); + is_mapped = false; + + for(auto e : edges) + { + int found_mapping = 0; + for(auto other_e : other_edges) + if(e->IsMappedShape(*other_e, trafo, tol)) + found_mapping++; + if(found_mapping != 1) + return false; + } + + return true; + } + void GeometryFace :: RestrictHTrig(Mesh& mesh, const PointGeomInfo& gi0, const PointGeomInfo& gi1, @@ -103,6 +199,123 @@ namespace netgen } }; + void NetgenGeometry :: Clear() + { + vertex_map.clear(); + edge_map.clear(); + face_map.clear(); + solid_map.clear(); + + vertices.SetSize0(); + edges.SetSize0(); + faces.SetSize0(); + solids.SetSize0(); + } + + void NetgenGeometry :: ProcessIdentifications() + { + for(auto i : Range(vertices)) + vertices[i]->nr = i; + for(auto i : Range(edges)) + edges[i]->nr = i; + for(auto i : Range(faces)) + faces[i]->nr = i; + for(auto i : Range(solids)) + solids[i]->nr = i; + + 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); + } + }; + + auto tol = 1e-8 * bounding_box.Diam(); + for(auto & f : faces) + for(auto & ident: f->identifications) + for(auto e : static_cast(ident.from)->edges) + for(auto e_other : static_cast(ident.to)->edges) + if(e->IsMappedShape(*e_other, ident.trafo, tol)) + e->identifications.Append( {e, e_other, ident.trafo, ident.type, ident.name} ); + + for(auto & e : edges) + for(auto & ident: e->identifications) + { + auto & from = static_cast(*ident.from); + auto & to = static_cast(*ident.to); + + GeometryVertex * pfrom[] = { &from.GetStartVertex(), &from.GetEndVertex() }; + GeometryVertex * pto[] = { &to.GetStartVertex(), &to.GetEndVertex() }; + + // swap points of other edge if necessary + Point<3> p_from0 = ident.trafo(from.GetStartVertex().GetPoint()); + Point<3> p_from1 = ident.trafo(from.GetEndVertex().GetPoint()); + Point<3> p_to0 = to.GetStartVertex().GetPoint(); + + if(Dist(p_from1, p_to0) < Dist(p_from0, p_to0)) + swap(pto[0], pto[1]); + + for(auto i : Range(2)) + pfrom[i]->identifications.Append( {pfrom[i], pto[i], ident.trafo, ident.type, ident.name} ); + } + + mirror_identifications(vertices); + mirror_identifications(edges); + mirror_identifications(faces); + + + 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 < s->primary->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; + } + if(other->primary->nr < s->primary->nr) + { + auto trafo = ident.trafo; + if(need_inverse) + trafo = trafo.CalcInverse(); + s->primary = other->primary; + s->primary_to_me.Combine(trafo, other->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 +454,285 @@ 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 = 10000; + 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 boundary = face.GetBoundary(facebndnr); - for(auto enr : Range(boundary)) - { - multithread.percent = 100. * ((double(enr)/boundary.Size() + facebndnr)/face.GetNBoundaries() + facenr)/faces.Size(); - const auto& oriented_edge = *boundary[enr]; - auto edgenr = GetEdgeIndex(oriented_edge); - const auto& edge = edges[edgenr]; - PointIndex startp, endp; - // throws if points are not found - startp = vert2meshpt.at(edge->GetStartVertex().GetHash()); - endp = vert2meshpt.at(edge->GetEndVertex().GetHash()); + auto nedges = edges.Size(); + Array> all_pnums(nedges); + Array> all_params(nedges); - // 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; + 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()); - 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); + // ignore collapsed edges + if(startp == endp && edge->GetLength() < 1e-10 * bounding_box.Diam()) + continue; - 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; - } + // ----------- Add Points to mesh and create segments ----- + auto & pnums = all_pnums[edgenr]; + auto & params = all_params[edgenr]; + Array> edge_points; + Array edge_params; - } while(i < nsubedges); + if(edge->primary == edge) + { + // check if start and end vertex are identified (if so, we only insert one segement and do z-refinement later) + bool is_identified_edge = false; + auto v0 = vertices[edge->GetStartVertex().nr].get(); + auto v1 = vertices[edge->GetEndVertex().nr].get(); + for(auto & ident : v0->identifications) + { + auto other = ident.from == v0 ? ident.to : ident.from; + if(other->nr == v1->nr && ident.type == Identifications::CLOSESURFACES) + { + is_identified_edge = true; + break; + } + } + if(is_identified_edge) + { + params.SetSize(2); params[0] = 0.; - params[nsubedges] = 1.; + params[1] = 1.; + } + else + { + DivideEdge(edge, mparam, mesh, edge_points, 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; - 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; + 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; + } - double eps = bounding_box.Diam() * 1e-8; + // 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]); + } - 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]); - } + params.SetSize(edge_params.Size()+2); + params[0] = 0.; + params.Last() = 1.; - 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; + for(auto i : Range(edge_params)) + params[i+1] = edge_params[i]; + } - // 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]); + pnums.SetSize(edge_points.Size() + 2); + pnums[0] = startp; + pnums.Last() = endp; - 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); - } + + for(auto i : Range(edge_points)) + { + auto pi = mesh.AddPoint(edge_points[i]); + tree.Insert(mesh[pi], pi); + pnums[i+1] = pi; + } + + 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; + seg.domin = edge->domin+1; + seg.domout = edge->domout+1; + 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 & 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) + { + 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, @@ -408,66 +741,239 @@ namespace netgen static Timer t1("Surface Meshing"); RegionTimer regt(t1); const char* savetask = multithread.task; multithread.task = "Mesh Surface"; + mesh.ClearFaceDescriptors(); + 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]; + FaceDescriptor fd(k+1, face.domin+1, face.domout+1, k+1); + mesh.AddFaceDescriptor(fd); + mesh.SetBCName(k, face.properties.GetName()); + if(face.primary == &face) + { + // check if this face connects two identified closesurfaces + bool is_connecting_closesurfaces = false; + auto & idents = mesh.GetIdentifications(); + std::set relevant_edges; + auto segments = face.GetBoundary(mesh); + for(const auto &s : segments) + relevant_edges.insert(s.edgenr-1); - 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); - } - } + Array is_point_in_tree(mesh.Points().Size()); + is_point_in_tree = false; + PointTree tree( bounding_box ); + for(const auto &s : segments) + for(auto pi : s.PNums()) + if(!is_point_in_tree[pi]) + { + tree.Insert(mesh[pi], pi); + is_point_in_tree[pi] = true; + } - // TODO Set max area 2* area of face + Array mapped_edges(edges.Size()); + constexpr int UNINITIALIZED = -2; + constexpr int NOT_MAPPED = -1; + mapped_edges = UNINITIALIZED; - auto noldsurfels = mesh.GetNSE(); + Transformation<3> trafo; + for(const auto &s : segments) + { + auto edgenr = s.edgenr-1; + auto & edge = *edges[edgenr]; + ShapeIdentification *edge_mapping; - static Timer t("GenerateMesh"); RegionTimer reg(t); - MESHING2_RESULT res = meshing.GenerateMesh(mesh, mparam, mparam.maxh, k+1); + // have edgenr first time, search for closesurface identification + if(mapped_edges[edgenr] == UNINITIALIZED) + { + mapped_edges[edgenr] = NOT_MAPPED; + for(auto & edge_ident : edge.identifications) + { + if(edge_ident.type == Identifications::CLOSESURFACES && + edge_ident.from->nr == edgenr && + relevant_edges.count(edge_ident.to->nr) > 0 + ) + { + trafo = edge_ident.trafo; + mapped_edges[edgenr] = edge_ident.to->nr; + is_connecting_closesurfaces = true; + break; + } + } + } - for(auto i : Range(noldsurfels, mesh.GetNSE())) - { - mesh.SurfaceElements()[i].SetIndex(k+1); - } - } + // this edge has a closesurface mapping to another -> make connecting quad + if(mapped_edges[edgenr] != NOT_MAPPED) + { + Element2d sel(4); + sel[0] = s[0]; + sel[1] = s[1]; + sel[2] = tree.Find(trafo(mesh[s[1]])); + sel[3] = tree.Find(trafo(mesh[s[0]])); + for(auto i : Range(4)) + sel.GeomInfo()[i] = face.Project(mesh[sel[i]]); + + sel.SetIndex(face.nr+1); + mesh.AddSurfaceElement(sel); + } + } + + if(!is_connecting_closesurfaces) + if(MeshFace(mesh, mparam, k, glob2loc)) + n_failed_faces++; + } + } + + 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; + } + + if (mparam.perfstepsend >= MESHCONST_OPTSURFACE) + { + mesh.CalcSurfacesOfNode(); + OptimizeSurface(mesh, mparam); + } + + 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); + } + + 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); + } + } + } + + 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; + } + + xbool do_invert = maybe; + + // now insert mapped surface elements + for(auto sei : mesh.SurfaceElements().Range()) + { + auto sel = mesh[sei]; + if(sel.GetIndex() != src.nr+1) + continue; + + if(do_invert.IsMaybe()) + { + auto n_src = src.GetNormal(mesh[sel[0]]); + auto n_dist = dst.GetNormal(mesh[sel[0]]); + Mat<3> normal_matrix; + CalcInverse(Trans(trafo.GetMatrix()), normal_matrix); + do_invert = n_src * (normal_matrix * n_dist) < 0.0; + } + 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]; + } + if(do_invert.IsTrue()) + 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 +1010,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 +1080,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..49bf8710 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -12,26 +12,82 @@ 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: - virtual ~GeometryEdge() {} - virtual const GeometryVertex& GetStartVertex() const = 0; - virtual const GeometryVertex& GetEndVertex() const = 0; + 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 + { + protected: + GeometryVertex *start, *end; + public: + int domin=-1, domout=-1; + + GeometryEdge( GeometryVertex &start_, GeometryVertex &end_ ) + : start(&start_), end(&end_) + {} + + virtual const GeometryVertex& GetStartVertex() const { return *start; } + virtual const GeometryVertex& GetEndVertex() const { return *end; } + virtual GeometryVertex& GetStartVertex() { return *start; } + virtual GeometryVertex& GetEndVertex() { return *end; } 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 +102,19 @@ 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() {} + Array edges; + 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. @@ -93,6 +153,8 @@ namespace netgen newgi = Project(newp); } + virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const override; + protected: void RestrictHTrig(Mesh& mesh, const PointGeomInfo& gi0, @@ -102,6 +164,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 +174,31 @@ namespace netgen Array> vertices; Array> edges; Array> faces; + Array> solids; Array, double>> restricted_h; Box<3> bounding_box; + int dimension = 3; + + std::map vertex_map; + std::map edge_map; + std::map face_map; + std::map solid_map; 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]; } + + void Clear(); + virtual int GenerateMesh (shared_ptr & mesh, MeshingParameters & mparam); void RestrictH(const Point<3>& pnt, double maxh) @@ -134,13 +215,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/meshclass.cpp b/libsrc/meshing/meshclass.cpp index efe51f4d..51a9c9f0 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -665,6 +665,13 @@ namespace netgen outfile << "geomtype\n" << int(geomtype) << "\n"; + outfile << "\n"; + outfile << "# surfnr\tdomin\tdomout\ttlosurf\tbcprop\n"; + outfile << "facedescriptors\n"; + outfile << GetNFD() << "\n"; + for(auto & fd : FaceDescriptors()) + outfile << fd.SurfNr() << ' ' << fd.DomainIn() << ' ' << fd.DomainOut() << ' ' << fd.TLOSurface() << ' ' << fd.BCProperty() << '\n'; + outfile << "\n"; outfile << "# surfnr bcnr domin domout np p1 p2 p3" @@ -1192,6 +1199,19 @@ namespace netgen geomtype = GEOM_TYPE(hi); } + if (strcmp (str, "facedescriptors") == 0) + { + int nfd; + infile >> nfd; + for(auto i : Range(nfd)) + { + int surfnr, domin, domout, tlosurf, bcprop; + infile >> surfnr >> domin >> domout >> tlosurf >> bcprop; + auto faceind = AddFaceDescriptor (FaceDescriptor(surfnr, domin, domout, tlosurf)); + GetFaceDescriptor(faceind).SetBCProperty(bcprop); + } + } + if (strcmp (str, "surfaceelements") == 0 || strcmp (str, "surfaceelementsgi")==0 || strcmp (str, "surfaceelementsuv") == 0) { @@ -6887,14 +6907,11 @@ namespace netgen int oldsize = bcnames.Size(); bcnames.SetSize (bcnr+1); // keeps contents for (int i = oldsize; i <= bcnr; i++) - bcnames[i] = nullptr; + bcnames[i] = new string("default"); } if ( bcnames[bcnr] ) delete bcnames[bcnr]; - if ( abcname != "default" ) - bcnames[bcnr] = new string ( abcname ); - else - bcnames[bcnr] = nullptr; + bcnames[bcnr] = new string ( abcname ); for (auto & fd : facedecoding) if (fd.BCProperty() <= bcnames.Size()) diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 597a61fa..68dcd347 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -156,8 +156,6 @@ namespace netgen if(!pi0.IsValid() || !pi1.IsValid()) continue; - if(pi1GetIdentifications(); + auto nmax = identifications.GetMaxNr(); + + bool have_closesurfaces = false; + for(auto i : Range(1,nmax+1)) + if(identifications.GetType(i) == Identifications::CLOSESURFACES) + have_closesurfaces = true; + if(!have_closesurfaces) + return; + + NgArray map; + for(auto identnr : Range(1,nmax+1)) + { + if(identifications.GetType(identnr) != Identifications::CLOSESURFACES) + continue; + + identifications.GetMap(identnr, map); + + for(auto & sel : mesh->SurfaceElements()) + { + bool is_mapped = true; + for(auto pi : sel.PNums()) + if(!PointIndex(map[pi]).IsValid()) + is_mapped = false; + + if(!is_mapped) + continue; + + // in case we have symmetric mapping (used in csg), only map in one direction + if(map[map[sel[0]]] == sel[0] && map[sel[0]] < sel[0]) + continue; + + // insert prism + auto np = sel.GetNP(); + Element el(2*np); + for(auto i : Range(np)) + { + el[i] = sel[i]; + el[i+np] = map[sel[i]]; + } + + el.SetIndex(md.domain); + mesh->AddVolumeElement(el); + } + } + } + void CloseOpenQuads( MeshingData & md) { auto & mesh = *md.mesh; @@ -488,6 +539,9 @@ namespace netgen if (md[i].mesh->CheckOverlappingBoundary()) throw NgException ("Stop meshing since boundary mesh is overlapping"); + // TODO: FillCloseSurface is not working with CSG closesurfaces + if(md[i].mesh->GetGeometry()->GetGeomType() == Mesh::GEOM_OCC) + FillCloseSurface( md[i] ); CloseOpenQuads( md[i] ); MeshDomain(md[i]); }); 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..8881123b 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -117,6 +117,15 @@ 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::implicitly_convertible(); + py::class_ (m, "MPI_Comm") #ifdef NG_MPI4PY .def(py::init([] (mpi4py_comm comm) @@ -944,19 +953,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_, GeometryVertex & start_, GeometryVertex & end_) + : GeometryEdge(start_, end_), + 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"); + + if(start != end) + { + // 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); + } + } + + 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() - s0)/(s1-s0); + 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..f34cee56 --- /dev/null +++ b/libsrc/occ/occ_edge.hpp @@ -0,0 +1,40 @@ +#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 + { + public: + T_Shape tedge; + TopoDS_Edge edge; + Handle(Geom_Curve) curve; + double s0, s1; + GProp_GProps props; + + public: + OCCEdge(TopoDS_Shape edge_, GeometryVertex & start_, GeometryVertex & end_); + + auto Shape() const { return edge; } + T_Shape TShape() const { return tedge; } + + 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..7b985f06 --- /dev/null +++ b/libsrc/occ/occ_face.cpp @@ -0,0 +1,252 @@ +#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]; + double s0, s1; + auto cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1); + + double s[2] = { seg.epgeominfo[0].dist, seg.epgeominfo[1].dist }; + + // dist is in [0,1], map parametrization to [s0, s1] + s[0] = s0 + s[0]*(s1-s0); + s[1] = s0 + s[1]*(s1-s0); + + // 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..210e358e --- /dev/null +++ b/libsrc/occ/occ_utils.cpp @@ -0,0 +1,40 @@ +#include +#include +#include + +#include "occ_utils.hpp" + +namespace netgen +{ + Point<3> occ2ng (Handle(TopoDS_TShape) shape) + { + return occ2ng( Handle(BRep_TVertex)::DownCast(shape)->Pnt() ); + } + + Transformation<3> occ2ng (const gp_Trsf & occ_trafo) + { + Transformation<3> trafo; + auto v = occ_trafo.TranslationPart(); + auto m = occ_trafo.VectorialPart(); + auto & tv = trafo.GetVector(); + auto & tm = trafo.GetMatrix(); + for(auto i : Range(3)) + { + tv[i] = v.Coord(i+1); + for(auto k : Range(3)) + tm(i,k) = m(i+1,k+1); + } + return trafo; + } + + 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..e51a2b51 --- /dev/null +++ b/libsrc/occ/occ_utils.hpp @@ -0,0 +1,283 @@ +#ifndef FILE_OCC_UTILS_INCLUDED +#define FILE_OCC_UTILS_INCLUDED + +#include +#include +#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)); + } + + DLL_HEADER Transformation<3> occ2ng (const gp_Trsf & t); + + 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 Wires() const + { + return SubShapes(TopAbs_WIRE); + } + 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 GetWires(const TopoDS_Shape & shape) + { + ListOfShapes sub; + for (TopExp_Explorer e(shape, TopAbs_WIRE); 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..643d1d9a 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,201 @@ 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, " / ", geom.GetNFaces(), " (plane space projection)"); + else + PrintMessage (2, "Face ", k, " / ", geom.GetNFaces(), " (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(); + static Timer t("MeshSurface: Find edges and points - Parameter"); RegionTimer r(t); + + Array gis(2*segments.Size()); + gis.SetSize (0); + + Box<2> uv_box(Box<2>::EMPTY_BOX); + for(auto & seg : segments) + for(auto i : Range(2)) + uv_box.Add( {seg.epgeominfo[i].u, seg.epgeominfo[i].v } ); + + BoxTree<2> uv_tree(uv_box); + Array found_points; + + for(auto & seg : segments) + { + PointGeomInfo gi[2]; + gi[0].trignum = gi[1].trignum = k; + gi[0].u = seg.epgeominfo[0].u; + gi[0].v = seg.epgeominfo[0].v; + gi[1].u = seg.epgeominfo[1].u; + gi[1].v = seg.epgeominfo[1].v; + + int locpnum[2] = {0, 0}; + + for (int j = 0; j < 2; j++) + { + Point<2> uv = {gi[j].u, gi[j].v}; + uv_tree.GetIntersecting(uv, uv, found_points); + + if(found_points.Size()) + locpnum[j] = found_points[0]; + else + { + PointIndex pi = seg[j]; + meshing.AddPoint (mesh.Point(pi), pi); + + gis.Append (gi[j]); + locpnum[j] = gis.Size(); + uv_tree.Insert(uv, locpnum[j]); + } + } + + meshing.AddBoundaryElement (locpnum[0], locpnum[1], gi[0], gi[1]); + } } + + + // 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; - for (int i = 0; i < mesh.GetNFD(); i++) - mesh.SetBCName (i, mesh.GetFaceDescriptor(i+1).GetBCName()); - multithread.task = savetask; - } + try { + static Timer t("GenerateMesh"); RegionTimer reg(t); + res = meshing.GenerateMesh (mesh, mparam_without_overlap, maxh, k); + } - void OCCOptimizeSurface(OCCGeometry & geom, Mesh & mesh, - const MeshingParameters & mparam) - { - const char * savetask = multithread.task; - multithread.task = "Optimizing surface"; - - static Timer timer_opt2d("Optimization 2D"); - timer_opt2d.Start(); - - for (int k = 1; k <= mesh.GetNFD(); k++) + catch (SingularMatrixException) { - // if (k != 42) continue; - // if (k != 36) continue; - - // (*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++) - { - // (*testout) << "optstep " << i << endl; - if (multithread.terminate) return; - - { - MeshOptimize2d meshopt(mesh); - meshopt.SetFaceIndex (k); - meshopt.SetImproveEdges (0); - meshopt.SetMetricWeight (mparam.elsizeweight); - meshopt.SetWriteStatus (0); - meshopt.EdgeSwapping (i > mparam.optsteps2d/2); - } - - if (multithread.terminate) return; - { - MeshOptimize2d meshopt(mesh); - meshopt.SetFaceIndex (k); - meshopt.SetImproveEdges (0); - meshopt.SetMetricWeight (mparam.elsizeweight); - meshopt.SetWriteStatus (0); - meshopt.ImproveMesh (mparam); - } - - { - MeshOptimize2d meshopt(mesh); - meshopt.SetFaceIndex (k); - meshopt.SetImproveEdges (0); - meshopt.SetMetricWeight (mparam.elsizeweight); - meshopt.SetWriteStatus (0); - meshopt.CombineImprove (); - } - - if (multithread.terminate) return; - { - MeshOptimize2d meshopt(mesh); - meshopt.SetFaceIndex (k); - meshopt.SetImproveEdges (0); - meshopt.SetMetricWeight (mparam.elsizeweight); - meshopt.SetWriteStatus (0); - meshopt.ImproveMesh (mparam); - } - } - + // (*myerr) << "Singular Matrix" << endl; + res = MESHING2_GIVEUP; } + catch (UVBoundsException) + { + // (*myerr) << "UV bounds exceeded" << endl; + res = MESHING2_GIVEUP; + } - mesh.CalcSurfacesOfNode(); - mesh.Compress(); - timer_opt2d.Stop(); + 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); - multithread.task = savetask; + 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..f6d38a2b 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, PARAMETERSPACE, true); + if(failed) + failed = OCCMeshFace(*this, mesh, glob2loc, mparam, nr, PLANESPACE, 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,125 @@ namespace netgen fsingular = esingular = vsingular = false; + NetgenGeometry::Clear(); + edge_map.clear(); + vertex_map.clear(); + face_map.clear(); + solid_map.clear(); - sprops.SetSize(somap.Extent()); - for (TopExp_Explorer e(shape, TopAbs_SOLID); e.More(); e.Next()) + // Add shapes + for(auto i1 : Range(1, vmap.Extent()+1)) { - auto s = e.Current(); - sprops[somap.FindIndex(s)-1] = &global_shape_properties[s.TShape()]; + auto v = vmap(i1); + auto tshape = v.TShape(); + if(vertex_map.count(tshape)!=0) + continue; + auto occ_vertex = make_unique(TopoDS::Vertex(v)); + occ_vertex->nr = vertices.Size(); + vertex_map[tshape] = occ_vertex->nr; + + 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 i1 : Range(1, emap.Extent()+1)) { - auto s = e.Current(); - fprops[fmap.FindIndex(s)-1] = &global_shape_properties[s.TShape()]; + auto e = emap(i1); + auto tshape = e.TShape(); + auto edge = TopoDS::Edge(e); + if(edge_map.count(tshape)!=0) + continue; + edge_map[tshape] = edges.Size(); + auto verts = GetVertices(e); + auto occ_edge = make_unique(edge, *vertices[vertex_map[verts[0].TShape()]], *vertices[vertex_map[verts[1].TShape()]] ); + 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 i1 : Range(1, fmap.Extent()+1)) { - auto s = e.Current(); - eprops[emap.FindIndex(s)-1] = &global_shape_properties[s.TShape()]; + auto f = fmap(i1); + auto tshape = f.TShape(); + if(face_map.count(tshape)==0) + { + + auto k = faces.Size(); + face_map[tshape] = k; + auto occ_face = make_unique(f); + + for(auto e : GetEdges(f)) + occ_face->edges.Append( edges[edge_map[e.TShape()]].get() ); + + if(global_shape_properties.count(tshape)>0) + occ_face->properties = global_shape_properties[tshape]; + faces.Append(std::move(occ_face)); + + if(dimension==2) + for(auto e : GetEdges(f)) + { + auto & edge = *edges[edge_map[e.TShape()]]; + if(e.Orientation() == TopAbs_REVERSED) + edge.domout = k; + else + edge.domin = k; + } + } } - */ - for (auto [nr,s] : Enumerate(emap)) - eprops[nr-1] = &global_shape_properties[s.TShape()]; + + + for(auto i1 : Range(1, somap.Extent()+1)) + { + auto s = somap(i1); + 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]) + { + if(shape_map.count(ident.from)==0 || shape_map.count(ident.to)==0) + continue; + + 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 ); + + bounding_box = ::netgen::GetBoundingBox( shape ); + ProcessIdentifications(); } @@ -1156,286 +1328,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 +1596,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 +1630,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 +1948,114 @@ namespace netgen return false; } + Point<3> GetCenter(const TopoDS_Shape & shape) + { + GProp_GProps props; + BRepGProp::LinearProperties(shape, props); + return occ2ng( props.CentreOfMass() ); + } + + 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); + auto verts_you = GetVertices(you); + + if(verts_me.size() != verts_you.size()) + return false; + + for (auto i : Range(verts_me.size())) + { + auto s = verts_me[i].TShape(); + if(vmap.count(s)>0) + continue; + auto p = trafo(occ2ng(s)); + tree.Insert( p, i ); + vmap[s] = nullptr; + } + + for (auto vert : verts_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) + return false; + } + return true; + } + + void Identify(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo) + { + gp_Trsf trafo; + if(opt_trafo) + { + trafo = *opt_trafo; + } + else + { + auto v = GetCenter(you) - GetCenter(me); + trafo.SetTranslation(gp_Vec(v[0], v[1], v[2])); + } + + ListOfShapes list_me, list_you; + list_me.push_back(me); + list_you.push_back(you); + Identify(list_me, list_you, name, type, trafo); + } + + void Identify(const ListOfShapes & me, const ListOfShapes & you, string name, Identifications::ID_TYPE type, gp_Trsf occ_trafo) + { + Transformation<3> trafo = occ2ng(occ_trafo); + + ListOfShapes id_me; + ListOfShapes id_you; + + if(auto faces_me = me.Faces(); faces_me.size()>0) + { + id_me = faces_me; + id_you = you.Faces(); + } + else if(auto edges_me = me.Edges(); edges_me.size()>0) + { + id_me = edges_me; + id_you = you.Edges(); + } + else + { + id_me = me.Vertices(); + id_you = you.Vertices(); + } + + for(auto shape_me : id_me) + for(auto shape_you : id_you) + { + if(!IsMappedShape(trafo, shape_me, shape_you)) + continue; + + OCCGeometry::identifications[shape_me.TShape()].push_back + (OCCIdentification { shape_me.TShape(), shape_you.TShape(), trafo, name, type }); + } + } + void OCCParameters :: Print(ostream & ost) const { ost << "OCC Parameters:" << endl @@ -2063,6 +2066,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 +2223,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 +2253,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..372bcc06 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -9,100 +9,30 @@ #ifdef OCCGEOMETRY +#include + #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; @@ -116,28 +46,7 @@ namespace netgen #define OCCGEOMETRYVISUALIZATIONFULLCHANGE 1 // Compute transformation matrices and redraw #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)); - } - - + bool IsMappedShape(const Transformation<3> & trafo, const TopoDS_Shape & me, const TopoDS_Shape & you); class EntityVisualizationCode { @@ -219,108 +128,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 +171,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 +191,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; } @@ -548,8 +346,11 @@ namespace netgen // void WriteOCC_STL(char * filename); 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; }; + + void Identify(const ListOfShapes & me, const ListOfShapes & you, string name, Identifications::ID_TYPE type, gp_Trsf occ_trafo); + void Identify(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo); void PrintContents (OCCGeometry * geom); @@ -564,12 +365,107 @@ 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 bool OCCMeshFace (const OCCGeometry & geom, Mesh & mesh, FlatArray glob2loc, + const MeshingParameters & mparam, int nr, int projecttype, bool delete_on_failure); - 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); + template + void PropagateIdentifications (TBuilder & builder, TopoDS_Shape shape, std::optional> trafo = nullopt) + { + std::map> mod_map; + std::map tshape_handled; + Transformation<3> trafo_inv; + if(trafo) + trafo_inv = trafo->CalcInverse(); + + for (auto typ : { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX }) + for (TopExp_Explorer e(shape, typ); e.More(); e.Next()) + { + auto tshape = e.Current().TShape(); + mod_map[tshape].insert(tshape); + tshape_handled[tshape] = false; + } + + for (auto typ : { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX }) + for (TopExp_Explorer e(shape, typ); e.More(); e.Next()) + { + auto tshape = e.Current().TShape(); + + for (auto mods : builder.Modified(e.Current())) + mod_map[tshape].insert(mods.TShape()); + } + + for (auto typ : { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX }) + for (TopExp_Explorer e(shape, typ); e.More(); e.Next()) + { + auto tshape = e.Current().TShape(); + + if(tshape_handled[tshape]) + continue; + tshape_handled[tshape] = true; + + if(OCCGeometry::identifications.count(tshape)==0) + continue; + + auto tshape_mapped = mod_map[tshape]; + + for(auto ident : OCCGeometry::identifications[tshape]) + { + // nothing happened + if(mod_map[ident.to].size()==1 && mod_map[ident.from].size() ==1) + continue; + + auto from = ident.from; + auto to = ident.to; + + for(auto from_mapped : mod_map[from]) + for(auto to_mapped : mod_map[to]) + { + if(from==from_mapped && to==to_mapped) + continue; + + TopoDS_Shape s_from; s_from.TShape(from_mapped); + TopoDS_Shape s_to; s_to.TShape(to_mapped); + Transformation<3> trafo_mapped = ident.trafo; + if(trafo) + { + Transformation<3> trafo_temp; + trafo_temp.Combine(ident.trafo, trafo_inv); + trafo_mapped.Combine(*trafo, trafo_temp); + } + + if(!IsMappedShape(trafo_mapped, s_from, s_to)) + continue; + + OCCIdentification id_new = ident; + id_new.to = to_mapped; + id_new.from = from_mapped; + id_new.trafo = trafo_mapped; + auto id_owner = from == tshape ? from_mapped : to_mapped; + OCCGeometry::identifications[id_owner].push_back(id_new); + } + } + } + } + + template + void PropagateProperties (TBuilder & builder, TopoDS_Shape shape, std::optional> trafo = nullopt) + { + bool have_identifications = false; + + for (auto typ : { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE }) + for (TopExp_Explorer e(shape, typ); e.More(); e.Next()) + { + auto tshape = e.Current().TShape(); + auto & prop = OCCGeometry::global_shape_properties[tshape]; + for (auto mods : builder.Modified(e.Current())) + OCCGeometry::global_shape_properties[mods.TShape()].Merge(prop); + have_identifications |= OCCGeometry::identifications.count(tshape) > 0; + } + if(have_identifications) + PropagateIdentifications(builder, shape, trafo); + } 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..bfb39cac 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; @@ -199,6 +166,10 @@ DLL_HEADER void ExportNgOCC(py::module &m) { self.SetFaceMaxH(fnr, meshsize); }, "Set maximum meshsize for face fnr. Face numbers are 0 based.") + .def("Draw", [](shared_ptr geo) + { + ng_geometry = geo; + }) .def("_visualizationData", [] (shared_ptr occ_geo) { std::vector vertices; @@ -294,12 +265,10 @@ DLL_HEADER void ExportNgOCC(py::module &m) geo->SetOCCParameters(occparam); auto mesh = make_shared(); mesh->SetGeometry(geo); + SetGlobalMesh(mesh); 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; }, py::arg("mp") = nullptr, 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..cdf338fe 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) { @@ -87,6 +49,16 @@ DLL_HEADER void ExportNgOCCBasic(py::module &m) .def("__sub__", [](gp_Pnt p1, gp_Pnt p2) { return gp_Vec(p2, p1); }) .def("__add__", [](gp_Pnt p, gp_Vec v) { return p.Translated(v); }) // gp_Pnt(p.X()+v.X(), p.Y()+v.Y(), p.Z()+v.Z()); }) .def("__sub__", [](gp_Pnt p, gp_Vec v) { return p.Translated(-v); }) // gp_Pnt(p.X()-v.X(), p.Y()-v.Y(), p.Z()-v.Z()); }) + .def("__getitem__", [](const gp_Pnt& p, int index) + { + if(index == 0) + return p.X(); + if(index == 1) + return p.Y(); + if(index == 2) + return p.Z(); + throw std::out_of_range("Point index must be in range [0,3)!"); + }) ; py::class_(m, "gp_Vec", "3d OCC vector") @@ -297,6 +269,7 @@ DLL_HEADER void ExportNgOCCBasic(py::module &m) py::class_(m, "gp_Trsf") .def(py::init<>()) .def("SetMirror", [] (gp_Trsf & trafo, const gp_Ax1 & ax) { trafo.SetMirror(ax); return trafo; }) + .def("Inverted", &gp_Trsf::Inverted) .def_static("Translation", [] (const gp_Vec & v) { gp_Trsf trafo; trafo.SetTranslation(v); return trafo; }) .def_static("Scale", [] (const gp_Pnt & p, double s) { gp_Trsf trafo; trafo.SetScale(p,s); return trafo; }) .def_static("Mirror", [] (const gp_Ax1 & ax) { gp_Trsf trafo; trafo.SetMirror(ax); return trafo; }) @@ -308,7 +281,9 @@ DLL_HEADER void ExportNgOCCBasic(py::module &m) { gp_Trsf trafo; trafo.SetTransformation(from, to); return trafo; }) .def(py::self * py::self) .def("__call__", [] (gp_Trsf & trafo, const TopoDS_Shape & shape) { - return BRepBuilderAPI_Transform(shape, trafo).Shape(); + BRepBuilderAPI_Transform builder(shape, trafo, true); + PropagateProperties(builder, shape, occ2ng(trafo)); + return builder.Shape(); }) .def("__str__", [](gp_Trsf & trafo) { diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index 126aa23e..6a819a0e 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -3,185 +3,76 @@ #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 +#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(); - // cout << "dist = " << dist << endl; - if (dist < mindist) - nearestshape = shape; - } - 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; @@ -312,23 +203,6 @@ py::object CastShape(const TopoDS_Shape & s) }; -template -void PropagateProperties (TBuilder & builder, TopoDS_Shape shape) -{ - // #ifdef OCC_HAVE_HISTORY - // Handle(BRepTools_History) history = builder.History(); - for (auto typ : { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE }) - for (TopExp_Explorer e(shape, typ); e.More(); e.Next()) - { - auto prop = OCCGeometry::global_shape_properties[e.Current().TShape()]; - // for (auto mods : history->Modified(e.Current())) - for (auto mods : builder.Modified(e.Current())) - OCCGeometry::global_shape_properties[mods.TShape()].Merge(prop); - } - // #endif -} - - class WorkPlane : public enable_shared_from_this { gp_Ax3 axes; @@ -819,36 +693,21 @@ 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("wires", GetWires, + "returns all sub-shapes of type 'WIRE'") + .def_property_readonly("vertices", GetVertices, + "returns all sub-shapes of type 'VERTEX'") + .def_property_readonly("bounding_box", [] ( const TopoDS_Shape &shape ) + { + auto box = GetBoundingBox(shape); + return py::make_tuple( ng2occ(box.PMin()), ng2occ(box.PMax()) ); + }, "returns bounding box (pmin, pmax)") .def("Properties", [] (const TopoDS_Shape & shape) { @@ -896,8 +755,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) // version 1: Transoformation gp_Trsf trafo; trafo.SetTranslation(v); - BRepBuilderAPI_Transform builder(shape, trafo); - PropagateProperties(builder, shape); + BRepBuilderAPI_Transform builder(shape, trafo, true); + PropagateProperties(builder, shape, occ2ng(trafo)); return builder.Shape(); // version 2: change location // ... @@ -908,8 +767,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) { gp_Trsf trafo; trafo.SetRotation(ax, ang*M_PI/180); - BRepBuilderAPI_Transform builder(shape, trafo); - PropagateProperties(builder, shape); + BRepBuilderAPI_Transform builder(shape, trafo, true); + PropagateProperties(builder, shape, occ2ng(trafo)); return builder.Shape(); }, py::arg("axis"), py::arg("ang"), "copy shape, and rotet copy by 'ang' degrees around 'axis'") @@ -918,8 +777,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) { gp_Trsf trafo; trafo.SetMirror(ax.Ax2()); - BRepBuilderAPI_Transform builder(shape, trafo); - PropagateProperties(builder, shape); + BRepBuilderAPI_Transform builder(shape, trafo, true); + PropagateProperties(builder, shape, occ2ng(trafo)); return builder.Shape(); }, py::arg("axes"), "copy shape, and mirror over plane defined by 'axes'") @@ -928,8 +787,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) { gp_Trsf trafo; trafo.SetMirror(ax); - BRepBuilderAPI_Transform builder(shape, trafo); - PropagateProperties(builder, shape); + BRepBuilderAPI_Transform builder(shape, trafo, true); + PropagateProperties(builder, shape, occ2ng(trafo)); return builder.Shape(); }, py::arg("axes"), "copy shape, and mirror around axis 'axis'") @@ -938,8 +797,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) { gp_Trsf trafo; trafo.SetScale(p, s); - BRepBuilderAPI_Transform builder(shape, trafo); - PropagateProperties(builder, shape); + BRepBuilderAPI_Transform builder(shape, trafo, true); + PropagateProperties(builder, shape, occ2ng(trafo)); return builder.Shape(); }, py::arg("p"), py::arg("s"), "copy shape, and scale copy by factor 's'") @@ -967,7 +826,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) return *name; else return string(); - }, [](const TopoDS_Shape & self, string name) { + }, [](const TopoDS_Shape & self, optional name) { OCCGeometry::global_shape_properties[self.TShape()].name = name; }, "'name' of shape") @@ -1138,16 +997,24 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) .def("Reversed", [](const TopoDS_Shape & shape) { return CastShape(shape.Reversed()); }) - .def("Extrude", [](const TopoDS_Shape & shape, double h) { + .def("Extrude", [](const TopoDS_Shape & shape, double h, + optional dir) { for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next()) { Handle(Geom_Surface) surf = BRep_Tool::Surface (TopoDS::Face(e.Current())); - gp_Vec du, dv; - gp_Pnt p; - surf->D1 (0,0,p,du,dv); - BRepPrimAPI_MakePrism builder(shape, h*du^dv); + gp_Vec edir; + if(dir.has_value()) + edir = *dir; + else + { + gp_Vec du, dv; + gp_Pnt p; + surf->D1 (0,0,p,du,dv); + edir = du^dv; + } + BRepPrimAPI_MakePrism builder(shape, h*edir, true); - for (auto typ : { TopAbs_EDGE, TopAbs_VERTEX }) + for (auto typ : { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX }) for (TopExp_Explorer e(shape, typ); e.More(); e.Next()) { auto prop = OCCGeometry::global_shape_properties[e.Current().TShape()]; @@ -1158,7 +1025,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) return builder.Shape(); } throw Exception("no face found for extrusion"); - }, py::arg("h"), "extrude shape to thickness 'h', shape must contain a plane surface") + }, py::arg("h"), py::arg("dir")=nullopt, "extrude shape to thickness 'h', shape must contain a plane surface, optionally give an extrusion direction") .def("Extrude", [] (const TopoDS_Shape & face, gp_Vec vec) { return BRepPrimAPI_MakePrism (face, vec).Shape(); @@ -1168,7 +1035,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) // for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next()) { // return BRepPrimAPI_MakeRevol (shape, A, D*M_PI/180).Shape(); - BRepPrimAPI_MakeRevol builder(shape, A, D*M_PI/180); + BRepPrimAPI_MakeRevol builder(shape, A, D*M_PI/180, true); for (auto typ : { TopAbs_EDGE, TopAbs_VERTEX }) for (TopExp_Explorer e(shape, typ); e.More(); e.Next()) @@ -1229,30 +1096,17 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) BRepMesh_IncrementalMesh (shape, deflection, true); }) - .def("Identify", [](const TopoDS_Shape & me, const TopoDS_Shape & you, string name) { - // only edges supported, by now - auto me_edge = TopoDS::Edge(me); - auto you_edge = TopoDS::Edge(you); - 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") + .def("Identify", py::overload_cast>(&Identify), + py::arg("other"), py::arg("name"), + py::arg("type")=Identifications::PERIODIC, py::arg("trafo")=nullopt, + "Identify shapes for periodic meshing") + .def("Distance", [](const TopoDS_Shape& self, + const TopoDS_Shape& other) + { + return BRepExtrema_DistShapeShape(self, other).Value(); + }) .def("Triangulation", [](const TopoDS_Shape & shape) { @@ -1582,6 +1436,14 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) auto ax = gp_Ax3(p, du^dv, du); return make_shared (ax); }) + .def("ProjectWire", [](const TopoDS_Face& face, + const TopoDS_Wire& wire) + { + BRepAlgo_NormalProjection builder(face); + builder.Add(wire); + builder.Build(); + return builder.Projection(); + }) ; py::class_ (m, "Solid"); @@ -1685,6 +1547,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) }) .def_property_readonly("solids", &ListOfShapes::Solids) .def_property_readonly("faces", &ListOfShapes::Faces) + .def_property_readonly("wires", &ListOfShapes::Wires) .def_property_readonly("edges", &ListOfShapes::Edges) .def_property_readonly("vertices", &ListOfShapes::Vertices) .def(py::self * py::self) @@ -1731,12 +1594,15 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) .def("Nearest", [] (ListOfShapes & shapes, gp_Pnt pnt) { return CastShape(shapes.Nearest(pnt)); }, py::arg("p"), "returns shape nearest to point 'p'") + .def("Nearest", [] (ListOfShapes & shapes, gp_Pnt2d pnt) + { return CastShape(shapes.Nearest( { pnt.X(), pnt.Y(), 0 })); }, + py::arg("p"), "returns shape nearest to point 'p'") .def_property("name", [](ListOfShapes& shapes) { throw Exception("Cannot get property of ListOfShapes, get the property from individual shapes!"); }, - [](ListOfShapes& shapes, std::string name) + [](ListOfShapes& shapes, optional name) { for(auto& shape : shapes) { @@ -1764,7 +1630,24 @@ 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") + .def("Identify", py::overload_cast(&Identify), + py::arg("other"), py::arg("name"), + py::arg("type")=Identifications::PERIODIC, py::arg("trafo"), + "Identify shapes for periodic meshing") + ; @@ -1873,13 +1756,13 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) "create box with opposite points 'p1' and 'p2'"); m.def("Prism", [] (const TopoDS_Shape & face, gp_Vec vec) { - return BRepPrimAPI_MakePrism (face, vec).Shape(); + return BRepPrimAPI_MakePrism (face, vec, true).Shape(); }, py::arg("face"), py::arg("v"), "extrude face along the vector 'v'"); m.def("Revolve", [] (const TopoDS_Shape & face,const gp_Ax1 &A, const double D) { - //comvert angle from deg to rad - return BRepPrimAPI_MakeRevol (face, A, D*M_PI/180).Shape(); + //convert angle from deg to rad + return BRepPrimAPI_MakeRevol (face, A, D*M_PI/180, true).Shape(); }); m.def("Pipe", [] (const TopoDS_Wire & spine, const TopoDS_Shape & profile, diff --git a/libsrc/occ/vsocc.cpp b/libsrc/occ/vsocc.cpp index 3f11f47d..9658eb77 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 @@ -87,7 +86,7 @@ namespace netgen // Added clipping planes to Geometry view SetClippingPlane(); - GLfloat matcoledge[] = { 0, 0, 1, 1}; + GLfloat matcoledge[] = { 0, 0, 0, 1}; GLfloat matcolhiedge[] = { 1, 0, 0, 1}; glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcoledge); diff --git a/nglib/CMakeLists.txt b/nglib/CMakeLists.txt index 6e296304..a8b4bac5 100644 --- a/nglib/CMakeLists.txt +++ b/nglib/CMakeLists.txt @@ -9,7 +9,6 @@ if(WIN32) $ $ - $ ) if(USE_GUI) set(nglib_objects ${nglib_objects} @@ -18,6 +17,9 @@ if(WIN32) $ ) endif(USE_GUI) + if(USE_OCC) + set(nglib_objects ${nglib_objects} $) + endif(USE_OCC) endif(WIN32) add_library(nglib SHARED nglib.cpp ${nglib_objects}) diff --git a/nglib/nglib.cpp b/nglib/nglib.cpp index 9eb1e8c3..9be8ddbe 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())) { @@ -910,11 +910,6 @@ namespace nglib // parameters structure mp->Transfer_Parameters(); - - // Only go into surface meshing if the face descriptors have already been added - if(!me->GetNFD()) - return NG_ERROR; - numpoints = me->GetNP(); // Initially set up only for surface meshing without any optimisation @@ -926,8 +921,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 diff --git a/tests/pytest/compare_results.py b/tests/pytest/compare_results.py index 84e9efe8..7f5b6b2b 100644 --- a/tests/pytest/compare_results.py +++ b/tests/pytest/compare_results.py @@ -73,6 +73,13 @@ for bad1,bad2, f1, f2 in zip(data['badness'], data2['badness'], data['file'], da if bad2>0 and bad2<0.9*bad1: print(f"file {f1} got better: {bad1} -> {bad2}") +for bad1,bad2, f1, f2 in zip(data['#trigs'], data2['#trigs'], data['file'], data2['file']): + assert f1==f2 + if bad2>0 and bad2>1.1*bad1: + print(f"file {f1} got worse: {bad1} -> {bad2}") + if bad2>0 and bad2<0.9*bad1: + print(f"file {f1} got better: {bad1} -> {bad2}") + n = len(data)+1 fig,ax = plt.subplots(figsize=(10,7)) for i,d in enumerate(['min trig angle','min tet angle','max trig angle','max tet angle']): diff --git a/tests/pytest/results.json b/tests/pytest/results.json index 595c6a87..5f079143 100644 --- a/tests/pytest/results.json +++ b/tests/pytest/results.json @@ -1590,8 +1590,8 @@ 0.0, 0.0 ], - "ne1d": 81, - "ne2d": 436, + "ne1d": 80, + "ne2d": 412, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -1606,7 +1606,7 @@ 0.0 ], "ne1d": 86, - "ne2d": 452, + "ne2d": 440, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -1621,7 +1621,7 @@ 0.0 ], "ne1d": 86, - "ne2d": 450, + "ne2d": 464, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -1635,8 +1635,8 @@ 0.0, 0.0 ], - "ne1d": 81, - "ne2d": 436, + "ne1d": 80, + "ne2d": 412, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -1650,8 +1650,8 @@ 0.0, 0.0 ], - "ne1d": 82, - "ne2d": 441, + "ne1d": 83, + "ne2d": 453, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -1665,8 +1665,8 @@ 0.0, 0.0 ], - "ne1d": 86, - "ne2d": 472, + "ne1d": 84, + "ne2d": 462, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -2790,8 +2790,38 @@ 0.0, 0.0 ], - "ne1d": 31, - "ne2d": 97, + "ne1d": 27, + "ne2d": 87, + "ne3d": 0, + "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", + "total_badness": 0.0 + }, + { + "angles_tet": [ + 0.0, + 0.0 + ], + "angles_trig": [ + 0.0, + 0.0 + ], + "ne1d": 28, + "ne2d": 80, + "ne3d": 0, + "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", + "total_badness": 0.0 + }, + { + "angles_tet": [ + 0.0, + 0.0 + ], + "angles_trig": [ + 0.0, + 0.0 + ], + "ne1d": 26, + "ne2d": 78, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -2806,7 +2836,7 @@ 0.0 ], "ne1d": 27, - "ne2d": 71, + "ne2d": 87, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -2820,8 +2850,8 @@ 0.0, 0.0 ], - "ne1d": 25, - "ne2d": 85, + "ne1d": 36, + "ne2d": 150, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -2835,38 +2865,8 @@ 0.0, 0.0 ], - "ne1d": 31, - "ne2d": 97, - "ne3d": 0, - "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", - "total_badness": 0.0 - }, - { - "angles_tet": [ - 0.0, - 0.0 - ], - "angles_trig": [ - 0.0, - 0.0 - ], - "ne1d": 41, - "ne2d": 173, - "ne3d": 0, - "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", - "total_badness": 0.0 - }, - { - "angles_tet": [ - 0.0, - 0.0 - ], - "angles_trig": [ - 0.0, - 0.0 - ], - "ne1d": 30, - "ne2d": 108, + "ne1d": 42, + "ne2d": 246, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -2883,7 +2883,7 @@ 0.0 ], "ne1d": 32, - "ne2d": 148, + "ne2d": 142, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -2913,7 +2913,7 @@ 0.0 ], "ne1d": 32, - "ne2d": 134, + "ne2d": 126, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -2928,7 +2928,7 @@ 0.0 ], "ne1d": 32, - "ne2d": 148, + "ne2d": 142, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -2942,8 +2942,8 @@ 0.0, 0.0 ], - "ne1d": 43, - "ne2d": 300, + "ne1d": 42, + "ne2d": 326, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -2957,8 +2957,8 @@ 0.0, 0.0 ], - "ne1d": 79, - "ne2d": 821, + "ne1d": 76, + "ne2d": 811, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -2975,7 +2975,7 @@ 0.0 ], "ne1d": 32, - "ne2d": 124, + "ne2d": 118, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -3005,7 +3005,7 @@ 0.0 ], "ne1d": 32, - "ne2d": 110, + "ne2d": 102, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -3020,7 +3020,7 @@ 0.0 ], "ne1d": 32, - "ne2d": 124, + "ne2d": 118, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -3034,8 +3034,8 @@ 0.0, 0.0 ], - "ne1d": 43, - "ne2d": 249, + "ne1d": 42, + "ne2d": 274, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 @@ -3049,8 +3049,8 @@ 0.0, 0.0 ], - "ne1d": 79, - "ne2d": 687, + "ne1d": 76, + "ne2d": 672, "ne3d": 0, "quality_histogram": "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]", "total_badness": 0.0 diff --git a/tests/pytest/test_occ.py b/tests/pytest/test_occ.py new file mode 100644 index 00000000..510775cb --- /dev/null +++ b/tests/pytest/test_occ.py @@ -0,0 +1,43 @@ +import pytest +import math +from pytest import approx +from pytest_check import check + +def check_volume(shape, volume, dim=3): + from netgen.occ import OCCGeometry + geo = OCCGeometry(shape, dim=dim) + + m = geo.GenerateMesh() + assert len(m.Elements2D()) > 0 + assert len(m.Elements1D()) > 0 + if dim==3: + assert len(m.Elements3D()) > 0 + + ngs = pytest.importorskip("ngsolve") + mesh = ngs.Mesh(m) + mesh.Curve(5) + with check: assert ngs.Integrate(1.0, mesh) == approx(volume) + +def test_rect_with_two_holes(): + occ = pytest.importorskip("netgen.occ") + + face = occ.WorkPlane().Rectangle(7,4) \ + .Circle(2,2,1).Reverse() \ + .Circle(5,2,1).Reverse().Face() + check_volume(face, 7*4-2*math.pi, 2) + +def test_unit_square(): + occ = pytest.importorskip("netgen.occ") + check_volume(occ.unit_square.shape, 1, dim=2) + +def test_box_and_cyl(): + occ = pytest.importorskip("netgen.occ") + box = occ.Box(occ.Pnt(0,0,0), occ.Pnt(1,1,1)) + check_volume(box, 1) + r = 0.3 + h = 0.5 + vcyl = r*r*math.pi*h + cyl = occ.Cylinder(occ.Pnt(1,0.5,0.5), occ.X, r=r, h=h) + check_volume(cyl, vcyl) + fused = box+cyl + check_volume(fused, 1+vcyl) diff --git a/tests/pytest/test_occ_identifications.py b/tests/pytest/test_occ_identifications.py new file mode 100644 index 00000000..d15abc5e --- /dev/null +++ b/tests/pytest/test_occ_identifications.py @@ -0,0 +1,76 @@ +import pytest + +from netgen.meshing import IdentificationType +idtype = IdentificationType.CLOSESURFACES + +def test_two_boxes(): + occ = pytest.importorskip("netgen.occ") + inner = occ.Box((0,0,0), (1,1,1)) + trafo = occ.gp_Trsf().Scale(inner.center, 1.1) + outer = trafo(inner) + + inner.Identify(outer, "", idtype, trafo) + shape = occ.Glue([outer-inner, inner]) + + geo = occ.OCCGeometry(shape) + mesh = geo.GenerateMesh(maxh=0.3) + have_prisms = False + + for el in mesh.Elements3D(): + if len(el.vertices)==6: + have_prisms = True + break + + assert have_prisms + +def test_two_circles(): + occ = pytest.importorskip("netgen.occ") + circ1 = occ.WorkPlane().Circle(1).Face() + trafo = occ.gp_Trsf().Scale(circ1.center, 1.1) + + circ2 = trafo(circ1) + circ1.edges[0].Identify(circ2.edges[0], "", idtype, trafo) + circ2 -= circ1 + shape = occ.Glue([circ1, circ2]) + + geo = occ.OCCGeometry(shape, 2) + mesh = geo.GenerateMesh(maxh=0.2) + have_quads = False + + for el in mesh.Elements2D(): + if len(el.vertices)==4: + have_quads = True + break + + assert have_quads + +def test_cut_identified_face(): + occ = pytest.importorskip("netgen.occ") + from netgen.occ import Z, Box, Cylinder, Glue, OCCGeometry + box = Box((-1,-1,0), (1,1,1)) + cyl = Cylinder( (0,0,0), Z, 0.5, 1 ) + + box.faces.Min(Z).Identify(box.faces.Max(Z), "", idtype) + shape = Glue([cyl, box]) + geo = OCCGeometry(shape) + mesh = geo.GenerateMesh(maxh=0.5) + + for el in mesh.Elements3D(): + assert len(el.vertices)==6 + +def test_identify_multiple_faces(): + occ = pytest.importorskip("netgen.occ") + from netgen.occ import Z, Box, Cylinder, Glue, OCCGeometry, gp_Trsf + box = Box((-1,-1,0), (1,1,1)) + cyl = Cylinder( (0,0,0), Z, 0.5, 1 ) + + shape = Glue([box, cyl]) + bot_faces = shape.faces[Z < 0.1] + top_faces = shape.faces[Z > 0.1] + bot_faces.Identify(top_faces, "", idtype, gp_Trsf.Translation((0,0,1))) + + geo = OCCGeometry(shape) + mesh = geo.GenerateMesh(maxh=0.3) + + for el in mesh.Elements3D(): + assert len(el.vertices)==6