src/SMESH_SWIG/smesh_tools.py

   1 #!/usr/bin/env python3
   2
   3 import salome
   4 import medcoupling as mc
   5
   6 import GEOM
   7 from salome.geom import geomBuilder
   8
   9 geompy = geomBuilder.New()
  10
  11 import  SMESH, SALOMEDS
  12 from salome.smesh import smeshBuilder
  13
  14 smesh = smeshBuilder.New()
  15
  16 from salome.kernel.logger import Logger
  17 logger = Logger("salome.smesh.smesh_tools")
  18 logger.setLevel("WARNING")
  19
  20 # prefix for groups with internal usage
  21 # i.e. used to transfer the faces and edges sub-shapes ids to the mesh
  22 __prefix = "____"
  23
  24 def __getIdsGrpDualFromOrig(mc_mesh_file, grp_name, mesh2d, grp_level):
  25     """ Identify the polygonal cells ids matching the original group on the
  26         original mesh (before dual mesh)
  27
  28     Args:
  29         mc_mesh_file (MEDFileMesh): mesh on which to read the group
  30         grp_name (string): name of the group to read
  31         mesh2d (MEDCouplingUMesh): mesh at lower level (-1 or -2) containing
  32                                    faces or segments cells
  33         grp_level (int): level on which to load the group (-1 or -2)
  34
  35     Returns:
  36         id_grp_poly (DataArrayInt64): ids of cells mathcing the group. None if
  37                                       the group has not been found.
  38         nodes_added (DataArrayInt64): new nodes added on the dual mesh
  39     """
  40     try:
  41         grp_tria = mc_mesh_file.getGroup(grp_level, grp_name)
  42     except:
  43         logger.debug("""No group found for %s at level %i.
  44                      It is normal behaviour for degenerated geom edge."""\
  45                      %(grp_name, grp_level))
  46         return None, None
  47     # Retrieve the nodes in group
  48     orig_grp_nodes = grp_tria.computeFetchedNodeIds()
  49     # Find all the cells lying on one of the nodes
  50     id_grp_poly = mesh2d.getCellIdsLyingOnNodes(orig_grp_nodes, False)
  51
  52     grp_poly = mesh2d[id_grp_poly]
  53     if grp_poly.getNumberOfCells() == 0:
  54         logger.debug("""No poly cell found for %s at level %i."""\
  55                      %(grp_name, grp_level))
  56         return None, None
  57
  58     # find the extra face cells, on the border of the group (lying on nodes,
  59     # but outside the group)
  60     id_poly_border = grp_poly.findCellIdsOnBoundary()
  61
  62     # ids of the cells in grp_poly
  63     id_poly = mc.DataArrayInt64.New(grp_poly.getNumberOfCells(), 1)
  64     id_poly.iota()
  65
  66     # cells that are really in the group
  67     id_to_keep = id_poly.buildSubstraction(id_poly_border)
  68
  69     id_grp_poly = id_grp_poly[id_to_keep]
  70
  71     id_grp_poly.setName(grp_name.strip())
  72
  73     # get nodes added on this group
  74     grp_poly = mesh2d[id_grp_poly]
  75     grp_nodes_poly = grp_poly.computeFetchedNodeIds()
  76     nodes_added = grp_nodes_poly.buildSubstraction(orig_grp_nodes)
  77
  78     return id_grp_poly, nodes_added
  79
  80 def __projectNodeOnSubshape(nodes, subshape, coords):
  81     """ Project nodes on a sub-shape (face or edge) and update the mesh
  82         coordinates
  83
  84     Args:
  85         nodes (DataArrayInt): nodes ids to project
  86         subshape (GEOM object): face or edge on which to project the nodes
  87         coords (DataArrayDouble): coordinates of the mesh to update. These
  88                                   coordinates are modified inside this function.
  89     """
  90     for i in nodes:
  91         x, y, z = coords[i].getValues()
  92         vertex = geompy.MakeVertex(x, y, z)
  93         try:
  94             prj = geompy.MakeProjection(vertex, subshape)
  95         except:
  96             logger.warning("Projection failed for %.5f %.5f %.5f but we continue with next node"%(x, y, z))
  97             continue
  98         new_coor = geompy.PointCoordinates(prj)
  99         # update its coordinates in the mesh
 100         coords[i] = new_coor
 101     pass
 102
 103 def __deleteObj(theObj):
 104     """ Delete object from a Study
 105
 106     Args:
 107         theObj (GEOM or SMESH object): object to remove from the study
 108     """
 109     aStudy = salome.myStudy
 110     aStudyBuilder = aStudy.NewBuilder()
 111     SO = aStudy.FindObjectIOR(aStudy.ConvertObjectToIOR(theObj))
 112     if SO is not None:
 113         aStudyBuilder.RemoveObjectWithChildren(SO)
 114     pass
 115
 116 def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True,
 117                            mesh_name="MESH"):
 118     """ Create a dual of the mesh in input_file into output_file
 119
 120     Args:
 121         mesh_ior (string): corba Id of the Tetrahedron mesh
 122         output_file (string): dual mesh file
 123         adapt_to_shape (bool): If True will project boundary points on shape
 124         mesh_name (string): Name of the dual Mesh
 125     """
 126     # Import mesh from file
 127     mesh = salome.orb.string_to_object(mesh_ior)
 128     if not mesh:
 129         raise Exception("Could not find mesh using id: ", mesh_ior)
 130
 131     shape = mesh.GetShapeToMesh()
 132
 133     if adapt_to_shape:
 134         faces = geompy.SubShapeAll(shape, geompy.ShapeType["FACE"])
 135         faces_ids = geompy.GetSubShapesIDs(shape, faces)
 136
 137         # Create group with each face
 138         # so that we don't need GetFaceNearPoint to get the face to project the
 139         # point to
 140         id2face = {}
 141         id2face_edges_ids = {}
 142         mesh_groups = []
 143         for face, face_id in zip(faces, faces_ids):
 144             gr_mesh = mesh.CreateGroupFromGEOM(SMESH.FACE,
 145                                                '%sface_%i'%(__prefix, face_id),
 146                                                face)
 147             id2face[face_id] = face
 148             mesh_groups.append(gr_mesh)
 149             # get the edges bounding this face
 150             # so that we can project the nodes on edges before nodes on faces
 151             face_edges = geompy.SubShapeAll(face, geompy.ShapeType["EDGE"])
 152             face_edges_ids = geompy.GetSubShapesIDs(shape, face_edges)
 153             id2face_edges_ids[face_id] = face_edges_ids
 154
 155         edges = geompy.SubShapeAll(shape, geompy.ShapeType["EDGE"])
 156         edges_ids = geompy.GetSubShapesIDs(shape, edges)
 157
 158         id2edge = {}
 159         for edge, edge_id in zip(edges, edges_ids):
 160             gr_mesh = mesh.CreateGroupFromGEOM(SMESH.EDGE,
 161                                                '%sedge_%i'%(__prefix, edge_id),
 162                                                edge)
 163             id2edge[edge_id] = edge
 164             mesh_groups.append(gr_mesh)
 165
 166     # We got a meshProxy so we need to convert pointer to MEDCoupling
 167     int_ptr = mesh.ExportMEDCoupling(True, True)
 168     dab = mc.FromPyIntPtrToDataArrayByte(int_ptr)
 169     mc_mesh_file = mc.MEDFileMesh.New(dab)
 170     tetras = mc_mesh_file[0]
 171     # End of SMESH -> MEDCoupling part for dualmesh
 172
 173     tetras = mc.MEDCoupling1SGTUMesh(tetras)
 174
 175     # Create the polyhedra from the tetrahedra (main goal of this function)
 176     polyh = tetras.computeDualMesh()
 177
 178     ## Adding skin + transfering groups on faces from tetras mesh
 179     mesh2d = polyh.buildUnstructured().computeSkin()
 180     mesh2d.setName(mesh_name)
 181
 182     polyh_coords = polyh.getCoords()
 183
 184     treated_edges = []
 185
 186     mc_groups = []
 187     for grp_name in mc_mesh_file.getGroupsOnSpecifiedLev(-1):
 188         # This group is created by the export
 189         if grp_name == "Group_Of_All_Faces":
 190             logger.debug("Skipping group: "+ grp_name)
 191             continue
 192         logger.debug("Transferring group: "+ grp_name)
 193
 194         # get the polygons ids on the dual mesh from the triangles group
 195         id_grp_poly, nodes_added_on_tri = \
 196             __getIdsGrpDualFromOrig(mc_mesh_file, grp_name, mesh2d, -1)
 197
 198         if id_grp_poly is not None and grp_name[:4] == __prefix:
 199             # This group is on a specific geom face
 200             face_id = grp_name.split("_")[-1]
 201             face_id = int(face_id)
 202             face = id2face[face_id]
 203
 204             # for each face, get the edges bounding it
 205             grp_poly = mesh2d[id_grp_poly]
 206             mesh1d = grp_poly.computeSkin()
 207
 208             face_edges_id = id2face_edges_ids[face_id]
 209             for edge_id in face_edges_id:
 210                 grp_seg_name = "%sedge_%i"%(__prefix, edge_id)
 211
 212                 # get the segs on the dual mesh from the segs group
 213                 id_grp_seg, nodes_added_on_segs = __getIdsGrpDualFromOrig(\
 214                                 mc_mesh_file, grp_seg_name, mesh1d, -2)
 215
 216                 # project new nodes on its geom edge
 217                 # (if the group exists on this edge and it has not already been
 218                 # treated)
 219                 if id_grp_seg is not None:
 220                     if edge_id not in treated_edges:
 221                         edge = id2edge[edge_id]
 222                         __projectNodeOnSubshape(nodes_added_on_segs,
 223                                                 edge, polyh_coords)
 224                     else:
 225                         treated_edges.append(edge_id)
 226
 227                     # remove these nodes from the nodes to project on face
 228                     nodes_added_on_tri = \
 229                        nodes_added_on_tri.buildSubstraction(nodes_added_on_segs)
 230
 231             # project new nodes on its geom face
 232             __projectNodeOnSubshape(nodes_added_on_tri, face, polyh_coords)
 233         else:
 234             # add the group to write it
 235             mc_groups.append(id_grp_poly)
 236
 237     # Creating output file
 238     logger.debug("Creating file with mesh: "+mesh_name)
 239     myfile = mc.MEDFileUMesh()
 240     myfile.setName(mesh_name)
 241     polyh.setName(mesh_name)
 242     myfile.setMeshAtLevel(0, polyh)
 243     myfile.setMeshAtLevel(-1, mesh2d)
 244
 245     for group in mc_groups:
 246         myfile.addGroup(-1, group)
 247
 248     logger.debug("Writing dual mesh in: "+output_file)
 249     myfile.write(output_file, 2)
 250
 251     if adapt_to_shape:
 252         # delete temporary groups
 253         for grp in mesh_groups:
 254             __deleteObj(grp)