Merge branch 'identify_occ_offset' into 'master'

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

See merge request ngsolve/netgen!653
This commit is contained in:
Schöberl, Joachim 2024-06-03 16:27:45 +02:00
commit 2c7912e5dc
6 changed files with 152 additions and 100 deletions

View File

@ -266,7 +266,7 @@ namespace netgen
for(auto & ident: f->identifications)
for(auto e : static_cast<GeometryFace*>(ident.from)->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} );
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<GeometryEdge*>(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<Transformation<3>> 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<double>::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<double>::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);
}

View File

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

View File

@ -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<Transformation<3>> trafo = nullopt;
string name;
Identifications::ID_TYPE type;
bool opposite_direction;
bool opposite_direction = false;
};
Standard_Integer BuildTriangulation( const TopoDS_Shape & shape );

View File

@ -2365,10 +2365,12 @@ namespace netgen
Array<Handle(StepRepr_RepresentationItem)> 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);
}

View File

@ -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);
}

View File

@ -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<string> 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")