test/SMESH_create_dual_mesh_tpipe.py

   1 #!/usr/bin/env python3
   2
   3 import sys
   4 import salome
   5 import math
   6
   7 salome.salome_init()
   8
   9 import GEOM
  10 from salome.geom import geomBuilder
  11
  12 geompy = geomBuilder.New()
  13
  14 # first cylinder
  15 r1 = 0.5
  16 h1 = 5
  17
  18 # second cylinder
  19 r2 = 0.3
  20 h2 = 3
  21
  22 length_piquage = 1.5
  23
  24 O = geompy.MakeVertex(0, 0, 0)
  25 OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
  26 OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
  27 OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
  28 geompy.addToStudy( O, 'O' )
  29 geompy.addToStudy( OX, 'OX' )
  30 geompy.addToStudy( OY, 'OY' )
  31 geompy.addToStudy( OZ, 'OZ' )
  32
  33 Cylinder_1 = geompy.MakeCylinderRH(r1, h1)
  34 Cylinder_2 = geompy.MakeCylinderRH(r2, h2)
  35 Rotation_1 = geompy.MakeRotation(Cylinder_2, OY, -90*math.pi/180.0)
  36 Translation_1 = geompy.MakeTranslation(Rotation_1, 0, 0, length_piquage)
  37
  38 tpipe = geompy.MakeFuseList([Cylinder_1, Translation_1], True, True)
  39 geompy.addToStudy( tpipe, 'tpipe' )
  40
  41 Inlet_z = geompy.GetFaceNearPoint(tpipe, O)
  42 geompy.addToStudyInFather( tpipe, Inlet_z, 'Inlet_z' )
  43
  44 p_inlet_x = geompy.MakeVertex(-h2, 0, length_piquage)
  45 Inlet_x = geompy.GetFaceNearPoint(tpipe, p_inlet_x)
  46 geompy.addToStudyInFather( tpipe, Inlet_x, 'Inlet_x' )
  47
  48 p_outlet = geompy.MakeVertex(0, 0, h1)
  49 Outlet = geompy.GetFaceNearPoint(tpipe, p_outlet)
  50 geompy.addToStudyInFather( tpipe, Outlet, 'Outlet' )
  51
  52 Wall = geompy.CreateGroup(tpipe, geompy.ShapeType["FACE"])
  53 faces = geompy.SubShapeAll(tpipe, geompy.ShapeType["FACE"])
  54 geompy.UnionList(Wall, faces)
  55 geompy.DifferenceList(Wall, [Inlet_x, Inlet_z, Outlet])
  56 geompy.addToStudyInFather( tpipe, Wall, 'Wall' )
  57
  58 p_corner = geompy.MakeVertex(-r2, 0, length_piquage+r2)
  59 corner = geompy.GetVertexNearPoint(tpipe, p_corner)
  60 geompy.addToStudyInFather( tpipe, corner, 'corner' )
  61
  62 geom_groups = [Inlet_x, Inlet_z, Outlet, Wall]
  63
  64 volumeGEOM = geompy.BasicProperties(tpipe)[2]
  65
  66 ###
  67 ### SMESH component
  68 ###
  69
  70 import  SMESH, SALOMEDS
  71 from salome.smesh import smeshBuilder
  72
  73 from salome.StdMeshers import StdMeshersBuilder
  74
  75 smesh = smeshBuilder.New()
  76
  77 # Coarse mesh with default hypothesis
  78 Mesh_1 = smesh.Mesh(tpipe, "Mesh_coarse")
  79
  80 Mesh_1.Triangle(algo=smeshBuilder.NETGEN_1D2D)
  81
  82 algo3d = Mesh_1.Tetrahedron(algo=smeshBuilder.NETGEN_3D)
  83 algo3d.MaxElementVolume(0.0002)
  84
  85 isDone = Mesh_1.Compute()
  86
  87 # Create groups
  88 d_geom_groups = {}
  89 d_mesh_groups = {}
  90 for geom_group in geom_groups:
  91   gr = Mesh_1.Group(geom_group)
  92   gr_name = gr.GetName()
  93   d_mesh_groups[gr_name] = gr
  94   d_geom_groups[gr_name] = geom_group
  95
  96 # Check tetra mesh volume
  97 tetra_volume = smesh.GetVolume(Mesh_1)
  98 shape_volume = geompy.BasicProperties(tpipe)[2]
  99 assert abs(tetra_volume-shape_volume)/shape_volume < 0.05
 100
 101 dual_Mesh_raw_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_raw_1', False)
 102
 103 # check polyhedrons
 104 assert dual_Mesh_raw_1.NbPolyhedrons() > 0
 105 assert dual_Mesh_raw_1.NbPolyhedrons() == dual_Mesh_raw_1.NbVolumes()
 106
 107 # Check dual mesh volume
 108 dual_raw_volume = dual_Mesh_raw_1.GetVolume()
 109 assert abs(dual_raw_volume-shape_volume)/shape_volume < 0.14
 110
 111 # Check groups
 112 dual_Mesh_raw_groups = dual_Mesh_raw_1.GetGroups()
 113 dual_Mesh_raw_group_names = dual_Mesh_raw_1.GetGroupNames()
 114
 115 assert len(dual_Mesh_raw_groups) == 4
 116
 117 ## Create dual mesh with projection on faces
 118 dual_Mesh_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_1', True)
 119
 120 # Check dual mesh volume
 121 dual_volume = dual_Mesh_1.GetVolume()
 122 print("shape_volume: ", shape_volume)
 123 print("tetra_volume: ", tetra_volume)
 124 print("dual_volume: ", dual_volume)
 125 print("dual_raw_volume: ", dual_raw_volume)
 126 assert (dual_volume >= dual_raw_volume)
 127
 128 assert abs(dual_volume-shape_volume)/shape_volume < 0.14
 129
 130 # Check groups
 131 dual_Mesh_groups = dual_Mesh_1.GetGroups()
 132 dual_Mesh_group_names = dual_Mesh_1.GetGroupNames()
 133 d_mesh_dual_groups = {}
 134 for gr_name, gr in zip(dual_Mesh_group_names, dual_Mesh_groups):
 135     gr_name = gr_name.strip()
 136     d_mesh_dual_groups[gr_name] = gr
 137
 138 for gr_dual_raw, gr_name in zip(dual_Mesh_raw_groups, dual_Mesh_raw_group_names):
 139   gr_name = gr_name.strip()
 140   gr_tri = d_mesh_groups[gr_name]
 141   gr_dual = d_mesh_dual_groups[gr_name]
 142   gr_geom = d_geom_groups[gr_name]
 143   area_gr_geom = geompy.BasicProperties(gr_geom)[1]
 144   area_gr_tri = smesh.GetArea(gr_tri)
 145   area_gr_dual_raw = smesh.GetArea(gr_dual_raw)
 146   area_gr_dual = smesh.GetArea(gr_dual)
 147   print(gr_name)
 148   print("Area geom: ", area_gr_geom)
 149   print("Area tri: ", area_gr_tri)
 150   print("Area dual raw:", area_gr_dual_raw)
 151   print("Area dual:", area_gr_dual)
 152   assert abs(area_gr_geom-area_gr_dual)/area_gr_geom < 0.015
 153   assert abs(area_gr_tri-area_gr_dual_raw)/area_gr_tri < 1e-3
 154
 155
 156 if salome.sg.hasDesktop():
 157   salome.sg.updateObjBrowser()