Merge branch 'fix_occ_internal_faces' into 'master'

Fix meshing of INTERNAL faces with Opencascade

See merge request ngsolve/netgen!591
This commit is contained in:
Hochsteger, Matthias 2023-08-16 18:51:00 +02:00
commit c5b1177151
3 changed files with 24 additions and 4 deletions

View File

@ -61,13 +61,19 @@ namespace netgen
edge_on_face[FORWARD].SetSize(n_edges);
edge_on_face[REVERSED].SetSize(n_edges);
for(auto edge_ : GetEdges(face))
// In case the face is INTERNAL, we need to orient it to FORWARD to get proper orientation for the edges
// (relative to the face) otherwise, all edges are also INTERNAL
auto oriented_face = TopoDS_Face(face);
if(oriented_face.Orientation() == TopAbs_INTERNAL)
oriented_face.Orientation(TopAbs_FORWARD);
for(auto edge_ : GetEdges(oriented_face))
{
auto edge = TopoDS::Edge(edge_);
auto edgenr = geom.GetEdge(edge).nr;
auto & orientation = edge_orientation[edgenr];
double s0, s1;
auto cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1);
auto cof = BRep_Tool::CurveOnSurface (edge, oriented_face, s0, s1);
if(edge.Orientation() == TopAbs_FORWARD || edge.Orientation() == TopAbs_INTERNAL)
{
curve_on_face[FORWARD][edgenr] = cof;
@ -84,7 +90,7 @@ namespace netgen
{
// add reversed edge
auto r_edge = TopoDS::Edge(edge.Reversed());
auto cof = BRep_Tool::CurveOnSurface (r_edge, face, s0, s1);
auto cof = BRep_Tool::CurveOnSurface (r_edge, oriented_face, s0, s1);
curve_on_face[REVERSED][edgenr] = cof;
orientation += REVERSED;
edge_on_face[REVERSED][edgenr] = r_edge;

View File

@ -1208,11 +1208,13 @@ namespace netgen
for(auto f : GetFaces(s))
{
auto & face = GetFace(f);
auto & face = static_cast<OCCFace&>(GetFace(f));
if(face.domin==-1)
face.domin = k;
else
face.domout = k;
if(face.Shape().Orientation() == TopAbs_INTERNAL)
face.domout = k;
}
}

View File

@ -41,3 +41,15 @@ def test_box_and_cyl():
check_volume(cyl, vcyl)
fused = box+cyl
check_volume(fused, 1+vcyl)
def test_internal_face():
occ = pytest.importorskip("netgen.occ")
box = occ.Box((0,0,0), (3, 1, 10))
face = occ.WorkPlane(occ.Axes((1.5,0,0), occ.X, occ.Y)).Rectangle(1, 6).Face()
shape = occ.Glue([box, face])
geo = occ.OCCGeometry(shape)
mesh = geo.GenerateMesh(maxh=0.5)
assert any(mesh.Elements2D().NumPy()['index'] == 8)