ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-24 21:10:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

Optional identification_name argument in Face::Offset to apply CLOSE_SURFACE identifications

Browse Source
This commit is contained in:
Matthias Hochsteger 2024-06-03 12:37:26 +02:00
parent 6d1c87f214
commit 571cbbe4df
6 changed files with 152 additions and 100 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

189
libsrc/meshing/basegeom.cpp
View File

@ -266,7 +266,7 @@ namespace netgen
for(auto & ident: f->identifications) for(auto & ident: f->identifications)
for(auto e : static_cast<GeometryFace*>(ident.from)->edges) for(auto e : static_cast<GeometryFace*>(ident.from)->edges)
for(auto e_other : static_cast<GeometryFace*>(ident.to)->edges) for(auto e_other : static_cast<GeometryFace*>(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} ); e->identifications.Append( {e, e_other, ident.trafo, ident.type, ident.name} );
for(auto & e : edges) for(auto & e : edges)
@ -278,9 +278,11 @@ namespace netgen
GeometryVertex * pfrom[] = { &from.GetStartVertex(), &from.GetEndVertex() }; GeometryVertex * pfrom[] = { &from.GetStartVertex(), &from.GetEndVertex() };
GeometryVertex * pto[] = { &to.GetStartVertex(), &to.GetEndVertex() }; GeometryVertex * pto[] = { &to.GetStartVertex(), &to.GetEndVertex() };
if(!ident.trafo) continue;
// swap points of other edge if necessary // swap points of other edge if necessary
Point<3> p_from0 = ident.trafo(from.GetStartVertex().GetPoint()); Point<3> p_from0 = (*ident.trafo)(from.GetStartVertex().GetPoint());
Point<3> p_from1 = ident.trafo(from.GetEndVertex().GetPoint()); Point<3> p_from1 = (*ident.trafo)(from.GetEndVertex().GetPoint());
Point<3> p_to0 = to.GetStartVertex().GetPoint(); Point<3> p_to0 = to.GetStartVertex().GetPoint();
if(Dist(p_from1, p_to0) < Dist(p_from0, p_to0)) if(Dist(p_from1, p_to0) < Dist(p_from0, p_to0))
@ -298,10 +300,7 @@ namespace netgen
auto find_primary = [&] (auto & shapes) auto find_primary = [&] (auto & shapes)
{ {
for(auto &s : shapes) for(auto &s : shapes)
{
s->primary = s.get(); s->primary = s.get();
s->primary_to_me = Transformation<3>{ Vec<3> {0,0,0} }; // init with identity
}
bool changed = true; bool changed = true;
@ -315,12 +314,19 @@ namespace netgen
auto other = need_inverse ? ident.to : ident.from; auto other = need_inverse ? ident.to : ident.from;
if(other->primary->nr < s->primary->nr) if(other->primary->nr < s->primary->nr)
{ {
auto trafo = ident.trafo;
if(need_inverse)
trafo = trafo.CalcInverse();
s->primary = other->primary; s->primary = other->primary;
s->primary_to_me.Combine(trafo, other->primary_to_me); if(ident.trafo)
changed = true; {
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); edge_params.SetSize(np-2);
for(auto i : Range(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; EdgePointGeomInfo gi;
edge->ProjectPoint(edge_points[i], &gi); edge->ProjectPoint(edge_points[i], &gi);
edge_params[i] = gi.dist; edge_params[i] = gi.dist;
@ -689,7 +697,8 @@ namespace netgen
{ {
for(size_t i : std::vector{0UL, pnums_primary.Size()-1}) 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; EdgePointGeomInfo gi;
edge->ProjectPoint(p_mapped, &gi); edge->ProjectPoint(p_mapped, &gi);
params[i] = gi.dist; params[i] = gi.dist;
@ -739,10 +748,16 @@ namespace netgen
if(ident.from == edge.get()) if(ident.from == edge.get())
{ {
auto & pnums = all_pnums[edge->nr]; auto & pnums = all_pnums[edge->nr];
if(pnums.Size() < 2) continue; // degenerated edge
// start and end vertex are already identified // start and end vertex are already identified
for(auto pi : pnums.Range(1, pnums.Size()-1)) 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<GeometryEdge*>(ident.to)->ProjectPoint(p_other, nullptr);
auto pi_other = tree.Find(p_other);
identifications.Add(pi, pi_other, ident.name, ident.type); identifications.Add(pi, pi_other, ident.name, ident.type);
} }
} }
@ -860,7 +875,7 @@ namespace netgen
constexpr int NOT_MAPPED = -1; constexpr int NOT_MAPPED = -1;
mapped_edges = UNINITIALIZED; mapped_edges = UNINITIALIZED;
Transformation<3> trafo; optional<Transformation<3>> trafo;
if(face.IsConnectingCloseSurfaces()) if(face.IsConnectingCloseSurfaces())
{ {
@ -902,8 +917,6 @@ namespace netgen
Element2d sel(4); Element2d sel(4);
sel[0] = s[0]; sel[0] = s[0];
sel[1] = s[1]; sel[1] = s[1];
sel[2] = tree.Find(trafo(mesh[s[1]]));
sel[3] = tree.Find(trafo(mesh[s[0]]));
auto gis = sel.GeomInfo(); auto gis = sel.GeomInfo();
for(auto i : Range(2)) for(auto i : Range(2))
{ {
@ -911,6 +924,21 @@ namespace netgen
gis[i].v = s.epgeominfo[i].v; 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 // find mapped segment to set PointGeomInfo correctly
Segment s_other; Segment s_other;
for(auto si_other : p2seg[sel[2]]) for(auto si_other : p2seg[sel[2]])
@ -1040,7 +1068,12 @@ namespace netgen
auto pi = seg[i]; auto pi = seg[i];
if(!is_point_in_tree[pi]) 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; is_point_in_tree[pi] = true;
} }
} }
@ -1068,6 +1101,7 @@ namespace netgen
} }
xbool do_invert = maybe; xbool do_invert = maybe;
if(!trafo) do_invert = true;
// now insert mapped surface elements // now insert mapped surface elements
for(auto sei : mesh.SurfaceElements().Range()) for(auto sei : mesh.SurfaceElements().Range())
@ -1076,14 +1110,6 @@ namespace netgen
if(sel.GetIndex() != src.nr+1) if(sel.GetIndex() != src.nr+1)
continue; 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; auto sel_new = sel;
sel_new.SetIndex(dst.nr+1); sel_new.SetIndex(dst.nr+1);
for(auto i : Range(sel.PNums())) for(auto i : Range(sel.PNums()))
@ -1091,62 +1117,75 @@ namespace netgen
auto pi = sel[i]; auto pi = sel[i];
if(!pmap[pi].IsValid()) 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]; sel_new[i] = pmap[pi];
mapto[{pi, dst.nr}] = pmap[pi]; mapto[{pi, dst.nr}] = pmap[pi];
mapto[{pmap[pi], src.nr}] = pi; mapto[{pmap[pi], src.nr}] = pi;
} }
if(do_invert.IsTrue()) if(do_invert.IsMaybe())
sel_new.Invert(); {
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())) for(auto i : Range(sel.PNums()))
{ {
auto pi = sel_new[i]; auto pi = sel_new[i];
if(uv_values.Range().Next() <= pi) if(uv_values.Range().Next() <= pi)
{ {
// new point (inner surface point) // new point (inner surface point)
PointGeomInfo gi; PointGeomInfo gi;
dst.CalcPointGeomInfo(mesh[sel_new[i]], gi); dst.CalcPointGeomInfo(mesh[sel_new[i]], gi);
sel_new.GeomInfo()[i] = gi; sel_new.GeomInfo()[i] = gi;
continue; continue;
} }
const auto & uvs = uv_values[pi]; const auto & uvs = uv_values[pi];
if(uvs.Size() == 1) if(uvs.Size() == 1)
{ {
// appears only once -> take uv values from edgepointgeominfo // appears only once -> take uv values from edgepointgeominfo
const auto & [u,v] = uvs[0]; const auto & [u,v] = uvs[0];
PointGeomInfo gi; PointGeomInfo gi;
gi.u = u; gi.u = u;
gi.v = v; gi.v = v;
sel_new.GeomInfo()[i] = gi; sel_new.GeomInfo()[i] = gi;
} }
else if(uvs.Size() > 1) else if(uvs.Size() > 1)
{ {
// multiple uv pairs -> project to close point and select closest uv pair // multiple uv pairs -> project to close point and select closest uv pair
double eps = 1e-3; double eps = 1e-3;
auto p = Point<3>((1.0-eps)*Vec<3>(mesh[sel_new.PNumMod(i+1)]) + 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+2)]) +
eps/2*Vec<3>(mesh[sel_new.PNumMod(i+3)])); eps/2*Vec<3>(mesh[sel_new.PNumMod(i+3)]));
PointGeomInfo gi_p, gi; PointGeomInfo gi_p, gi;
dst.CalcPointGeomInfo(p, gi_p); dst.CalcPointGeomInfo(p, gi_p);
gi.trignum = gi_p.trignum; gi.trignum = gi_p.trignum;
double min_dist = numeric_limits<double>::max(); double min_dist = numeric_limits<double>::max();
for(const auto & [u,v] : uvs) for(const auto & [u,v] : uvs)
{ {
double dist = (gi_p.u-u)*(gi_p.u-u) + (gi_p.v-v)*(gi_p.v-v); double dist = (gi_p.u-u)*(gi_p.u-u) + (gi_p.v-v)*(gi_p.v-v);
if(dist < min_dist) if(dist < min_dist)
{ {
min_dist = dist; min_dist = dist;
gi.u = u; gi.u = u;
gi.v = v; gi.v = v;
} }
} }
sel_new.GeomInfo()[i] = gi; sel_new.GeomInfo()[i] = gi;
} }
else else
throw Exception(string(__FILE__) + ":"+ToString(__LINE__) + " shouldn't come here"); throw Exception(string(__FILE__) + ":"+ToString(__LINE__) + " shouldn't come here");
} }
mesh.AddSurfaceElement(sel_new); mesh.AddSurfaceElement(sel_new);
} }

