Code cleaning in dual mesh

This commit is contained in:
Christophe Bourcier 2022-11-16 10:44:02 +01:00
parent 21ac1ad91f
commit 11d9d2377c
3 changed files with 51 additions and 34 deletions

View File

@ -789,7 +789,6 @@ class smeshBuilder( SMESH._objref_SMESH_Gen, object ):
""" """
if isinstance( mesh, Mesh ): if isinstance( mesh, Mesh ):
mesh = mesh.GetMesh() mesh = mesh.GetMesh()
print("calling createdualmesh from Python")
dualMesh = SMESH._objref_SMESH_Gen.CreateDualMesh(self, mesh, meshName, adaptToShape) dualMesh = SMESH._objref_SMESH_Gen.CreateDualMesh(self, mesh, meshName, adaptToShape)
return Mesh(self, self.geompyD, dualMesh) return Mesh(self, self.geompyD, dualMesh)

View File

@ -1,12 +1,7 @@
#!/usr/bin/env python3 #!/usr/bin/env python3
import sys
import salome import salome
import medcoupling as mc import medcoupling as mc
from math import pi
import numpy as np
#salome.salome_init()
import GEOM import GEOM
from salome.geom import geomBuilder from salome.geom import geomBuilder
@ -20,30 +15,34 @@ smesh = smeshBuilder.New()
from salome.kernel.logger import Logger from salome.kernel.logger import Logger
logger = Logger("salome.smesh.smesh_tools") logger = Logger("salome.smesh.smesh_tools")
logger.setLevel("DEBUG") logger.setLevel("WARNING")
# prefix for groups with internal usage # prefix for groups with internal usage
# i.e. used to transfer the faces and edges sub-shapes ids to the mesh # i.e. used to transfer the faces and edges sub-shapes ids to the mesh
__prefix = "____" __prefix = "____"
def __getIdsGrpDualFromOrig(mc_mesh_file, grp_name, mesh2d, grp_level): def __getIdsGrpDualFromOrig(mc_mesh_file, grp_name, mesh2d, grp_level):
""" Identify the polygonal cells ids matching the original group on the original mesh (before dual mesh) """ Identify the polygonal cells ids matching the original group on the
original mesh (before dual mesh)
Args: Args:
mc_mesh_file (MEDFileMesh): mesh on which to read the group mc_mesh_file (MEDFileMesh): mesh on which to read the group
grp_name (string): name of the group to read grp_name (string): name of the group to read
mesh2d (MEDCouplingUMesh): mesh at lower level (-1 or -2) containing faces or segments cells mesh2d (MEDCouplingUMesh): mesh at lower level (-1 or -2) containing
faces or segments cells
grp_level (int): level on which to load the group (-1 or -2) grp_level (int): level on which to load the group (-1 or -2)
Returns: Returns:
id_grp_poly (DataArrayInt64): ids of cells mathcing the group. None if the group has not been found. id_grp_poly (DataArrayInt64): ids of cells mathcing the group. None if
the group has not been found.
nodes_added (DataArrayInt64): new nodes added on the dual mesh nodes_added (DataArrayInt64): new nodes added on the dual mesh
""" """
try: try:
grp_tria = mc_mesh_file.getGroup(grp_level, grp_name) grp_tria = mc_mesh_file.getGroup(grp_level, grp_name)
except: except:
logger.debug("""No group found for %s at level %i. logger.debug("""No group found for %s at level %i.
It is normal behaviour for degenerated geom edge."""%(grp_name, grp_level)) It is normal behaviour for degenerated geom edge."""\
%(grp_name, grp_level))
return None, None return None, None
# Retrieve the nodes in group # Retrieve the nodes in group
orig_grp_nodes = grp_tria.computeFetchedNodeIds() orig_grp_nodes = grp_tria.computeFetchedNodeIds()
@ -52,14 +51,16 @@ def __getIdsGrpDualFromOrig(mc_mesh_file, grp_name, mesh2d, grp_level):
grp_poly = mesh2d[id_grp_poly] grp_poly = mesh2d[id_grp_poly]
if grp_poly.getNumberOfCells() == 0: if grp_poly.getNumberOfCells() == 0:
logger.debug("""No poly cell found for %s at level %i."""%(grp_name, grp_level)) logger.debug("""No poly cell found for %s at level %i."""\
%(grp_name, grp_level))
return None, None return None, None
# find the extra face cells, on the border of the group (lying on nodes, but outside the group) # find the extra face cells, on the border of the group (lying on nodes,
# but outside the group)
id_poly_border = grp_poly.findCellIdsOnBoundary() id_poly_border = grp_poly.findCellIdsOnBoundary()
# ids of the cells in grp_poly # ids of the cells in grp_poly
id_poly=mc.DataArrayInt64.New(grp_poly.getNumberOfCells(), 1) id_poly = mc.DataArrayInt64.New(grp_poly.getNumberOfCells(), 1)
id_poly.iota() id_poly.iota()
# cells that are really in the group # cells that are really in the group
@ -77,12 +78,14 @@ def __getIdsGrpDualFromOrig(mc_mesh_file, grp_name, mesh2d, grp_level):
return id_grp_poly, nodes_added return id_grp_poly, nodes_added
def __projectNodeOnSubshape(nodes, subshape, coords): def __projectNodeOnSubshape(nodes, subshape, coords):
""" Project nodes on a sub-shape (face or edge) and update the mesh coordinates """ Project nodes on a sub-shape (face or edge) and update the mesh
coordinates
Args: Args:
nodes (DataArrayInt): nodes ids to project nodes (DataArrayInt): nodes ids to project
subshape (GEOM object): face or edge on which to project the nodes subshape (GEOM object): face or edge on which to project the nodes
coords (DataArrayDouble): coordinates of the mesh to update. These coordinates are modified inside this function. coords (DataArrayDouble): coordinates of the mesh to update. These
coordinates are modified inside this function.
""" """
for i in nodes: for i in nodes:
x, y, z = coords[i].getValues() x, y, z = coords[i].getValues()
@ -110,12 +113,15 @@ def __deleteObj(theObj):
aStudyBuilder.RemoveObjectWithChildren(SO) aStudyBuilder.RemoveObjectWithChildren(SO)
pass pass
def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name="MESH"): def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True,
mesh_name="MESH"):
""" Create a dual of the mesh in input_file into output_file """ Create a dual of the mesh in input_file into output_file
Args: Args:
mesh_ior (string): corba Id of the Tetrahedron mesh mesh_ior (string): corba Id of the Tetrahedron mesh
output_file (string): dual mesh file output_file (string): dual mesh file
adapt_to_shape (bool): If True will project boundary points on shape
mesh_name (string): Name of the dual Mesh
""" """
# Import mesh from file # Import mesh from file
mesh = salome.orb.string_to_object(mesh_ior) mesh = salome.orb.string_to_object(mesh_ior)
@ -129,12 +135,15 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name
faces_ids = geompy.GetSubShapesIDs(shape, faces) faces_ids = geompy.GetSubShapesIDs(shape, faces)
# Create group with each face # Create group with each face
# so that we don't need GetFaceNearPoint to get the face to project the point to # so that we don't need GetFaceNearPoint to get the face to project the
# point to
id2face = {} id2face = {}
id2face_edges_ids = {} id2face_edges_ids = {}
mesh_groups = [] mesh_groups = []
for face, face_id in zip(faces, faces_ids): for face, face_id in zip(faces, faces_ids):
gr_mesh = mesh.CreateGroupFromGEOM(SMESH.FACE, '%sface_%i'%(__prefix, face_id), face) gr_mesh = mesh.CreateGroupFromGEOM(SMESH.FACE,
'%sface_%i'%(__prefix, face_id),
face)
id2face[face_id] = face id2face[face_id] = face
mesh_groups.append(gr_mesh) mesh_groups.append(gr_mesh)
# get the edges bounding this face # get the edges bounding this face
@ -148,7 +157,9 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name
id2edge = {} id2edge = {}
for edge, edge_id in zip(edges, edges_ids): for edge, edge_id in zip(edges, edges_ids):
gr_mesh = mesh.CreateGroupFromGEOM(SMESH.EDGE, '%sedge_%i'%(__prefix, edge_id), edge) gr_mesh = mesh.CreateGroupFromGEOM(SMESH.EDGE,
'%sedge_%i'%(__prefix, edge_id),
edge)
id2edge[edge_id] = edge id2edge[edge_id] = edge
mesh_groups.append(gr_mesh) mesh_groups.append(gr_mesh)
@ -170,8 +181,6 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name
polyh_coords = polyh.getCoords() polyh_coords = polyh.getCoords()
edges_group_names = mc_mesh_file.getGroupsOnSpecifiedLev(-2)
treated_edges = [] treated_edges = []
mc_groups = [] mc_groups = []
@ -183,7 +192,8 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name
logger.debug("Transferring group: "+ grp_name) logger.debug("Transferring group: "+ grp_name)
# get the polygons ids on the dual mesh from the triangles group # get the polygons ids on the dual mesh from the triangles group
id_grp_poly, nodes_added_on_tri = __getIdsGrpDualFromOrig(mc_mesh_file, grp_name, mesh2d, -1) id_grp_poly, nodes_added_on_tri = \
__getIdsGrpDualFromOrig(mc_mesh_file, grp_name, mesh2d, -1)
if id_grp_poly is not None and grp_name[:4] == __prefix: if id_grp_poly is not None and grp_name[:4] == __prefix:
# This group is on a specific geom face # This group is on a specific geom face
@ -200,19 +210,23 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name
grp_seg_name = "%sedge_%i"%(__prefix, edge_id) grp_seg_name = "%sedge_%i"%(__prefix, edge_id)
# get the segs on the dual mesh from the segs group # get the segs on the dual mesh from the segs group
id_grp_seg, nodes_added_on_segs = __getIdsGrpDualFromOrig(mc_mesh_file, grp_seg_name, mesh1d, -2) id_grp_seg, nodes_added_on_segs = __getIdsGrpDualFromOrig(\
mc_mesh_file, grp_seg_name, mesh1d, -2)
# project new nodes on its geom edge # project new nodes on its geom edge
# (if the group exists on this edge and it has not already been treated) # (if the group exists on this edge and it has not already been
# treated)
if id_grp_seg is not None: if id_grp_seg is not None:
if edge_id not in treated_edges: if edge_id not in treated_edges:
edge = id2edge[edge_id] edge = id2edge[edge_id]
__projectNodeOnSubshape(nodes_added_on_segs, edge, polyh_coords) __projectNodeOnSubshape(nodes_added_on_segs,
edge, polyh_coords)
else: else:
treated_edges.append(edge_id) treated_edges.append(edge_id)
# remove these nodes from the nodes to project on face # remove these nodes from the nodes to project on face
nodes_added_on_tri = nodes_added_on_tri.buildSubstraction(nodes_added_on_segs) nodes_added_on_tri = \
nodes_added_on_tri.buildSubstraction(nodes_added_on_segs)
# project new nodes on its geom face # project new nodes on its geom face
__projectNodeOnSubshape(nodes_added_on_tri, face, polyh_coords) __projectNodeOnSubshape(nodes_added_on_tri, face, polyh_coords)
@ -236,5 +250,5 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name
if adapt_to_shape: if adapt_to_shape:
# delete temporary groups # delete temporary groups
for gr in mesh_groups: for grp in mesh_groups:
__deleteObj(gr) __deleteObj(grp)

