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

fix PropagateIdentifications after Glue

Browse Source
This commit is contained in:
mhochsteger@cerbsim.com 2021-12-15 20:05:18 +01:00
parent bf261d533f
commit b7448ce2c1
9 changed files with 178 additions and 129 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

82
libsrc/meshing/basegeom.cpp
View File

@ -56,18 +56,18 @@ namespace netgen
return false;
auto & e = *other_ptr;
if(tol < Dist(GetCenter(), e.GetCenter()))
if(tol < Dist(trafo(GetCenter()), e.GetCenter()))
return false;
auto v0 = GetStartVertex().GetPoint();
auto v1 = GetEndVertex().GetPoint();
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 = GetPoint(0.5);
v1 = trafo(GetPoint(0.5));
w1 = other_ptr->GetPoint(0.5);
}
@ -75,6 +75,38 @@ namespace netgen
(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<const GeometryFace*>(&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<bool> 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,
@ -169,6 +201,15 @@ namespace netgen
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))
@ -181,11 +222,39 @@ namespace netgen
}
};
auto tol = 1e-8 * bounding_box.Diam();
for(auto & f : faces)
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))
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<GeometryEdge&>(*ident.from);
auto & to = static_cast<GeometryEdge&>(*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 = ident.trafo(to.GetStartVertex().GetPoint());
if(Dist(p_from1, p_to0) < Dist(p_from1, 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);
// todo: propagate identifications faces -> edges -> vertices
auto find_primary = [&] (auto & shapes)
{
@ -729,6 +798,9 @@ namespace netgen
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);
}

16
libsrc/meshing/basegeom.hpp
View File

@ -70,11 +70,19 @@ namespace netgen
class DLL_HEADER GeometryEdge : public GeometryShape
{
protected:
GeometryVertex *start, *end;
public:
int domin=-1, domout=-1;
virtual const GeometryVertex& GetStartVertex() const = 0;
virtual const GeometryVertex& GetEndVertex() const = 0;
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;
@ -100,11 +108,13 @@ namespace netgen
class DLL_HEADER GeometryFace : public GeometryShape
{
public:
Array<GeometryEdge*> edges;
int domin=-1, domout=-1;
virtual Point<3> GetCenter() const = 0;
virtual size_t GetNBoundaries() const = 0;
virtual Array<Segment> 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.
@ -143,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,

1
libsrc/meshing/meshfunc.cpp
View File

@ -210,6 +210,7 @@ namespace netgen
el[i] = sel[i];
el[i+np] = map[sel[i]];
}
el.SetIndex(md.domain);
mesh->AddVolumeElement(el);
}

22
libsrc/occ/occ_edge.cpp
View File

@ -7,8 +7,9 @@
namespace netgen
{
OCCEdge::OCCEdge(TopoDS_Shape edge_)
: tedge(edge_.TShape()),
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);
@ -18,26 +19,13 @@ namespace netgen
if(verts.size() != 2)
throw Exception("OCC edge does not have 2 vertices");
start = OCCVertex(verts[0]);
end = OCCVertex(verts[1]);
// swap start/end if necessary
double d00 = Dist(GetPoint(0), start.GetPoint());
double d01 = Dist(GetPoint(0), end.GetPoint());
double d00 = Dist(GetPoint(0), start->GetPoint());
double d01 = Dist(GetPoint(0), end->GetPoint());
if(d01 < d00)
swap(start, end);
}
const GeometryVertex& OCCEdge::GetStartVertex() const
{
return start;
}
const GeometryVertex& OCCEdge::GetEndVertex() const
{
return end;
}
double OCCEdge::GetLength() const
{
return props.Mass();

8
libsrc/occ/occ_edge.hpp
View File

@ -14,6 +14,7 @@ namespace netgen
{
class OCCEdge : public GeometryEdge
{
public:
T_Shape tedge;
TopoDS_Edge edge;
Handle(Geom_Curve) curve;
@ -21,16 +22,11 @@ namespace netgen
GProp_GProps props;
public:
OCCVertex start;
OCCVertex end;
OCCEdge(TopoDS_Shape edge_);
OCCEdge(TopoDS_Shape edge_, GeometryVertex & start_, GeometryVertex & end_);
auto Shape() const { return edge; }
T_Shape TShape() const { return tedge; }
const GeometryVertex& GetStartVertex() const override;
const GeometryVertex& GetEndVertex() const override;
double GetLength() const override;
Point<3> GetCenter() const override;
Point<3> GetPoint(double t) const override;

1
libsrc/occ/occ_face.cpp
View File

@ -98,7 +98,6 @@ namespace netgen
// auto cof = curve_on_face[ORIENTATION][edgenr];
auto edge = edge_on_face[ORIENTATION][edgenr];
OCCEdge gedge(edge);
double s0, s1;
auto cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1);

98
libsrc/occ/occgeom.cpp
View File

@ -1113,8 +1113,10 @@ namespace netgen
auto tshape = v.TShape();
if(vertex_map.count(tshape)!=0)
continue;
vertex_map[tshape] = vertices.Size();
auto occ_vertex = make_unique<OCCVertex>(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));
@ -1129,10 +1131,9 @@ namespace netgen
if(BRep_Tool::Degenerated(edge))
continue;
edge_map[tshape] = edges.Size();
auto occ_edge = make_unique<OCCEdge>(edge);
auto verts = GetVertices(e);
auto occ_edge = make_unique<OCCEdge>(edge, *vertices[vertex_map[verts[0].TShape()]], *vertices[vertex_map[verts[1].TShape()]] );
occ_edge->properties = global_shape_properties[tshape];
occ_edge->start.nr = vertex_map[occ_edge->start.TShape()];
occ_edge->end.nr = vertex_map[occ_edge->end.TShape()];
edges.Append(std::move(occ_edge));
}
@ -1145,6 +1146,10 @@ namespace netgen
auto k = faces.Size();
face_map[tshape] = k;
auto occ_face = make_unique<OCCFace>(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));
@ -1194,6 +1199,9 @@ namespace netgen
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(),
@ -1208,9 +1216,8 @@ namespace netgen
add_identifications( edges, edge_map );
add_identifications( faces, face_map );
ProcessIdentifications();
bounding_box = ::netgen::GetBoundingBox( shape );
ProcessIdentifications();
}
@ -1940,32 +1947,6 @@ namespace netgen
return occ2ng( props.CentreOfMass() );
}
void OCCGeometry :: IdentifyEdges(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional<gp_Trsf> opt_trafo)
{
Transformation<3> trafo;
if(opt_trafo)
trafo = occ2ng(*opt_trafo);
else
{
auto cme = GetCenter(me);
auto cyou = GetCenter(you);
trafo = Transformation<3>{cyou-cme};
}
identifications[me.TShape()].push_back( {me.TShape(), you.TShape(), trafo, name, type} );
auto vme = GetVertices(me);
auto vyou = GetVertices(you);
Point<3> pme0 = trafo(occ2ng(vme[0]));
Point<3> pme1 = trafo(occ2ng(vme[1]));
Point<3> pyou = occ2ng(vyou[0]);
bool do_swap = Dist(pme1, pyou) < Dist(pme0, pyou);
for(auto i : Range(2))
identifications[vme[i].TShape()].push_back( {vme[i].TShape(), vyou[do_swap ? 1-i : i].TShape(), trafo, name, type} );
}
bool IsMappedShape(const Transformation<3> & trafo, const TopoDS_Shape & me, const TopoDS_Shape & you)
{
if(me.ShapeType() != you.ShapeType()) return false;
@ -1998,7 +1979,6 @@ namespace netgen
vmap[s] = nullptr;
}
bool all_verts_mapped = true;
for (auto vert : verts_you)
{
auto s = vert.TShape();
@ -2010,16 +1990,18 @@ namespace netgen
return true;
});
if(!vert_mapped)
{
all_verts_mapped = false;
break;
}
return false;
}
return all_verts_mapped;
return true;
}
void OCCGeometry :: IdentifyFaces(const TopoDS_Shape & solid, const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional<gp_Trsf> opt_trafo)
void OCCGeometry :: Identify(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional<gp_Trsf> opt_trafo)
{
auto type_me = me.ShapeType();
auto type_you = you.ShapeType();
if(type_me != type_you)
throw NgException ("Identify: cannot identify different shape types");
Transformation<3> trafo;
if(opt_trafo)
{
@ -2032,25 +2014,33 @@ namespace netgen
trafo = Transformation<3>(cyou-cme);
}
auto faces_me = GetFaces(me);
auto faces_you = GetFaces(you);
ListOfShapes id_me;
ListOfShapes id_you;
for(auto face_me : faces_me)
for(auto face_you : faces_you)
if(auto faces_me = GetFaces(me); faces_me.size()>0)
{
id_me = faces_me;
id_you = GetFaces(you);
}
else if(auto edges_me = GetEdges(me); edges_me.size()>0)
{
id_me = edges_me;
id_you = GetEdges(you);
}
else
{
id_me = GetVertices(me);
id_you = GetVertices(you);
}
for(auto shape_me : id_me)
for(auto shape_you : id_you)
{
if(!IsMappedShape(trafo, face_me, face_you))
if(!IsMappedShape(trafo, shape_me, shape_you))
continue;
identifications[face_me.TShape()].push_back
(OCCIdentification { face_me.TShape(), face_you.TShape(), trafo, name, type });
auto edges_me = GetEdges(me);
auto edges_you = GetEdges(you);
for (auto e_me : edges_me)
for (auto e_you : edges_you)
if(IsMappedShape(trafo, e_me, e_you))
IdentifyEdges(e_me, e_you, name, type, opt_trafo);
identifications[shape_me.TShape()].push_back
(OCCIdentification { shape_me.TShape(), shape_you.TShape(), trafo, name, type });
}
}

