From 99a1ddda17be7665c684f9a3b76c3c5cb58a0bfe Mon Sep 17 00:00:00 2001 From: Christophe Bourcier Date: Fri, 4 Nov 2022 17:57:54 +0100 Subject: [PATCH] Fix group of faces identification from dual mesh skin --- src/SMESH_SWIG/smesh_tools.py | 24 ++--- test/SMESH_create_dual_mesh_tpipe.py | 134 +++++++++++++++++++++++++++ test/tests.set | 1 + 3 files changed, 147 insertions(+), 12 deletions(-) create mode 100644 test/SMESH_create_dual_mesh_tpipe.py diff --git a/src/SMESH_SWIG/smesh_tools.py b/src/SMESH_SWIG/smesh_tools.py index e475ab24e..24ccf390f 100644 --- a/src/SMESH_SWIG/smesh_tools.py +++ b/src/SMESH_SWIG/smesh_tools.py @@ -73,18 +73,18 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name grp_poly = mesh2d[id_grp_poly] - # We use the interpolation to remove the element that are not really in - # the group (the ones that are next to a nodes nut not in the group - # will have the sum of their column in the enterpolation matrix equal - # to zero) - rem = mc.MEDCouplingRemapper() - - rem.prepare(grp_poly, grp_tria, "P0P0") - m = rem.getCrudeCSRMatrix() - _, id_to_keep = np.where(m.sum(dtype=np.int64, axis=0) >= 1e-07) - - id_grp_poly = id_grp_poly[id_to_keep.tolist()] - id_grp_poly.setName(grp_name) + # find the extra face cells, on the border of the group (lying on nodes, but outside the group) + id_poly_border = grp_poly.findCellIdsOnBoundary() + + # ids of the cells in grp_poly + id_poly=mc.DataArrayInt64.New(grp_poly.getNumberOfCells(), 1) + id_poly.iota() + + # cells that are really in the group + id_to_keep = id_poly.buildSubstraction(id_poly_border) + + id_grp_poly = id_grp_poly[id_to_keep] + id_grp_poly.setName(grp_name.strip()) myfile.addGroup(-1, id_grp_poly) diff --git a/test/SMESH_create_dual_mesh_tpipe.py b/test/SMESH_create_dual_mesh_tpipe.py new file mode 100644 index 000000000..ebb36e4d2 --- /dev/null +++ b/test/SMESH_create_dual_mesh_tpipe.py @@ -0,0 +1,134 @@ +#!/usr/bin/env python3 + +import sys +import salome +import math + +salome.salome_init() + +import GEOM +from salome.geom import geomBuilder + +geompy = geomBuilder.New() + +# first cylinder +r1 = 0.5 +h1 = 5 + +# second cylinder +r2 = 0.3 +h2 = 3 + +length_piquage = 1.5 + +O = geompy.MakeVertex(0, 0, 0) +OX = geompy.MakeVectorDXDYDZ(1, 0, 0) +OY = geompy.MakeVectorDXDYDZ(0, 1, 0) +OZ = geompy.MakeVectorDXDYDZ(0, 0, 1) +geompy.addToStudy( O, 'O' ) +geompy.addToStudy( OX, 'OX' ) +geompy.addToStudy( OY, 'OY' ) +geompy.addToStudy( OZ, 'OZ' ) + +Cylinder_1 = geompy.MakeCylinderRH(r1, h1) +Cylinder_2 = geompy.MakeCylinderRH(r2, h2) +Rotation_1 = geompy.MakeRotation(Cylinder_2, OY, -90*math.pi/180.0) +Translation_1 = geompy.MakeTranslation(Rotation_1, 0, 0, length_piquage) + +tpipe = geompy.MakeFuseList([Cylinder_1, Translation_1], True, True) +geompy.addToStudy( tpipe, 'tpipe' ) + +Inlet_z = geompy.GetFaceNearPoint(tpipe, O) +geompy.addToStudyInFather( tpipe, Inlet_z, 'Inlet_z' ) + +p_inlet_x = geompy.MakeVertex(-h2, 0, length_piquage) +Inlet_x = geompy.GetFaceNearPoint(tpipe, p_inlet_x) +geompy.addToStudyInFather( tpipe, Inlet_x, 'Inlet_x' ) + +p_outlet = geompy.MakeVertex(0, 0, h1) +Outlet = geompy.GetFaceNearPoint(tpipe, p_outlet) +geompy.addToStudyInFather( tpipe, Outlet, 'Outlet' ) + +Wall = geompy.CreateGroup(tpipe, geompy.ShapeType["FACE"]) +faces = geompy.SubShapeAll(tpipe, geompy.ShapeType["FACE"]) +geompy.UnionList(Wall, faces) +geompy.DifferenceList(Wall, [Inlet_x, Inlet_z, Outlet]) +geompy.addToStudyInFather( tpipe, Wall, 'Wall' ) + +p_corner = geompy.MakeVertex(-r2, 0, length_piquage+r2) +corner = geompy.GetVertexNearPoint(tpipe, p_corner) +geompy.addToStudyInFather( tpipe, corner, 'corner' ) + +geom_groups = [Inlet_x, Inlet_z, Outlet, Wall] + +volumeGEOM = geompy.BasicProperties(tpipe)[2] + +### +### SMESH component +### + +import SMESH, SALOMEDS +from salome.smesh import smeshBuilder + +from salome.StdMeshers import StdMeshersBuilder + +smesh = smeshBuilder.New() + +# Coarse mesh with default hypothesis +Mesh_1 = smesh.Mesh(tpipe, "Mesh_coarse") + +Mesh_1.Triangle(algo=smeshBuilder.NETGEN_1D2D) + +Mesh_1.Tetrahedron(algo=smeshBuilder.NETGEN_3D) +isDone = Mesh_1.Compute() + +# Create groups +d_mesh_groups = {} +for geom_group in geom_groups: + gr = Mesh_1.Group(geom_group) + gr_name = gr.GetName() + d_mesh_groups[gr_name] = gr + +# Check tetra mesh volume +mesh_volume = smesh.GetVolume(Mesh_1) +shape_volume = geompy.BasicProperties(tpipe)[2] +assert abs(mesh_volume-shape_volume)/shape_volume < 0.05 + +dual_Mesh_raw_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_raw_1', False) + +# check polyhedrons +assert dual_Mesh_raw_1.NbPolyhedrons() > 0 +assert dual_Mesh_raw_1.NbPolyhedrons() == dual_Mesh_raw_1.NbVolumes() + +# Check dual mesh volume +dual_raw_volume = dual_Mesh_raw_1.GetVolume() +assert abs(dual_raw_volume-shape_volume)/shape_volume < 0.11 + +# Check groups +dual_Mesh_raw_groups = dual_Mesh_raw_1.GetGroups() +dual_Mesh_raw_group_names = dual_Mesh_raw_1.GetGroupNames() + +assert len(dual_Mesh_raw_groups) == 4 + +for gr_dual, gr_name in zip(dual_Mesh_raw_groups, dual_Mesh_raw_group_names): + gr_name = gr_name.strip() + gr_tri = d_mesh_groups[gr_name] + area_gr_tri = smesh.GetArea(gr_tri) + area_gr_dual = smesh.GetArea(gr_dual) + print(gr_name) + print("Area tri: ", area_gr_tri) + print("Area dual:", area_gr_dual) + assert abs(area_gr_tri-area_gr_dual)/area_gr_tri < 1e-3 + +#dual_Mesh_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_1', True) + + +#dual_volume = dual_Mesh_1.GetVolume() +#dual_raw_volume = dual_Mesh_raw_1.GetVolume() +#print("dual_volume: ", dual_volume) +#print("dual_raw_volume: ", dual_raw_volume) + +#assert (dual_volume >= dual_raw_volume) + +if salome.sg.hasDesktop(): + salome.sg.updateObjBrowser() diff --git a/test/tests.set b/test/tests.set index fddb00249..21491d807 100644 --- a/test/tests.set +++ b/test/tests.set @@ -63,6 +63,7 @@ SET(BAD_TESTS SMESH_test2.py SMESH_test4.py SMESH_create_dual_mesh_adapt.py + SMESH_create_dual_mesh_tpipe.py ) IF(NOT WIN32) LIST(APPEND BAD_TESTS -- 2.39.2