5 import medcoupling as mc
12 from salome.geom import geomBuilder
14 geompy = geomBuilder.New()
16 import SMESH, SALOMEDS
17 from salome.smesh import smeshBuilder
19 smesh = smeshBuilder.New()
21 from salome.kernel.logger import Logger
22 logger = Logger("salome.smesh.smesh_tools")
23 logger.setLevel("DEBUG")
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
29 mesh_ior (string): corba Id of the Tetrahedron mesh
30 output_file (string): dual mesh file
32 # Import mesh from file
33 mesh = salome.orb.string_to_object(mesh_ior)
35 raise Exception("Could not find mesh using id: ", mesh_ior)
37 shape = mesh.GetShapeToMesh()
39 # Creating output file
40 logger.debug("Creating file with mesh: "+mesh_name)
41 myfile = mc.MEDFileUMesh()
42 myfile.setName(mesh_name)
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
52 tetras = mc.MEDCoupling1SGTUMesh(tetras)
53 polyh = tetras.computeDualMesh()
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)
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)
66 logger.debug("Transferring group: "+ grp_name)
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)
74 grp_poly = mesh2d[id_grp_poly]
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
80 rem = mc.MEDCouplingRemapper()
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)
86 id_grp_poly = id_grp_poly[id_to_keep.tolist()]
87 id_grp_poly.setName(grp_name)
89 myfile.addGroup(-1, id_grp_poly)
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] )
100 logger.debug("Adapting to shape")
101 ptsAddedCoo = ptsAddedMesh.getCoords()
102 ptsAddedCooModified = ptsAddedCoo[:]
104 # Matching faces with their ids
105 faces = geompy.ExtractShapes(shape, geompy.ShapeType["FACE"], True)
108 id2face[face.GetSubShapeIndices()[0]] = face
110 ## Projecting each points added by the dual mesh on the surface it is
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
121 polyh.getCoords()[ptsAdded] = ptsAddedCooModified
123 polyh.setName(mesh_name)
124 myfile.setMeshAtLevel(0, polyh)
126 logger.debug("Writting dual mesh in :"+output_file)
127 myfile.write(output_file, 2)