4 import medcoupling as mc
7 from salome.geom import geomBuilder
9 geompy = geomBuilder.New()
11 import SMESH, SALOMEDS
12 from salome.smesh import smeshBuilder
14 smesh = smeshBuilder.New()
16 from salome.kernel.logger import Logger
17 logger = Logger("salome.smesh.smesh_tools")
18 logger.setLevel("WARNING")
20 # prefix for groups with internal usage
21 # i.e. used to transfer the faces and edges sub-shapes ids to the mesh
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)
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)
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
41 grp_tria = mc_mesh_file.getGroup(grp_level, grp_name)
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))
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)
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))
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()
62 # ids of the cells in grp_poly
63 id_poly = mc.DataArrayInt64.New(grp_poly.getNumberOfCells(), 1)
66 # cells that are really in the group
67 id_to_keep = id_poly.buildSubstraction(id_poly_border)
69 id_grp_poly = id_grp_poly[id_to_keep]
71 id_grp_poly.setName(grp_name.strip())
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)
78 return id_grp_poly, nodes_added
80 def __projectNodeOnSubshape(nodes, subshape, coords):
81 """ Project nodes on a sub-shape (face or edge) and update the mesh
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.
91 x, y, z = coords[i].getValues()
92 vertex = geompy.MakeVertex(x, y, z)
94 prj = geompy.MakeProjection(vertex, subshape)
96 logger.warning("Projection failed for %.5f %.5f %.5f but we continue with next node"%(x, y, z))
98 new_coor = geompy.PointCoordinates(prj)
99 # update its coordinates in the mesh
103 def __deleteObj(theObj):
104 """ Delete object from a Study
107 theObj (GEOM or SMESH object): object to remove from the study
109 aStudy = salome.myStudy
110 aStudyBuilder = aStudy.NewBuilder()
111 SO = aStudy.FindObjectIOR(aStudy.ConvertObjectToIOR(theObj))
113 aStudyBuilder.RemoveObjectWithChildren(SO)
116 def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True,
118 """ Create a dual of the mesh in input_file into output_file
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
126 # Import mesh from file
127 mesh = salome.orb.string_to_object(mesh_ior)
129 raise Exception("Could not find mesh using id: ", mesh_ior)
131 shape = mesh.GetShapeToMesh()
134 faces = geompy.SubShapeAll(shape, geompy.ShapeType["FACE"])
135 faces_ids = geompy.GetSubShapesIDs(shape, faces)
137 # Create group with each face
138 # so that we don't need GetFaceNearPoint to get the face to project the
141 id2face_edges_ids = {}
143 for face, face_id in zip(faces, faces_ids):
144 gr_mesh = mesh.CreateGroupFromGEOM(SMESH.FACE,
145 '%sface_%i'%(__prefix, face_id),
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
155 edges = geompy.SubShapeAll(shape, geompy.ShapeType["EDGE"])
156 edges_ids = geompy.GetSubShapesIDs(shape, edges)
159 for edge, edge_id in zip(edges, edges_ids):
160 gr_mesh = mesh.CreateGroupFromGEOM(SMESH.EDGE,
161 '%sedge_%i'%(__prefix, edge_id),
163 id2edge[edge_id] = edge
164 mesh_groups.append(gr_mesh)
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
173 tetras = mc.MEDCoupling1SGTUMesh(tetras)
175 # Create the polyhedra from the tetrahedra (main goal of this function)
176 polyh = tetras.computeDualMesh()
178 ## Adding skin + transfering groups on faces from tetras mesh
179 mesh2d = polyh.buildUnstructured().computeSkin()
180 mesh2d.setName(mesh_name)
182 polyh_coords = polyh.getCoords()
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)
192 logger.debug("Transferring group: "+ grp_name)
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)
198 if id_grp_poly is not None and grp_name[:4] == __prefix and adapt_to_shape:
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]
204 # for each face, get the edges bounding it
205 grp_poly = mesh2d[id_grp_poly]
206 mesh1d = grp_poly.computeSkin()
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)
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)
216 # project new nodes on its geom edge
217 # (if the group exists on this edge and it has not already been
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,
225 treated_edges.append(edge_id)
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)
231 # project new nodes on its geom face
232 __projectNodeOnSubshape(nodes_added_on_tri, face, polyh_coords)
234 # add the group to write it
235 mc_groups.append(id_grp_poly)
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)
245 for group in mc_groups:
246 myfile.addGroup(-1, group)
248 logger.debug("Writing dual mesh in: "+output_file)
249 myfile.write(output_file, 2)
252 # delete temporary groups
253 for grp in mesh_groups: