From aedb674c28199faebd54341cc71ad581f19971bc Mon Sep 17 00:00:00 2001 From: Christophe Bourcier Date: Mon, 14 Nov 2022 13:59:56 +0100 Subject: [PATCH] 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 Also: - Code cleaning in dual mesh - Remove tmp folder of dual mesh - Fix for Windows (thanks to Nabil) --- src/SMESH_I/SMESH_Gen_i.cxx | 9 +- src/SMESH_SWIG/smeshBuilder.py | 1 - src/SMESH_SWIG/smesh_tools.py | 261 ++++++++++++++++++++------- test/SMESH_create_dual_mesh_adapt.py | 2 + test/SMESH_create_dual_mesh_tpipe.py | 53 ++++-- 5 files changed, 241 insertions(+), 85 deletions(-) diff --git a/src/SMESH_I/SMESH_Gen_i.cxx b/src/SMESH_I/SMESH_Gen_i.cxx index 5c72ecf45..32265c514 100644 --- a/src/SMESH_I/SMESH_Gen_i.cxx +++ b/src/SMESH_I/SMESH_Gen_i.cxx @@ -2859,14 +2859,15 @@ SMESH::SMESH_Mesh_ptr SMESH_Gen_i::CreateDualMesh(SMESH::SMESH_IDSource_ptr mesh ats = "False"; std::string cmd="import salome.smesh.smesh_tools as smt\n"; - cmd +="smt.smesh_create_dual_mesh(\"" + mesh_ior + "\", \"" + + cmd +="smt.smesh_create_dual_mesh(\"" + mesh_ior + "\", r\"" + dual_mesh_file.string() + "\", mesh_name=\"" + mesh_name + "\", adapt_to_shape=" + ats + ")"; MESSAGE(cmd); PyObject *py_main = PyImport_AddModule("__main__"); PyObject *py_dict = PyModule_GetDict(py_main); + PyObject *local_dict = PyDict_New(); - PyRun_String(cmd.c_str(), Py_file_input, py_dict, py_dict); + PyRun_String(cmd.c_str(), Py_file_input, py_dict, local_dict); if (PyErr_Occurred()) { // Restrieving python error @@ -2931,6 +2932,10 @@ SMESH::SMESH_Mesh_ptr SMESH_Gen_i::CreateDualMesh(SMESH::SMESH_IDSource_ptr mesh MESSAGE("Dump of groups"); SMESH::ListOfGroups_var groups = newMesh->GetGroups(); +#ifndef _DEBUG_ + fs::remove_all(tmp_folder); +#endif + return newMesh._retn(); } diff --git a/src/SMESH_SWIG/smeshBuilder.py b/src/SMESH_SWIG/smeshBuilder.py index ec65b9fa9..972a4507c 100644 --- a/src/SMESH_SWIG/smeshBuilder.py +++ b/src/SMESH_SWIG/smeshBuilder.py @@ -804,7 +804,6 @@ class smeshBuilder( SMESH._objref_SMESH_Gen, object ): """ if isinstance( mesh, Mesh ): mesh = mesh.GetMesh() - print("calling createdualmesh from Python") dualMesh = SMESH._objref_SMESH_Gen.CreateDualMesh(self, mesh, meshName, adaptToShape) return Mesh(self, self.geompyD, dualMesh) diff --git a/src/SMESH_SWIG/smesh_tools.py b/src/SMESH_SWIG/smesh_tools.py index 24ccf390f..19d9fb6fd 100644 --- a/src/SMESH_SWIG/smesh_tools.py +++ b/src/SMESH_SWIG/smesh_tools.py @@ -1,12 +1,7 @@ #!/usr/bin/env python3 -import sys import salome import medcoupling as mc -from math import pi -import numpy as np - -#salome.salome_init() import GEOM from salome.geom import geomBuilder @@ -20,14 +15,113 @@ smesh = smeshBuilder.New() from salome.kernel.logger import Logger logger = Logger("salome.smesh.smesh_tools") -logger.setLevel("DEBUG") +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 -def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name="MESH"): + 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(mesh_ior) @@ -36,11 +130,38 @@ 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 +171,19 @@ 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() + 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 +191,64 @@ 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] - - # 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()) - - myfile.addGroup(-1, id_grp_poly) - - # 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] ) - - 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 + # 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: + # 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) 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 grp in mesh_groups: + __deleteObj(grp) diff --git a/test/SMESH_create_dual_mesh_adapt.py b/test/SMESH_create_dual_mesh_adapt.py index ccabb5aa5..e24cafbde 100644 --- a/test/SMESH_create_dual_mesh_adapt.py +++ b/test/SMESH_create_dual_mesh_adapt.py @@ -85,6 +85,8 @@ dual_Mesh_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_1', True) #Comparing volumes dual_volume = dual_Mesh_1.GetVolume() dual_raw_volume = dual_Mesh_raw_1.GetVolume() +tetra_volume = Mesh_1.GetVolume() +print("tetra_volume: ", tetra_volume) print("dual_volume: ", dual_volume) print("dual_raw_volume: ", dual_raw_volume) diff --git a/test/SMESH_create_dual_mesh_tpipe.py b/test/SMESH_create_dual_mesh_tpipe.py index ebb36e4d2..28bc26a4f 100644 --- a/test/SMESH_create_dual_mesh_tpipe.py +++ b/test/SMESH_create_dual_mesh_tpipe.py @@ -79,20 +79,24 @@ Mesh_1 = smesh.Mesh(tpipe, "Mesh_coarse") Mesh_1.Triangle(algo=smeshBuilder.NETGEN_1D2D) -Mesh_1.Tetrahedron(algo=smeshBuilder.NETGEN_3D) +algo3d = Mesh_1.Tetrahedron(algo=smeshBuilder.NETGEN_3D) +algo3d.MaxElementVolume(0.0002) + 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) +tetra_volume = smesh.GetVolume(Mesh_1) shape_volume = geompy.BasicProperties(tpipe)[2] -assert abs(mesh_volume-shape_volume)/shape_volume < 0.05 +assert abs(tetra_volume-shape_volume)/shape_volume < 0.05 dual_Mesh_raw_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_raw_1', False) @@ -102,7 +106,7 @@ assert dual_Mesh_raw_1.NbPolyhedrons() == dual_Mesh_raw_1.NbVolumes() # Check dual mesh volume dual_raw_volume = dual_Mesh_raw_1.GetVolume() -assert abs(dual_raw_volume-shape_volume)/shape_volume < 0.11 +assert abs(dual_raw_volume-shape_volume)/shape_volume < 0.14 # Check groups dual_Mesh_raw_groups = dual_Mesh_raw_1.GetGroups() @@ -110,25 +114,44 @@ 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("shape_volume: ", shape_volume) +print("tetra_volume: ", tetra_volume) +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.14 + +# 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() -- 2.30.2