doc/examples/viewing_meshes_ex03.py

   1 #!/usr/bin/env python
   2
   3 ###
   4 ### This file shows how to get the information that are displayed when using Mesh Information for elements/nodes
   5 ###
   6
   7 import sys
   8 import salome
   9
  10 salome.salome_init()
  11 import salome_notebook
  12 notebook = salome_notebook.NoteBook()
  13 sys.path.insert(0, r'/local00/home/B61570/work_in_progress/mesh_info')
  14
  15 ###
  16 ### GEOM component
  17 ###
  18
  19 import GEOM
  20 from salome.geom import geomBuilder
  21 import math
  22 import SALOMEDS
  23
  24
  25 geompy = geomBuilder.New()
  26
  27 Box_1 = geompy.MakeBoxDXDYDZ(200, 200, 200)
  28 geompy.addToStudy( Box_1, 'Box_1' )
  29 bottom = geompy.CreateGroup(Box_1, geompy.ShapeType["FACE"])
  30 geompy.UnionIDs(bottom, [31])
  31
  32 ###
  33 ### SMESH component
  34 ###
  35
  36 import  SMESH, SALOMEDS
  37 from salome.smesh import smeshBuilder
  38
  39 smesh = smeshBuilder.New()
  40
  41
  42 ## Tetra
  43 NETGEN_3D_Parameters_1 = smesh.CreateHypothesisByAverageLength( 'NETGEN_Parameters', 'NETGENEngine', 25, 0 )
  44 Mesh_tetra = smesh.Mesh(Box_1,'Mesh_tetra')
  45 status = Mesh_tetra.AddHypothesis( Box_1, NETGEN_3D_Parameters_1 )
  46 NETGEN_1D_2D_3D = Mesh_tetra.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
  47 isDone = Mesh_tetra.Compute()
  48 if not isDone:
  49     raise ("Could not compute mesh: "+Mesh_tetra.GetName())
  50
  51 ## Tetra
  52 Mesh_quadratic = smesh.Mesh(Box_1,'Mesh_quadratic')
  53 NETGEN_1D_2D_3D_1 = Mesh_quadratic.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
  54 NETGEN_3D_Parameters_2 = NETGEN_1D_2D_3D_1.Parameters()
  55 NETGEN_3D_Parameters_2.SetMaxSize( 34.641 )
  56 NETGEN_3D_Parameters_2.SetMinSize( 0.34641 )
  57 NETGEN_3D_Parameters_2.SetSecondOrder( 1 )
  58 isDone = Mesh_quadratic.Compute()
  59 if not isDone:
  60     raise ("Could not compute mesh: "+Mesh_quadratic.GetName())
  61
  62
  63 # Hexa mesh
  64 Mesh_hexa = smesh.Mesh(Box_1,'Mesh_hexa')
  65 Regular_1D = Mesh_hexa.Segment()
  66 Number_of_Segments_1 = Regular_1D.NumberOfSegments(15)
  67 Quadrangle_2D = Mesh_hexa.Quadrangle(algo=smeshBuilder.QUADRANGLE)
  68 Hexa_3D = Mesh_hexa.Hexahedron(algo=smeshBuilder.Hexa)
  69 isDone = Mesh_hexa.Compute()
  70 if not isDone:
  71     raise ("Could not compute mesh: "+Mesh_hexa.GetName())
  72
  73 # Poly Mesh
  74 Mesh_poly = smesh.CreateDualMesh(Mesh_tetra, 'dual_Mesh_1', True)
  75
  76 # Prism mesh
  77
  78 Mesh_prism = smesh.Mesh(Box_1,'Mesh_prism')
  79 Regular_1D_1 = Mesh_prism.Segment()
  80 Number_of_Segments_2 = Regular_1D_1.NumberOfSegments(15)
  81 NETGEN_1D_2D = Mesh_prism.Triangle(algo=smeshBuilder.NETGEN_1D2D,geom=bottom)
  82 NETGEN_2D_Parameters_1 = NETGEN_1D_2D.Parameters()
  83 NETGEN_2D_Parameters_1.SetMaxSize(35)
  84 NETGEN_2D_Parameters_1.SetMinSize(0.3)
  85 Prism_3D = Mesh_prism.Prism()
  86 isDone = Mesh_prism.Compute()
  87 if not isDone:
  88     raise ("Could not compute mesh: "+Mesh_prism.GetName())
  89
  90 # Pyramid mesh
  91 Mesh_pyramids = smesh.Mesh(Box_1,'Mesh_pyramids')
  92 Regular_1D_2 = Mesh_pyramids.Segment()
  93 Number_of_Segments_3 = Regular_1D_2.NumberOfSegments(15)
  94 Quadrangle_2D_1 = Mesh_pyramids.Quadrangle(algo=smeshBuilder.QUADRANGLE)
  95 NETGEN_3D = Mesh_pyramids.Tetrahedron()
  96 bottom_1 = Mesh_pyramids.GroupOnGeom(bottom,'bottom',SMESH.FACE)
  97 isDone = Mesh_pyramids.Compute()
  98 if not isDone:
  99     raise ("Could not compute mesh: "+Mesh_pyramids.GetName())
 100
 101
 102 ## Set names of Mesh objects
 103 smesh.SetName(NETGEN_1D_2D_3D.GetAlgorithm(), 'NETGEN 1D-2D-3D')
 104 smesh.SetName(NETGEN_3D_Parameters_1, 'NETGEN 3D Parameters_1')
 105 smesh.SetName(Mesh_tetra.GetMesh(), 'Mesh_tetra')
 106 smesh.SetName(Mesh_hexa.GetMesh(), 'Mesh_hexa')
 107
 108 if salome.sg.hasDesktop():
 109   salome.sg.updateObjBrowser()
 110
 111 # Look in SMESH_GUI/SMESHGUI_MeshInfo.cxx +1666 for list of what is
 112
 113 def face_info(mesh, elem_id):
 114     """
 115     Print equivalent of Mesh Information for a face
 116     """
 117     elem_type = mesh.GetElementGeomType(elem_id)
 118
 119     conn = mesh.GetElemNodes(elem_id)
 120
 121     nb_nodes = len(conn)
 122
 123     position = mesh.GetElementPosition(elem_id)
 124     pos = f"{position.shapeType} #{position.shapeID}"
 125
 126     grav_center = mesh.BaryCenter(elem_id)
 127
 128     normal = mesh.GetFaceNormal(elem_id, normalized=True)
 129
 130     aspect_ratio = mesh.GetAspectRatio(elem_id)
 131     #aspect_ratio = mesh.FunctorValue(SMESH.FT_AspectRatio, elem_id, isElem=True)
 132
 133     warping = mesh.GetWarping(elem_id)
 134     #warping = mesh.FunctorValue(SMESH.FT_Warping, elem_id, isElem=True)
 135
 136     min_angle = mesh.GetMinimumAngle(elem_id)
 137     #min_angle = mesh.FunctorValue(SMESH.FT_MinimumAngle, elem_id, isElem=True)
 138
 139     taper = mesh.GetTaper(elem_id)
 140     #taper = mesh.FunctorValue(SMESH.FT_Taper, elem_id, isElem=True)
 141
 142     skew = mesh.GetSkew(elem_id)
 143     #skew = mesh.FunctorValue(SMESH.FT_Skew, elem_id, isElem=True)
 144
 145     area = mesh.GetArea(elem_id)
 146     #area = mesh.FunctorValue(SMESH.FT_Area, elem_id, isElem=True)
 147
 148     diameter = mesh.GetMaxElementLength(elem_id)
 149     #diameter = mesh.FunctorValue(SMESH.FT_MaxElementLength2D, elem_id, isElem=True)
 150
 151     min_length = mesh.FunctorValue(SMESH.FT_Length2D, elem_id, isElem=True)
 152
 153     string = f"""
 154 Id: {elem_id}
 155 Type: {elem_type}
 156 Nb Nodes: {nb_nodes}
 157 Connectivity: {conn}
 158 Position: {pos}
 159 Gravity center:
 160  - X: {grav_center[0]}
 161  - Y: {grav_center[1]}
 162  - Z: {grav_center[2]}
 163 Normal:
 164  - X: {normal[0]}
 165  - Y: {normal[1]}
 166  - Z: {normal[2]}
 167 Quality:
 168  - Aspect Ratio: {aspect_ratio}
 169  - Warping: {warping}
 170  - Minimum Angle: {min_angle}
 171  - Taper: {taper}
 172  - Skew: {skew}
 173  - Area: {area}
 174  - Element Diameter 2D: {diameter}
 175  - Minimum Edge Length: {min_length}
 176 """
 177     print(string)
 178
 179 def volume_info(mesh, elem_id):
 180     """
 181     Print equivalent of Mesh Information for a volume
 182     """
 183     elem_type = mesh.GetElementGeomType(elem_id)
 184
 185     if elem_type in [SMESH.Entity_Polyhedra, SMESH.Entity_Quad_Polyhedra]:
 186         iface = 1
 187         face_conn = [12]
 188         conn = []
 189         while face_conn != []:
 190             face_conn = mesh.GetElemFaceNodes(elem_id, iface)
 191             iface += 1
 192             conn.append(face_conn)
 193             nb_nodes = len(mesh.GetElemNodes(elem_id))
 194     else:
 195         conn = mesh.GetElemNodes(elem_id)
 196         nb_nodes = len(conn)
 197
 198
 199     position = mesh.GetElementPosition(elem_id)
 200     pos = f"{position.shapeType} #{position.shapeID}"
 201
 202     grav_center = mesh.BaryCenter(elem_id)
 203
 204     aspect_ratio = mesh.GetAspectRatio(elem_id)
 205     #aspect_ratio = mesh.FunctorValue(SMESH.FT_AspectRatio3D, elem_id, isElem=True)
 206
 207     volume = mesh.GetVolume(elem_id)
 208     #volume = mesh.FunctorValue(SMESH.FT_Volume3D, elem_id, isElem=True)
 209
 210     jacob = mesh.GetScaledJacobian(elem_id)
 211     #jacob = mesh.FunctorValue(SMESH.FT_ScaledJacobian, elem_id, isElem=True)
 212
 213     diameter = mesh.GetMaxElementLength(elem_id)
 214     #diameter = mesh.FunctorValue(SMESH.FT_MaxElementLength3D, elem_id, isElem=True)
 215
 216     min_length = mesh.FunctorValue(SMESH.FT_Length3D, elem_id, isElem=True)
 217
 218     string = f"""
 219 Id: {elem_id}
 220 Type: {elem_type}
 221 Nb Nodes: {nb_nodes}
 222 Connectivity: {conn}
 223 Position: {pos}
 224 Gravity center:
 225  - X: {grav_center[0]}
 226  - Y: {grav_center[1]}
 227  - Z: {grav_center[2]}
 228 Quality:
 229  - Aspect Ratio 3D: {aspect_ratio}
 230  - Volume: {volume}
 231  - Scaled Jacobian: {jacob}
 232  - Element Diameter 3D: {diameter}
 233  - Minimum Edge Length: {min_length}
 234 """
 235     print(string)
 236
 237 def node_info(mesh, node_id):
 238
 239     coord = mesh.GetNodeXYZ(node_id)
 240
 241     conn_edge = mesh.GetNodeInverseElements(node_id, SMESH.EDGE)
 242     conn_face = mesh.GetNodeInverseElements(node_id, SMESH.FACE)
 243     conn_vol = mesh.GetNodeInverseElements(node_id, SMESH.VOLUME)
 244
 245     position = mesh.GetNodePosition(node_id)
 246     pos = f"{position.shapeType} #{position.shapeID}"
 247
 248     vec = [None, None]
 249     vec[0:len(position.params)] = position.params
 250
 251     string = f"""
 252 Id: {node_id}
 253 Coordinates:
 254 - X: {coord[0]}
 255 - Y: {coord[1]}
 256 - Z: {coord[2]}
 257 Connectivity
 258 - Edges: {conn_edge}
 259 - Faces: {conn_face}
 260 - Volumes: {conn_vol}
 261 Position: {pos}
 262 - U: {vec[0]}
 263 - V: {vec[1]}
 264 """
 265     print(string)
 266
 267 ###
 268 # Volume
 269 ##
 270
 271 # Tetrahedron
 272 volume_info(Mesh_tetra, 3000)
 273 # Hexahedron
 274 volume_info(Mesh_hexa, 3000)
 275 # Polyhedron
 276 volume_info(Mesh_poly, 3000)
 277 # Prism
 278 volume_info(Mesh_prism, 1400)
 279 # Pyramids
 280 volume_info(Mesh_pyramids, 8176)
 281 # Quadratic tetra
 282 volume_info(Mesh_quadratic, 1180)
 283
 284 ###
 285 # Face
 286 ##
 287
 288 ## Triangle
 289 face_info(Mesh_tetra, 147)
 290 #Quadrangle
 291 face_info(Mesh_hexa, 1464)
 292 # Polygon
 293 face_info(Mesh_poly, 771)
 294 # Quadratic triangle
 295 face_info(Mesh_quadratic, 138)
 296
 297
 298 ###
 299 # Node
 300 ###
 301 # U & V
 302 node_info(Mesh_tetra, 152)
 303 # U
 304 node_info(Mesh_tetra, 32)
 305 # None
 306 node_info(Mesh_tetra, 2)