4
libsrc/occ/occgeom.hpp
View File

@ -44,6 +44,7 @@ namespace netgen
#define OCCGEOMETRYVISUALIZATIONFULLCHANGE 1 // Compute transformation matrices and redraw
#define OCCGEOMETRYVISUALIZATIONHALFCHANGE 2 // Redraw
bool IsMappedShape(const Transformation<3> & trafo, const TopoDS_Shape & me, const TopoDS_Shape & you);
class EntityVisualizationCode
{
@ -341,8 +342,7 @@ namespace netgen
bool ErrorInSurfaceMeshing ();
// void WriteOCC_STL(char * filename);
static void IdentifyEdges(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional<gp_Trsf> opt_trafo = nullopt);
static void IdentifyFaces(const TopoDS_Shape & solid,const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional<gp_Trsf> opt_trafo = nullopt);
static void Identify(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional<gp_Trsf> opt_trafo = nullopt);
private:
//bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;

75
libsrc/occ/python_occ_shapes.cpp
View File

@ -199,14 +199,14 @@ py::object CastShape(const TopoDS_Shape & s)
template <class TBuilder>
void PropagateIdentifications (TBuilder & builder, TopoDS_Shape shape)
{
std::map<T_Shape, T_Shape> mod_map;
std::map<T_Shape, std::set<T_Shape>> mod_map;
std::map<T_Shape, bool> tshape_handled;
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] = tshape;
mod_map[tshape].insert(tshape);
tshape_handled[tshape] = false;
}
@ -215,11 +215,8 @@ void PropagateIdentifications (TBuilder & builder, TopoDS_Shape shape)
{
auto tshape = e.Current().TShape();
auto & modified = builder.Modified(e.Current());
if(modified.Size()!=1)
continue;
mod_map[tshape] = modified.First().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 })
@ -236,21 +233,33 @@ void PropagateIdentifications (TBuilder & builder, TopoDS_Shape shape)
auto tshape_mapped = mod_map[tshape];
for(auto & ident : OCCGeometry::identifications[tshape])
for(auto ident : OCCGeometry::identifications[tshape])
{
// update existing identification
if(tshape == tshape_mapped)
{
ident.to = mod_map[ident.to];
ident.from = mod_map[ident.from];
}
else
{
OCCIdentification id_new = ident;
id_new.to = mod_map[id_new.to];
id_new.from = mod_map[id_new.from];
OCCGeometry::identifications[mod_map[tshape]].push_back(id_new);
}
// 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);
if(!IsMappedShape(ident.trafo, s_from, s_to))
continue;
OCCIdentification id_new = ident;
id_new.to = to_mapped;
id_new.from = from_mapped;
auto id_owner = from == tshape ? from_mapped : to_mapped;
OCCGeometry::identifications[id_owner].push_back(id_new);
}
}
}
}
@ -1071,27 +1080,9 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
BRepMesh_IncrementalMesh (shape, deflection, true);
})
.def("Identify", [](const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE idtype) {
// only edges supported, by now
auto type = me.ShapeType();
auto tyou = you.ShapeType();
if(type != tyou)
throw NgException ("Identify: cannot identify different shape types");
switch(type)
{
case TopAbs_VERTEX:
case TopAbs_EDGE:
OCCGeometry::IdentifyEdges(me, you, name, idtype);
break;
default:
throw NgException ("Identify: unsupported shape type");
break;
}
}, py::arg("other"), py::arg("name"), py::arg("type")=Identifications::PERIODIC, "Identify shapes for periodic meshing")
.def("Identify", OCCGeometry::IdentifyFaces, "Identify faces",
py::arg("from"), py::arg("to"), py::arg("name"), py::arg("type")=Identifications::PERIODIC, py::arg("trafo")=nullopt)
.def("Identify", OCCGeometry::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)