#!/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
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)
+ """ 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
+ 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.
+ 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))
+ 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()
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))
+ 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)
+ # 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 = mc.DataArrayInt64.New(grp_poly.getNumberOfCells(), 1)
id_poly.iota()
# cells that are really in the group
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
+ """ 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.
+ 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()
aStudyBuilder.RemoveObjectWithChildren(SO)
pass
-def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name="MESH"):
+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)
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
+ # 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)
+ 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
id2edge = {}
for edge, edge_id in zip(edges, edges_ids):
- gr_mesh = mesh.CreateGroupFromGEOM(SMESH.EDGE, '%sedge_%i'%(__prefix, edge_id), edge)
+ gr_mesh = mesh.CreateGroupFromGEOM(SMESH.EDGE,
+ '%sedge_%i'%(__prefix, edge_id),
+ edge)
id2edge[edge_id] = edge
mesh_groups.append(gr_mesh)
polyh_coords = polyh.getCoords()
- edges_group_names = mc_mesh_file.getGroupsOnSpecifiedLev(-2)
-
treated_edges = []
mc_groups = []
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)
+ 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
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)
+ 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 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)
+ __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)
+ 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)
if adapt_to_shape:
# delete temporary groups
- for gr in mesh_groups:
- __deleteObj(gr)
+ for grp in mesh_groups:
+ __deleteObj(grp)
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[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)
# 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.12
# Check groups
dual_Mesh_raw_groups = dual_Mesh_raw_1.GetGroups()
# 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.11
+assert abs(dual_volume-shape_volume)/shape_volume < 0.12
# Check groups
dual_Mesh_groups = dual_Mesh_1.GetGroups()
