src/SMESH_SWIG/smesh_tools.py

   1 #!/usr/bin/env python3
   2
   3 import sys
   4 import salome
   5 import medcoupling as mc
   6 from math import pi
   7 import numpy as np
   8
   9 #salome.salome_init()
  10
  11 import GEOM
  12 from salome.geom import geomBuilder
  13
  14 geompy = geomBuilder.New()
  15
  16 import  SMESH, SALOMEDS
  17 from salome.smesh import smeshBuilder
  18
  19 smesh = smeshBuilder.New()
  20
  21 from salome.kernel.logger import Logger
  22 logger = Logger("salome.smesh.smesh_tools")
  23 logger.setLevel("DEBUG")
  24
  25 def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name="MESH"):
  26     """ Create a dual of the mesh in input_file into output_file
  27
  28     Args:
  29         mesh_ior (string): corba Id of the Tetrahedron mesh
  30         output_file (string): dual mesh file
  31     """
  32     # Import mesh from file
  33     mesh = salome.orb.string_to_object(mesh_ior)
  34     if not mesh:
  35         raise Exception("Could not find mesh using id: ", mesh_ior)
  36
  37     shape = mesh.GetShapeToMesh()
  38
  39     # Creating output file
  40     logger.debug("Creating file with mesh: "+mesh_name)
  41     myfile = mc.MEDFileUMesh()
  42     myfile.setName(mesh_name)
  43
  44
  45     # We got a meshProxy so we need to convert pointer to MEDCoupling
  46     int_ptr = mesh.ExportMEDCoupling(True, True)
  47     dab = mc.FromPyIntPtrToDataArrayByte(int_ptr)
  48     mc_mesh_file = mc.MEDFileMesh.New(dab)
  49     tetras = mc_mesh_file[0]
  50     # End of SMESH -> MEDCoupling part for dualmesh
  51
  52     tetras = mc.MEDCoupling1SGTUMesh(tetras)
  53     polyh = tetras.computeDualMesh()
  54
  55     ## Adding skin + transfering groups on faces from tetras mesh
  56     mesh2d = polyh.buildUnstructured().computeSkin()
  57     mesh2d.setName(mesh_name)
  58     myfile.setMeshAtLevel(-1, mesh2d)
  59
  60
  61     for grp_name in mc_mesh_file.getGroupsOnSpecifiedLev(-1):
  62         # This group is created by the export
  63         if grp_name == "Group_Of_All_Faces":
  64             logger.debug("Skipping group: "+ grp_name)
  65             continue
  66         logger.debug("Transferring group: "+ grp_name)
  67
  68         grp_tria = mc_mesh_file.getGroup(-1, grp_name)
  69         # Retrieve the nodes in group
  70         grp_nodes = grp_tria.computeFetchedNodeIds()
  71         # Find all the cells lying on one of the nodes
  72         id_grp_poly = mesh2d.getCellIdsLyingOnNodes(grp_nodes, False)
  73
  74         grp_poly = mesh2d[id_grp_poly]
  75
  76         # We use the interpolation to remove the element that are not really in
  77         # the group (the ones that are next to a nodes nut not in the group
  78         # will have the sum of their column in the enterpolation matrix equal
  79         # to zero)
  80         rem = mc.MEDCouplingRemapper()
  81
  82         rem.prepare(grp_poly, grp_tria, "P0P0")
  83         m = rem.getCrudeCSRMatrix()
  84         _, id_to_keep = np.where(m.sum(dtype=np.int64, axis=0) >= 1e-07)
  85
  86         id_grp_poly = id_grp_poly[id_to_keep.tolist()]
  87         id_grp_poly.setName(grp_name)
  88
  89         myfile.addGroup(-1, id_grp_poly)
  90
  91     # Getting list of new points added on the skin
  92     skin = tetras.buildUnstructured().computeSkin()
  93     skin_polyh = polyh.buildUnstructured().computeSkin()
  94     allNodesOnSkinPolyh = skin_polyh.computeFetchedNodeIds()
  95     allNodesOnSkin = skin.computeFetchedNodeIds()
  96     ptsAdded = allNodesOnSkinPolyh.buildSubstraction(allNodesOnSkin)
  97     ptsAddedMesh = mc.MEDCouplingUMesh.Build0DMeshFromCoords( skin_polyh.getCoords()[ptsAdded] )
  98
  99     if adapt_to_shape:
 100         logger.debug("Adapting to shape")
 101         ptsAddedCoo = ptsAddedMesh.getCoords()
 102         ptsAddedCooModified = ptsAddedCoo[:]
 103
 104         # Matching faces with their ids
 105         faces = geompy.ExtractShapes(shape, geompy.ShapeType["FACE"], True)
 106         id2face = {}
 107         for face in faces:
 108             id2face[face.GetSubShapeIndices()[0]] = face
 109
 110         ## Projecting each points added by the dual mesh on the surface it is
 111         # associated with
 112         for i, tup in enumerate(ptsAddedCooModified):
 113             vertex = geompy.MakeVertex(*tuple(tup))
 114             shapes = geompy.GetShapesNearPoint(shape, vertex,
 115                                                geompy.ShapeType["FACE"])
 116             prj = geompy.MakeProjection(vertex,
 117                                         id2face[shapes.GetSubShapeIndices()[0]])
 118             new_coor = geompy.PointCoordinates(prj)
 119             ptsAddedCooModified[i] = new_coor
 120
 121         polyh.getCoords()[ptsAdded] = ptsAddedCooModified
 122
 123     polyh.setName(mesh_name)
 124     myfile.setMeshAtLevel(0, polyh)
 125
 126     logger.debug("Writting dual mesh in :"+output_file)
 127     myfile.write(output_file, 2)