From: Christophe Bourcier Date: Mon, 14 Nov 2022 12:59:56 +0000 (+0100) Subject: Improve projection of dual mesh nodes: X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=21ac1ad91f0f1fb7eb7da3e90e2cb8ba0e2672c9;p=modules%2Fsmesh.git Improve projection of dual mesh nodes: - first project nodes on edges, then project on faces - keep the relation from the sub-shapes to the mesh by groups --- diff --git a/src/SMESH_SWIG/smesh_tools.py b/src/SMESH_SWIG/smesh_tools.py index 24ccf390f..563d10ba8 100644 --- a/src/SMESH_SWIG/smesh_tools.py +++ b/src/SMESH_SWIG/smesh_tools.py @@ -22,6 +22,94 @@ from salome.kernel.logger import Logger logger = Logger("salome.smesh.smesh_tools") logger.setLevel("DEBUG") +# 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 @@ -36,11 +124,33 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name shape = mesh.GetShapeToMesh() - # Creating output file - logger.debug("Creating file with mesh: "+mesh_name) - myfile = mc.MEDFileUMesh() - myfile.setName(mesh_name) + 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) @@ -50,14 +160,21 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name # 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() ## Adding skin + transfering groups on faces from tetras mesh mesh2d = polyh.buildUnstructured().computeSkin() mesh2d.setName(mesh_name) - myfile.setMeshAtLevel(-1, mesh2d) + polyh_coords = polyh.getCoords() + + edges_group_names = mc_mesh_file.getGroupsOnSpecifiedLev(-2) + 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": @@ -65,63 +182,59 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name continue logger.debug("Transferring group: "+ grp_name) - grp_tria = mc_mesh_file.getGroup(-1, grp_name) - # Retrieve the nodes in group - grp_nodes = grp_tria.computeFetchedNodeIds() - # Find all the cells lying on one of the nodes - id_grp_poly = mesh2d.getCellIdsLyingOnNodes(grp_nodes, False) - - grp_poly = mesh2d[id_grp_poly] + # 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) - # find the extra face cells, on the border of the group (lying on nodes, but outside the group) - id_poly_border = grp_poly.findCellIdsOnBoundary() + if id_grp_poly is not None and grp_name[:4] == __prefix: + # This group is on a specific geom face + face_id = grp_name.split("_")[-1] + face_id = int(face_id) + face = id2face[face_id] - # ids of the cells in grp_poly - id_poly=mc.DataArrayInt64.New(grp_poly.getNumberOfCells(), 1) - id_poly.iota() + # for each face, get the edges bounding it + grp_poly = mesh2d[id_grp_poly] + mesh1d = grp_poly.computeSkin() - # cells that are really in the group - id_to_keep = id_poly.buildSubstraction(id_poly_border) + face_edges_id = id2face_edges_ids[face_id] + for edge_id in face_edges_id: + grp_seg_name = "%sedge_%i"%(__prefix, edge_id) - id_grp_poly = id_grp_poly[id_to_keep] - id_grp_poly.setName(grp_name.strip()) + # 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) - myfile.addGroup(-1, id_grp_poly) + # 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) - # Getting list of new points added on the skin - 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] ) + # remove these nodes from the nodes to project on face + nodes_added_on_tri = nodes_added_on_tri.buildSubstraction(nodes_added_on_segs) - if adapt_to_shape: - logger.debug("Adapting to shape") - ptsAddedCoo = ptsAddedMesh.getCoords() - ptsAddedCooModified = ptsAddedCoo[:] - - # Matching faces with their ids - faces = geompy.ExtractShapes(shape, geompy.ShapeType["FACE"], True) - id2face = {} - for face in faces: - id2face[face.GetSubShapeIndices()[0]] = face - - ## Projecting each points added by the dual mesh on the surface it is - # associated with - for i, tup in enumerate(ptsAddedCooModified): - vertex = geompy.MakeVertex(*tuple(tup)) - shapes = geompy.GetShapesNearPoint(shape, vertex, - geompy.ShapeType["FACE"]) - prj = geompy.MakeProjection(vertex, - id2face[shapes.GetSubShapeIndices()[0]]) - new_coor = geompy.PointCoordinates(prj) - ptsAddedCooModified[i] = new_coor - - polyh.getCoords()[ptsAdded] = ptsAddedCooModified + # 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) myfile.setMeshAtLevel(0, polyh) + myfile.setMeshAtLevel(-1, mesh2d) + + for group in mc_groups: + myfile.addGroup(-1, group) - logger.debug("Writting dual mesh in :"+output_file) + logger.debug("Writing dual mesh in: "+output_file) myfile.write(output_file, 2) + + if adapt_to_shape: + # delete temporary groups + for gr in mesh_groups: + __deleteObj(gr) diff --git a/test/SMESH_create_dual_mesh_tpipe.py b/test/SMESH_create_dual_mesh_tpipe.py index ebb36e4d2..3d121876d 100644 --- a/test/SMESH_create_dual_mesh_tpipe.py +++ b/test/SMESH_create_dual_mesh_tpipe.py @@ -83,11 +83,13 @@ Mesh_1.Tetrahedron(algo=smeshBuilder.NETGEN_3D) isDone = Mesh_1.Compute() # Create groups +d_geom_groups = {} d_mesh_groups = {} for geom_group in geom_groups: gr = Mesh_1.Group(geom_group) gr_name = gr.GetName() d_mesh_groups[gr_name] = gr + d_geom_groups[gr_name] = geom_group # Check tetra mesh volume mesh_volume = smesh.GetVolume(Mesh_1) @@ -110,25 +112,42 @@ dual_Mesh_raw_group_names = dual_Mesh_raw_1.GetGroupNames() assert len(dual_Mesh_raw_groups) == 4 -for gr_dual, gr_name in zip(dual_Mesh_raw_groups, dual_Mesh_raw_group_names): +## Create dual mesh with projection on faces +dual_Mesh_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_1', True) + +# Check dual mesh volume +dual_volume = dual_Mesh_1.GetVolume() +print("dual_volume: ", dual_volume) +print("dual_raw_volume: ", dual_raw_volume) +assert (dual_volume >= dual_raw_volume) + +assert abs(dual_volume-shape_volume)/shape_volume < 0.11 + +# Check groups +dual_Mesh_groups = dual_Mesh_1.GetGroups() +dual_Mesh_group_names = dual_Mesh_1.GetGroupNames() +d_mesh_dual_groups = {} +for gr_name, gr in zip(dual_Mesh_group_names, dual_Mesh_groups): + gr_name = gr_name.strip() + d_mesh_dual_groups[gr_name] = gr + +for gr_dual_raw, gr_name in zip(dual_Mesh_raw_groups, dual_Mesh_raw_group_names): gr_name = gr_name.strip() gr_tri = d_mesh_groups[gr_name] + gr_dual = d_mesh_dual_groups[gr_name] + gr_geom = d_geom_groups[gr_name] + area_gr_geom = geompy.BasicProperties(gr_geom)[1] area_gr_tri = smesh.GetArea(gr_tri) + area_gr_dual_raw = smesh.GetArea(gr_dual_raw) area_gr_dual = smesh.GetArea(gr_dual) print(gr_name) + print("Area geom: ", area_gr_geom) print("Area tri: ", area_gr_tri) + print("Area dual raw:", area_gr_dual_raw) print("Area dual:", area_gr_dual) - assert abs(area_gr_tri-area_gr_dual)/area_gr_tri < 1e-3 - -#dual_Mesh_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_1', True) - - -#dual_volume = dual_Mesh_1.GetVolume() -#dual_raw_volume = dual_Mesh_raw_1.GetVolume() -#print("dual_volume: ", dual_volume) -#print("dual_raw_volume: ", dual_raw_volume) + assert abs(area_gr_geom-area_gr_dual)/area_gr_geom < 0.015 + assert abs(area_gr_tri-area_gr_dual_raw)/area_gr_tri < 1e-3 -#assert (dual_volume >= dual_raw_volume) if salome.sg.hasDesktop(): salome.sg.updateObjBrowser()