4
libsrc/meshing/basegeom.hpp
View File

@ -54,7 +54,7 @@ namespace netgen
{ {
GeometryShape * from; GeometryShape * from;
GeometryShape * to; GeometryShape * to;
Transformation<3> trafo; optional<Transformation<3>> trafo;
Identifications::ID_TYPE type; Identifications::ID_TYPE type;
string name = ""; string name = "";
}; };
@ -67,7 +67,7 @@ namespace netgen
ShapeProperties properties; ShapeProperties properties;
Array<ShapeIdentification> identifications; Array<ShapeIdentification> identifications;
GeometryShape * primary; GeometryShape * primary;
Transformation<3> primary_to_me; optional<Transformation<3>> primary_to_me = nullopt;
virtual ~GeometryShape() {} virtual ~GeometryShape() {}
virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const; virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const;

7
libsrc/occ/occ_utils.hpp
View File

@ -65,15 +65,14 @@ namespace netgen
DLL_HEADER Box<3> GetBoundingBox( const TopoDS_Shape & shape ); DLL_HEADER Box<3> GetBoundingBox( const TopoDS_Shape & shape );
class OCCIdentification struct OCCIdentification
{ {
public:
TopoDS_Shape from; TopoDS_Shape from;
TopoDS_Shape to; TopoDS_Shape to;
Transformation<3> trafo; optional<Transformation<3>> trafo = nullopt;
string name; string name;
Identifications::ID_TYPE type; Identifications::ID_TYPE type;
bool opposite_direction; bool opposite_direction = false;
}; };
Standard_Integer BuildTriangulation( const TopoDS_Shape & shape ); Standard_Integer BuildTriangulation( const TopoDS_Shape & shape );