View File

@ -79,7 +79,9 @@ Mesh_1 = smesh.Mesh(tpipe, "Mesh_coarse")
Mesh_1.Triangle(algo=smeshBuilder.NETGEN_1D2D) Mesh_1.Triangle(algo=smeshBuilder.NETGEN_1D2D)
Mesh_1.Tetrahedron(algo=smeshBuilder.NETGEN_3D) algo3d = Mesh_1.Tetrahedron(algo=smeshBuilder.NETGEN_3D)
algo3d.MaxElementVolume(0.0002)
isDone = Mesh_1.Compute() isDone = Mesh_1.Compute()
# Create groups # Create groups
@ -92,9 +94,9 @@ for geom_group in geom_groups:
d_geom_groups[gr_name] = geom_group d_geom_groups[gr_name] = geom_group
# Check tetra mesh volume # Check tetra mesh volume
mesh_volume = smesh.GetVolume(Mesh_1) tetra_volume = smesh.GetVolume(Mesh_1)
shape_volume = geompy.BasicProperties(tpipe)[2] shape_volume = geompy.BasicProperties(tpipe)[2]
assert abs(mesh_volume-shape_volume)/shape_volume < 0.05 assert abs(tetra_volume-shape_volume)/shape_volume < 0.05
dual_Mesh_raw_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_raw_1', False) dual_Mesh_raw_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_raw_1', False)
@ -104,7 +106,7 @@ assert dual_Mesh_raw_1.NbPolyhedrons() == dual_Mesh_raw_1.NbVolumes()
# Check dual mesh volume # Check dual mesh volume
dual_raw_volume = dual_Mesh_raw_1.GetVolume() dual_raw_volume = dual_Mesh_raw_1.GetVolume()
assert abs(dual_raw_volume-shape_volume)/shape_volume < 0.11 assert abs(dual_raw_volume-shape_volume)/shape_volume < 0.12
# Check groups # Check groups
dual_Mesh_raw_groups = dual_Mesh_raw_1.GetGroups() dual_Mesh_raw_groups = dual_Mesh_raw_1.GetGroups()
@ -117,11 +119,13 @@ dual_Mesh_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_1', True)
# Check dual mesh volume # Check dual mesh volume
dual_volume = dual_Mesh_1.GetVolume() dual_volume = dual_Mesh_1.GetVolume()
print("shape_volume: ", shape_volume)
print("tetra_volume: ", tetra_volume)
print("dual_volume: ", dual_volume) print("dual_volume: ", dual_volume)
print("dual_raw_volume: ", dual_raw_volume) print("dual_raw_volume: ", dual_raw_volume)
assert (dual_volume >= dual_raw_volume) assert (dual_volume >= dual_raw_volume)
assert abs(dual_volume-shape_volume)/shape_volume < 0.11 assert abs(dual_volume-shape_volume)/shape_volume < 0.12
# Check groups # Check groups
dual_Mesh_groups = dual_Mesh_1.GetGroups() dual_Mesh_groups = dual_Mesh_1.GetGroups()