X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH_SWIG%2Fsmesh_tools.py;h=b964f55eabbac9e79592727ae2ea2fb3c13fc1cd;hp=38cd978e0748d8f879d968544707d91a557d0d71;hb=HEAD;hpb=08485b06fa309a40c3769a9113c2b88dd26e103d diff --git a/src/SMESH_SWIG/smesh_tools.py b/src/SMESH_SWIG/smesh_tools.py index 38cd978e0..b964f55ea 100644 --- a/src/SMESH_SWIG/smesh_tools.py +++ b/src/SMESH_SWIG/smesh_tools.py @@ -1,11 +1,7 @@ #!/usr/bin/env python3 -import sys import salome import medcoupling as mc -from math import pi - -#salome.salome_init() import GEOM from salome.geom import geomBuilder @@ -17,58 +13,242 @@ from salome.smesh import smeshBuilder smesh = smeshBuilder.New() -def create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name="MESH"): +from salome.kernel.logger import Logger +logger = Logger("salome.smesh.smesh_tools") +logger.setLevel("WARNING") + +# prefix for groups with internal usage +# i.e. used to transfer the faces and edges sub-shapes ids to the mesh +__prefix = "____" + +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) + + Args: + mc_mesh_file (MEDFileMesh): mesh on which to read the group + grp_name (string): name of the group to read + 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) + + Returns: + 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 + """ + try: + grp_tria = mc_mesh_file.getGroup(grp_level, grp_name) + except: + logger.debug("""No group found for %s at level %i. + It is normal behaviour for degenerated geom edge."""\ + %(grp_name, grp_level)) + return None, None + # Retrieve the nodes in group + orig_grp_nodes = grp_tria.computeFetchedNodeIds() + # Find all the cells lying on one of the nodes + id_grp_poly = mesh2d.getCellIdsLyingOnNodes(orig_grp_nodes, False) + + grp_poly = mesh2d[id_grp_poly] + if grp_poly.getNumberOfCells() == 0: + logger.debug("""No poly cell found for %s at level %i."""\ + %(grp_name, grp_level)) + return None, None + + # find the extra face cells, on the border of the group (lying on nodes, + # but outside the group) + id_poly_border = grp_poly.findCellIdsOnBoundary() + + # ids of the cells in grp_poly + id_poly = mc.DataArrayInt64.New(grp_poly.getNumberOfCells(), 1) + id_poly.iota() + + # cells that are really in the group + id_to_keep = id_poly.buildSubstraction(id_poly_border) + + id_grp_poly = id_grp_poly[id_to_keep] + + id_grp_poly.setName(grp_name.strip()) + + # get nodes added on this group + grp_poly = mesh2d[id_grp_poly] + grp_nodes_poly = grp_poly.computeFetchedNodeIds() + nodes_added = grp_nodes_poly.buildSubstraction(orig_grp_nodes) + + return id_grp_poly, nodes_added + +def __projectNodeOnSubshape(nodes, subshape, coords): + """ Project nodes on a sub-shape (face or edge) and update the mesh + coordinates + + Args: + nodes (DataArrayInt): nodes ids to project + 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. + """ + for i in nodes: + x, y, z = coords[i].getValues() + vertex = geompy.MakeVertex(x, y, z) + try: + prj = geompy.MakeProjection(vertex, subshape) + except: + logger.warning("Projection failed for %.5f %.5f %.5f but we continue with next node"%(x, y, z)) + continue + new_coor = geompy.PointCoordinates(prj) + # update its coordinates in the mesh + coords[i] = new_coor + pass + +def __deleteObj(theObj): + """ Delete object from a Study + + Args: + theObj (GEOM or SMESH object): object to remove from the study + """ + aStudy = salome.myStudy + aStudyBuilder = aStudy.NewBuilder() + SO = aStudy.FindObjectIOR(aStudy.ConvertObjectToIOR(theObj)) + if SO is not None: + aStudyBuilder.RemoveObjectWithChildren(SO) + pass + +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 Args: mesh_ior (string): corba Id of the Tetrahedron mesh 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 - # mesh = salome.orb.string_to_object(salome.salome_study.myStudy.FindObjectID(mesh_id).GetIOR()) mesh = salome.orb.string_to_object(mesh_ior) if not mesh: raise Exception("Could not find mesh using id: ", mesh_ior) shape = mesh.GetShapeToMesh() + if adapt_to_shape: + faces = geompy.SubShapeAll(shape, geompy.ShapeType["FACE"]) + faces_ids = geompy.GetSubShapesIDs(shape, faces) + + # Create group with each face + # so that we don't need GetFaceNearPoint to get the face to project the + # point to + id2face = {} + id2face_edges_ids = {} + mesh_groups = [] + for face, face_id in zip(faces, faces_ids): + gr_mesh = mesh.CreateGroupFromGEOM(SMESH.FACE, + '%sface_%i'%(__prefix, face_id), + face) + id2face[face_id] = face + mesh_groups.append(gr_mesh) + # get the edges bounding this face + # so that we can project the nodes on edges before nodes on faces + face_edges = geompy.SubShapeAll(face, geompy.ShapeType["EDGE"]) + face_edges_ids = geompy.GetSubShapesIDs(shape, face_edges) + id2face_edges_ids[face_id] = face_edges_ids + + edges = geompy.SubShapeAll(shape, geompy.ShapeType["EDGE"]) + edges_ids = geompy.GetSubShapesIDs(shape, edges) + + id2edge = {} + for edge, edge_id in zip(edges, edges_ids): + gr_mesh = mesh.CreateGroupFromGEOM(SMESH.EDGE, + '%sedge_%i'%(__prefix, edge_id), + edge) + id2edge[edge_id] = edge + mesh_groups.append(gr_mesh) + # We got a meshProxy so we need to convert pointer to MEDCoupling int_ptr = mesh.ExportMEDCoupling(True, True) dab = mc.FromPyIntPtrToDataArrayByte(int_ptr) - tetras = mc.MEDFileMesh.New(dab)[0] + mc_mesh_file = mc.MEDFileMesh.New(dab) + tetras = mc_mesh_file[0] # End of SMESH -> MEDCoupling part for dualmesh tetras = mc.MEDCoupling1SGTUMesh(tetras) + + # Create the polyhedra from the tetrahedra (main goal of this function) polyh = tetras.computeDualMesh() - skin = tetras.buildUnstructured().computeSkin() - skin_polyh = polyh.buildUnstructured().computeSkin() - allNodesOnSkinPolyh = skin_polyh.computeFetchedNodeIds() - allNodesOnSkin = skin.computeFetchedNodeIds() - ptsAdded = allNodesOnSkinPolyh.buildSubstraction(allNodesOnSkin) - ptsAddedMesh = mc.MEDCouplingUMesh.Build0DMeshFromCoords( skin_polyh.getCoords()[ptsAdded] ) - if adapt_to_shape: - ptsAddedCoo = ptsAddedMesh.getCoords() - ptsAddedCooModified = ptsAddedCoo[:] - - # We need the geometry for that - # TODO : Loop on faces identify points associated to which face - faces = geompy.ExtractShapes(shape, geompy.ShapeType["FACE"], True) - #assert( len(faces) == 1 ) - ## projection des points ajoutés par le dual sur la surface - #for i,tup in enumerate(ptsAddedCooModified): - # vertex = geompy.MakeVertex(*tuple(tup)) - # prj = geompy.MakeProjection(vertex, faces) - # newCoor = geompy.PointCoordinates( prj ) - # ptsAddedCooModified[i] = newCoor - ## assign coordinates with projected ones - #polyh.getCoords()[ptsAdded] = ptsAddedCooModified - - print("Writing dual mesh in ", output_file) + ## Adding skin + transfering groups on faces from tetras mesh + mesh2d = polyh.buildUnstructured().computeSkin() + mesh2d.setName(mesh_name) + + polyh_coords = polyh.getCoords() + + treated_edges = [] + + mc_groups = [] + for grp_name in mc_mesh_file.getGroupsOnSpecifiedLev(-1): + # This group is created by the export + if grp_name == "Group_Of_All_Faces": + logger.debug("Skipping group: "+ grp_name) + continue + logger.debug("Transferring group: "+ grp_name) + + # 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) + + if id_grp_poly is not None and grp_name[:4] == __prefix and adapt_to_shape: + # This group is on a specific geom face + face_id = grp_name.split("_")[-1] + face_id = int(face_id) + face = id2face[face_id] + + # for each face, get the edges bounding it + grp_poly = mesh2d[id_grp_poly] + mesh1d = grp_poly.computeSkin() + + face_edges_id = id2face_edges_ids[face_id] + for edge_id in face_edges_id: + grp_seg_name = "%sedge_%i"%(__prefix, edge_id) + + # 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) + + # project new nodes on its geom edge + # (if the group exists on this edge and it has not already been + # treated) + if id_grp_seg is not None: + if edge_id not in treated_edges: + edge = id2edge[edge_id] + __projectNodeOnSubshape(nodes_added_on_segs, + edge, polyh_coords) + else: + treated_edges.append(edge_id) + + # remove these nodes from the nodes to project on face + nodes_added_on_tri = \ + nodes_added_on_tri.buildSubstraction(nodes_added_on_segs) + + # project new nodes on its geom face + __projectNodeOnSubshape(nodes_added_on_tri, face, polyh_coords) + else: + # add the group to write it + mc_groups.append(id_grp_poly) + + # Creating output file + logger.debug("Creating file with mesh: "+mesh_name) + myfile = mc.MEDFileUMesh() + myfile.setName(mesh_name) polyh.setName(mesh_name) - polyh.write(output_file) - - + myfile.setMeshAtLevel(0, polyh) + myfile.setMeshAtLevel(-1, mesh2d) + for group in mc_groups: + myfile.addGroup(-1, group) + logger.debug("Writing dual mesh in: "+output_file) + myfile.write(output_file, 2) + if adapt_to_shape: + # delete temporary groups + for grp in mesh_groups: + __deleteObj(grp)