13
libsrc/occ/occgeom.cpp
View File

@ -2365,10 +2365,12 @@ namespace netgen
Array<Handle(StepRepr_RepresentationItem)> items; Array<Handle(StepRepr_RepresentationItem)> items;
items.Append(MakeReal(ident.from == shape ? 1 : 0)); items.Append(MakeReal(ident.from == shape ? 1 : 0));
items.Append(to); 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)) for(auto i : Range(9))
items.Append(MakeReal(m(i))); items.Append(MakeReal(m(i)));
auto & v = ident.trafo.GetVector(); auto & v = trafo.GetVector();
for(auto i : Range(3)) for(auto i : Range(3))
items.Append(MakeReal(v(i))); items.Append(MakeReal(v(i)));
items.Append(MakeInt(ident.type)); items.Append(MakeInt(ident.type));
@ -2407,12 +2409,15 @@ namespace netgen
ident.to = shape_origin; ident.to = shape_origin;
} }
auto & m = ident.trafo.GetMatrix(); Transformation<3> trafo;
auto & m = trafo.GetMatrix();
for(auto i : Range(9)) for(auto i : Range(9))
m(i) = ReadReal(id_item->ItemElementValue(3+i)); m(i) = ReadReal(id_item->ItemElementValue(3+i));
auto & v = ident.trafo.GetVector(); auto & v = trafo.GetVector();
for(auto i : Range(3)) for(auto i : Range(3))
v(i) = ReadReal(id_item->ItemElementValue(12+i)); 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))); ident.type = Identifications::ID_TYPE(ReadInt(id_item->ItemElementValue(15)));
result.push_back(ident); result.push_back(ident);
} }

6
libsrc/occ/occgeom.hpp
View File

@ -516,11 +516,13 @@ namespace netgen
if(from.IsSame(from_mapped) && to.IsSame(to_mapped)) if(from.IsSame(from_mapped) && to.IsSame(to_mapped))
continue; continue;
Transformation<3> trafo_mapped = ident.trafo; if(!ident.trafo) continue;
Transformation<3> trafo_mapped = *ident.trafo;
if(trafo) if(trafo)
{ {
Transformation<3> trafo_temp; Transformation<3> trafo_temp;
trafo_temp.Combine(ident.trafo, trafo_inv); trafo_temp.Combine(*ident.trafo, trafo_inv);
trafo_mapped.Combine(*trafo, trafo_temp); trafo_mapped.Combine(*trafo, trafo_temp);
} }

33
libsrc/occ/python_occ_shapes.cpp
View File

@ -1152,46 +1152,53 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
.def("Offset", [](const TopoDS_Shape & shape, .def("Offset", [](const TopoDS_Shape & shape,
double offset, double tol, bool intersection, double offset, double tol, bool intersection,
string joinT, bool removeIntEdges) { string joinT, bool removeIntEdges, optional<string> identification_name) {
BRepOffsetAPI_MakeOffsetShape maker; BRepOffsetAPI_MakeOffsetShape maker;
GeomAbs_JoinType joinType; GeomAbs_JoinType joinType;
if(joinT == "arc") if(joinT == "arc")
joinType = GeomAbs_Arc; joinType = GeomAbs_Arc;
else if(joinT == "intersection") else if(joinT == "intersection")
joinType = GeomAbs_Intersection; joinType = GeomAbs_Intersection;
else if(joinT == "tangent")
joinType = GeomAbs_Tangent;
else else
throw Exception("Only joinTypes 'arc' and 'intersection' exist!"); throw Exception("Only joinTypes 'arc', 'intersection' and 'tangent' exist!");
maker.PerformByJoin(shape, offset, tol, maker.PerformByJoin(shape, offset, tol,
BRepOffset_Skin, intersection, BRepOffset_Skin, intersection,
false, joinType, removeIntEdges); false, joinType, removeIntEdges);
// PropagateProperties (maker, shape); // PropagateProperties (maker, shape);
// bool have_identifications = false;
for (auto typ : { TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX }) for (auto typ : { TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX })
for (TopExp_Explorer e(shape, typ); e.More(); e.Next()) for (TopExp_Explorer e(shape, typ); e.More(); e.Next())
{ {
auto s = e.Current(); auto s = e.Current();
// have_identifications |= OCCGeometry::HaveIdentifications(s);
if(!OCCGeometry::HaveProperties(s))
continue;
auto prop = OCCGeometry::GetProperties(s); auto prop = OCCGeometry::GetProperties(s);
for (auto mods : maker.Generated(s)) for (auto mods : maker.Generated(s))
{ {
// OCCGeometry::GetProperties(mods).Merge(prop); if(OCCGeometry::HaveProperties(s))
auto & props = OCCGeometry::GetProperties(mods); {
props.Merge(prop); auto & new_props = OCCGeometry::GetProperties(mods);
if (prop.name) props.name = string("offset_")+(*prop.name); 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(); return maker.Shape();
}, py::arg("offset"), py::arg("tol"), }, py::arg("offset"), py::arg("tol"),
py::arg("intersection") = false,py::arg("joinType")="arc", py::arg("intersection") = false,py::arg("joinType")="arc",
py::arg("removeIntersectingEdges") = false, py::arg("removeIntersectingEdges") = false,
py::arg("identification_name") = nullopt,
"makes shell-like solid from faces") "makes shell-like solid from faces")