From 571cbbe4df9f90453040a080b0477e97038f58d9 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Mon, 3 Jun 2024 12:37:26 +0200 Subject: [PATCH] Optional identification_name argument in Face::Offset to apply CLOSE_SURFACE identifications --- libsrc/meshing/basegeom.cpp | 189 +++++++++++++++++++------------ libsrc/meshing/basegeom.hpp | 4 +- libsrc/occ/occ_utils.hpp | 7 +- libsrc/occ/occgeom.cpp | 13 ++- libsrc/occ/occgeom.hpp | 6 +- libsrc/occ/python_occ_shapes.cpp | 33 +++--- 6 files changed, 152 insertions(+), 100 deletions(-) diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp index cab14535..119986ec 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -266,7 +266,7 @@ namespace netgen 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)) + if(ident.trafo && e->IsMappedShape(*e_other, *ident.trafo, tol)) e->identifications.Append( {e, e_other, ident.trafo, ident.type, ident.name} ); for(auto & e : edges) @@ -278,9 +278,11 @@ namespace netgen GeometryVertex * pfrom[] = { &from.GetStartVertex(), &from.GetEndVertex() }; GeometryVertex * pto[] = { &to.GetStartVertex(), &to.GetEndVertex() }; + if(!ident.trafo) continue; + // 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_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)) @@ -298,10 +300,7 @@ namespace netgen 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; @@ -315,12 +314,19 @@ namespace netgen auto other = need_inverse ? ident.to : ident.from; 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; + if(ident.trafo) + { + auto trafo = *ident.trafo; + if(need_inverse) + trafo = trafo.CalcInverse(); + if(!s->primary_to_me) + s->primary_to_me = Transformation<3>( Vec<3>{0., 0., 0.} ); + if(!other->primary_to_me) + other->primary_to_me = Transformation<3>( Vec<3>{0., 0., 0.} ); + s->primary_to_me->Combine(trafo, *other->primary_to_me); + changed = true; + } } } } @@ -658,7 +664,9 @@ namespace netgen edge_params.SetSize(np-2); for(auto i : Range(np-2)) { - edge_points[i] = trafo(mesh[pnums_primary[i+1]]); + edge_points[i] = mesh[pnums_primary[i+1]]; + if(trafo) + edge_points[i] = (*trafo)(edge_points[i]); EdgePointGeomInfo gi; edge->ProjectPoint(edge_points[i], &gi); edge_params[i] = gi.dist; @@ -689,7 +697,8 @@ namespace netgen { for(size_t i : std::vector{0UL, pnums_primary.Size()-1}) { - auto p_mapped = trafo(mesh[pnums_primary[i]]); + auto p_mapped = mesh[pnums_primary[i]]; + if(trafo) p_mapped = (*trafo)(p_mapped); EdgePointGeomInfo gi; edge->ProjectPoint(p_mapped, &gi); params[i] = gi.dist; @@ -739,10 +748,16 @@ namespace netgen if(ident.from == edge.get()) { auto & pnums = all_pnums[edge->nr]; + if(pnums.Size() < 2) continue; // degenerated edge // 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])); + Point<3> p_other = mesh[pi]; + if(ident.trafo) + p_other = (*ident.trafo)(mesh[pi]); + else + static_cast(ident.to)->ProjectPoint(p_other, nullptr); + auto pi_other = tree.Find(p_other); identifications.Add(pi, pi_other, ident.name, ident.type); } } @@ -860,7 +875,7 @@ namespace netgen constexpr int NOT_MAPPED = -1; mapped_edges = UNINITIALIZED; - Transformation<3> trafo; + optional> trafo; if(face.IsConnectingCloseSurfaces()) { @@ -902,8 +917,6 @@ namespace netgen 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]])); auto gis = sel.GeomInfo(); for(auto i : Range(2)) { @@ -911,6 +924,21 @@ namespace netgen gis[i].v = s.epgeominfo[i].v; } + Point<3> p2 = mesh[s[1]]; + Point<3> p3 = mesh[s[0]]; + if(trafo) + { + p2 = (*trafo)(p2); + p3 = (*trafo)(p3); + } + else + { + edges[mapped_edges[edgenr]]->ProjectPoint(p2, nullptr); + edges[mapped_edges[edgenr]]->ProjectPoint(p3, nullptr); + } + sel[2] = tree.Find(p2); + sel[3] = tree.Find(p3); + // find mapped segment to set PointGeomInfo correctly Segment s_other; for(auto si_other : p2seg[sel[2]]) @@ -1040,7 +1068,12 @@ namespace netgen auto pi = seg[i]; if(!is_point_in_tree[pi]) { - tree.Insert(trafo(mesh[pi]), pi); + auto p = mesh[pi]; + if(trafo) + p = (*trafo)(p); + else + dst.Project(p); + tree.Insert(p, pi); is_point_in_tree[pi] = true; } } @@ -1068,6 +1101,7 @@ namespace netgen } xbool do_invert = maybe; + if(!trafo) do_invert = true; // now insert mapped surface elements for(auto sei : mesh.SurfaceElements().Range()) @@ -1076,14 +1110,6 @@ namespace netgen if(sel.GetIndex() != src.nr+1) continue; - if(do_invert.IsMaybe()) - { - auto n_src = src.GetNormal(mesh[sel[0]]); - auto n_dist = dst.GetNormal(trafo(mesh[sel[0]])); - Mat<3> normal_matrix; - CalcInverse(Trans(trafo.GetMatrix()), normal_matrix); - do_invert = (normal_matrix * n_src) * n_dist < 0.0; - } auto sel_new = sel; sel_new.SetIndex(dst.nr+1); for(auto i : Range(sel.PNums())) @@ -1091,62 +1117,75 @@ namespace netgen auto pi = sel[i]; if(!pmap[pi].IsValid()) { - pmap[pi] = mesh.AddPoint(trafo(mesh[pi]), 1, SURFACEPOINT); + auto p = mesh[pi]; + if(trafo) + p = (*trafo)(p); + else + dst.Project(p); + pmap[pi] = mesh.AddPoint(p, 1, SURFACEPOINT); } sel_new[i] = pmap[pi]; mapto[{pi, dst.nr}] = pmap[pi]; mapto[{pmap[pi], src.nr}] = pi; } - if(do_invert.IsTrue()) - sel_new.Invert(); + if(do_invert.IsMaybe()) + { + auto n_src = src.GetNormal(mesh[sel[0]]); + auto n_dist = dst.GetNormal(mesh[sel_new[0]]); + Mat<3> normal_matrix; + CalcInverse(Trans(trafo->GetMatrix()), normal_matrix); + do_invert = (normal_matrix * n_src) * n_dist < 0.0; + } + if(do_invert.IsTrue()) + sel_new.Invert(); - for(auto i : Range(sel.PNums())) - { - auto pi = sel_new[i]; - if(uv_values.Range().Next() <= pi) - { - // new point (inner surface point) - PointGeomInfo gi; - dst.CalcPointGeomInfo(mesh[sel_new[i]], gi); - sel_new.GeomInfo()[i] = gi; - continue; - } + for(auto i : Range(sel.PNums())) + { + auto pi = sel_new[i]; + if(uv_values.Range().Next() <= pi) + { + // new point (inner surface point) + PointGeomInfo gi; + dst.CalcPointGeomInfo(mesh[sel_new[i]], gi); + sel_new.GeomInfo()[i] = gi; + continue; + } - const auto & uvs = uv_values[pi]; - if(uvs.Size() == 1) - { - // appears only once -> take uv values from edgepointgeominfo - const auto & [u,v] = uvs[0]; - PointGeomInfo gi; - gi.u = u; - gi.v = v; - sel_new.GeomInfo()[i] = gi; - } - else if(uvs.Size() > 1) - { - // multiple uv pairs -> project to close point and select closest uv pair - double eps = 1e-3; - auto p = Point<3>((1.0-eps)*Vec<3>(mesh[sel_new.PNumMod(i+1)]) + - eps/2*Vec<3>(mesh[sel_new.PNumMod(i+2)]) + - eps/2*Vec<3>(mesh[sel_new.PNumMod(i+3)])); - PointGeomInfo gi_p, gi; - dst.CalcPointGeomInfo(p, gi_p); - gi.trignum = gi_p.trignum; - double min_dist = numeric_limits::max(); - for(const auto & [u,v] : uvs) - { - double dist = (gi_p.u-u)*(gi_p.u-u) + (gi_p.v-v)*(gi_p.v-v); - if(dist < min_dist) - { - min_dist = dist; - gi.u = u; - gi.v = v; - } - } - sel_new.GeomInfo()[i] = gi; - } - else - throw Exception(string(__FILE__) + ":"+ToString(__LINE__) + " shouldn't come here"); + const auto & uvs = uv_values[pi]; + if(uvs.Size() == 1) + { + // appears only once -> take uv values from edgepointgeominfo + const auto & [u,v] = uvs[0]; + PointGeomInfo gi; + gi.u = u; + gi.v = v; + sel_new.GeomInfo()[i] = gi; + } + else if(uvs.Size() > 1) + { + // multiple uv pairs -> project to close point and select closest uv pair + double eps = 1e-3; + auto p = Point<3>((1.0-eps)*Vec<3>(mesh[sel_new.PNumMod(i+1)]) + + eps/2*Vec<3>(mesh[sel_new.PNumMod(i+2)]) + + eps/2*Vec<3>(mesh[sel_new.PNumMod(i+3)])); + PointGeomInfo gi_p, gi; + dst.CalcPointGeomInfo(p, gi_p); + gi.trignum = gi_p.trignum; + double min_dist = numeric_limits::max(); + for(const auto & [u,v] : uvs) + { + double dist = (gi_p.u-u)*(gi_p.u-u) + (gi_p.v-v)*(gi_p.v-v); + if(dist < min_dist) + { + min_dist = dist; + gi.u = u; + gi.v = v; + } + } + sel_new.GeomInfo()[i] = gi; + } + else + throw Exception(string(__FILE__) + ":"+ToString(__LINE__) + " shouldn't come here"); } mesh.AddSurfaceElement(sel_new); } diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index f625baec..8d5a0ee5 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -54,7 +54,7 @@ namespace netgen { GeometryShape * from; GeometryShape * to; - Transformation<3> trafo; + optional> trafo; Identifications::ID_TYPE type; string name = ""; }; @@ -67,7 +67,7 @@ namespace netgen ShapeProperties properties; Array identifications; GeometryShape * primary; - Transformation<3> primary_to_me; + optional> primary_to_me = nullopt; virtual ~GeometryShape() {} virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const; diff --git a/libsrc/occ/occ_utils.hpp b/libsrc/occ/occ_utils.hpp index 36a2f3d9..8696cb55 100644 --- a/libsrc/occ/occ_utils.hpp +++ b/libsrc/occ/occ_utils.hpp @@ -65,15 +65,14 @@ namespace netgen DLL_HEADER Box<3> GetBoundingBox( const TopoDS_Shape & shape ); - class OCCIdentification + struct OCCIdentification { - public: TopoDS_Shape from; TopoDS_Shape to; - Transformation<3> trafo; + optional> trafo = nullopt; string name; Identifications::ID_TYPE type; - bool opposite_direction; + bool opposite_direction = false; }; Standard_Integer BuildTriangulation( const TopoDS_Shape & shape ); diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index 821cbf8b..bfc96cb0 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -2365,10 +2365,12 @@ namespace netgen Array items; items.Append(MakeReal(ident.from == shape ? 1 : 0)); items.Append(to); - auto & m = ident.trafo.GetMatrix(); + Transformation<3> trafo; + if(ident.trafo) trafo = *ident.trafo; + auto & m = trafo.GetMatrix(); for(auto i : Range(9)) items.Append(MakeReal(m(i))); - auto & v = ident.trafo.GetVector(); + auto & v = trafo.GetVector(); for(auto i : Range(3)) items.Append(MakeReal(v(i))); items.Append(MakeInt(ident.type)); @@ -2407,12 +2409,15 @@ namespace netgen ident.to = shape_origin; } - auto & m = ident.trafo.GetMatrix(); + Transformation<3> trafo; + auto & m = trafo.GetMatrix(); for(auto i : Range(9)) m(i) = ReadReal(id_item->ItemElementValue(3+i)); - auto & v = ident.trafo.GetVector(); + auto & v = trafo.GetVector(); for(auto i : Range(3)) v(i) = ReadReal(id_item->ItemElementValue(12+i)); + if(FlatVector(9, &trafo.GetMatrix()(0,0)).L2Norm() != .0 && trafo.GetVector().Length2() != .0) + ident.trafo = trafo; ident.type = Identifications::ID_TYPE(ReadInt(id_item->ItemElementValue(15))); result.push_back(ident); } diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index 22bd1428..fa00163a 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -516,11 +516,13 @@ namespace netgen if(from.IsSame(from_mapped) && to.IsSame(to_mapped)) continue; - Transformation<3> trafo_mapped = ident.trafo; + if(!ident.trafo) continue; + Transformation<3> trafo_mapped = *ident.trafo; + if(trafo) { Transformation<3> trafo_temp; - trafo_temp.Combine(ident.trafo, trafo_inv); + trafo_temp.Combine(*ident.trafo, trafo_inv); trafo_mapped.Combine(*trafo, trafo_temp); } diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index ed92b955..8d728873 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -1152,46 +1152,53 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) .def("Offset", [](const TopoDS_Shape & shape, double offset, double tol, bool intersection, - string joinT, bool removeIntEdges) { + string joinT, bool removeIntEdges, optional identification_name) { BRepOffsetAPI_MakeOffsetShape maker; GeomAbs_JoinType joinType; if(joinT == "arc") joinType = GeomAbs_Arc; else if(joinT == "intersection") joinType = GeomAbs_Intersection; + else if(joinT == "tangent") + joinType = GeomAbs_Tangent; else - throw Exception("Only joinTypes 'arc' and 'intersection' exist!"); + throw Exception("Only joinTypes 'arc', 'intersection' and 'tangent' exist!"); maker.PerformByJoin(shape, offset, tol, BRepOffset_Skin, intersection, false, joinType, removeIntEdges); // PropagateProperties (maker, shape); - - // bool have_identifications = false; for (auto typ : { TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX }) for (TopExp_Explorer e(shape, typ); e.More(); e.Next()) { auto s = e.Current(); - // have_identifications |= OCCGeometry::HaveIdentifications(s); - if(!OCCGeometry::HaveProperties(s)) - continue; auto prop = OCCGeometry::GetProperties(s); for (auto mods : maker.Generated(s)) { - // OCCGeometry::GetProperties(mods).Merge(prop); - auto & props = OCCGeometry::GetProperties(mods); - props.Merge(prop); - if (prop.name) props.name = string("offset_")+(*prop.name); + if(OCCGeometry::HaveProperties(s)) + { + auto & new_props = OCCGeometry::GetProperties(mods); + new_props.Merge(prop); + if (prop.name) new_props.name = string("offset_")+(*prop.name); + } + if(identification_name) + { + OCCIdentification ident; + ident.from = s; + ident.to = mods; + ident.name = *identification_name; + ident.type = Identifications::CLOSESURFACES; + OCCGeometry::GetIdentifications(s).push_back(ident); + } } } - // if(have_identifications) - // PropagateIdentifications(builder, shape, trafo); return maker.Shape(); }, py::arg("offset"), py::arg("tol"), py::arg("intersection") = false,py::arg("joinType")="arc", py::arg("removeIntersectingEdges") = false, + py::arg("identification_name") = nullopt, "makes shell-like solid from faces")