tests/layers_imprinting.py

   1 # -*- coding: utf-8 -*-
   2 # Copyright (C) 2017-2022  CEA/DEN, EDF R&D
   3 #
   4 # This library is free software; you can redistribute it and/or
   5 # modify it under the terms of the GNU Lesser General Public
   6 # License as published by the Free Software Foundation; either
   7 # version 2.1 of the License, or (at your option) any later version.
   8 #
   9 # This library is distributed in the hope that it will be useful,
  10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 # Lesser General Public License for more details.
  13 #
  14 # You should have received a copy of the GNU Lesser General Public
  15 # License along with this library; if not, write to the Free Software
  16 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 #
  18 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 #
  20
  21 import sys
  22 import salome
  23
  24 salome.salome_init()
  25
  26 ###
  27 ### GEOM component
  28 ###
  29
  30 import GEOM
  31 from salome.geom import geomBuilder
  32 import math
  33 import SALOMEDS
  34
  35
  36 geompy = geomBuilder.New()
  37
  38 # first cylinder
  39 r1 = 0.5
  40 h1 = 5
  41
  42 # second cylinder
  43 r2 = 0.3
  44 h2 = 3
  45
  46 length_piquage = 1.5
  47
  48 O = geompy.MakeVertex(0, 0, 0)
  49 OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
  50 OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
  51 OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
  52 geompy.addToStudy( O, 'O' )
  53 geompy.addToStudy( OX, 'OX' )
  54 geompy.addToStudy( OY, 'OY' )
  55 geompy.addToStudy( OZ, 'OZ' )
  56
  57 Cylinder_1 = geompy.MakeCylinderRH(r1, h1)
  58 Cylinder_2 = geompy.MakeCylinderRH(r2, h2)
  59 Rotation_1 = geompy.MakeRotation(Cylinder_2, OY, -90*math.pi/180.0)
  60 Translation_1 = geompy.MakeTranslation(Rotation_1, 0, 0, length_piquage)
  61
  62 piquage = geompy.MakeFuseList([Cylinder_1, Translation_1], True, True)
  63 geompy.addToStudy( piquage, 'piquage' )
  64
  65 Inlet_z = geompy.GetFaceNearPoint(piquage, O)
  66 geompy.addToStudyInFather( piquage, Inlet_z, 'Inlet_z' )
  67
  68 p_inlet_x = geompy.MakeVertex(-h2, 0, length_piquage)
  69 Inlet_x = geompy.GetFaceNearPoint(piquage, p_inlet_x)
  70 geompy.addToStudyInFather( piquage, Inlet_x, 'Inlet_x' )
  71
  72 p_outlet = geompy.MakeVertex(0, 0, h1)
  73 Outlet = geompy.GetFaceNearPoint(piquage, p_outlet)
  74 geompy.addToStudyInFather( piquage, Outlet, 'Outlet' )
  75
  76 Wall = geompy.CreateGroup(piquage, geompy.ShapeType["FACE"])
  77 faces = geompy.SubShapeAll(piquage, geompy.ShapeType["FACE"])
  78 geompy.UnionList(Wall, faces)
  79 geompy.DifferenceList(Wall, [Inlet_x, Inlet_z, Outlet])
  80 geompy.addToStudyInFather( piquage, Wall, 'Wall' )
  81
  82 p_corner = geompy.MakeVertex(-r2, 0, length_piquage+r2)
  83 corner = geompy.GetVertexNearPoint(piquage, p_corner)
  84 geompy.addToStudyInFather( piquage, corner, 'corner' )
  85
  86 geom_groups = [Inlet_x, Inlet_z, Outlet, Wall]
  87 d_geom_groups = {}
  88
  89 ###
  90 ### SMESH component
  91 ###
  92
  93 import  SMESH, SALOMEDS
  94 from salome.smesh import smeshBuilder
  95
  96 from salome.StdMeshers import StdMeshersBuilder
  97 from salome.HYBRIDPlugin import HYBRIDPluginBuilder
  98
  99 smesh = smeshBuilder.New()
 100
 101
 102 # Viscous layers without imprinting
 103 # =================================
 104
 105 Mesh_1 = smesh.Mesh(piquage, "Mesh_without_imprinting")
 106 NETGEN_2D = Mesh_1.Triangle(algo=smeshBuilder.NETGEN_1D2D)
 107 NETGEN_2D_Parameters = NETGEN_2D.Parameters()
 108
 109 NETGEN_2D_Parameters.SetMinSize( 0.01 )
 110 NETGEN_2D_Parameters.SetLocalSizeOnShape(corner, 0.01)
 111 NETGEN_2D_Parameters.SetFineness( 5 )
 112 NETGEN_2D_Parameters.SetGrowthRate( 0.1 )
 113 NETGEN_2D_Parameters.SetNbSegPerEdge( 2 )
 114 NETGEN_2D_Parameters.SetNbSegPerRadius( 3 )
 115
 116 wall_faces = geompy.SubShapeAll(Wall, geompy.ShapeType["FACE"])
 117 wall_ids = geompy.GetSubShapesIDs(piquage, wall_faces)
 118
 119 MG_Hybrid = Mesh_1.Tetrahedron(algo=smeshBuilder.HYBRID)
 120 MG_Hybrid_Parameters_1 = MG_Hybrid.Parameters()
 121 MG_Hybrid_Parameters_1.SetLayersOnAllWrap( 0 )
 122 MG_Hybrid_Parameters_1.SetElementGeneration( HYBRIDPluginBuilder.Generation_Tetra_Dominant )
 123 MG_Hybrid_Parameters_1.SetHeightFirstLayer( 0.01 )
 124 MG_Hybrid_Parameters_1.SetBoundaryLayersProgression( 1.1 )
 125 MG_Hybrid_Parameters_1.SetNbOfBoundaryLayers( 3 )
 126 MG_Hybrid_Parameters_1.SetFacesWithLayers( wall_ids )
 127
 128
 129 isDone = Mesh_1.Compute()
 130
 131 if not isDone:
 132   raise Exception("Error when computing Mesh_without_imprinting")
 133
 134 d_groups_1 = {}
 135
 136 for geom_group in geom_groups:
 137   name = geom_group.GetName()
 138   gr = Mesh_1.Group(geom_group)
 139   d_groups_1[name] = gr
 140   d_geom_groups[name] = geom_group
 141
 142 assert Mesh_1.NbQuadrangles() == 0
 143
 144 shape_volume = geompy.BasicProperties(piquage)[2]
 145
 146 # Compare whole mesh volume
 147 volume_error_1 = abs(shape_volume-Mesh_1.GetVolume())/shape_volume
 148
 149 assert volume_error_1 < 0.02
 150
 151 # Viscous layers with imprinting
 152 # ==============================
 153
 154 Mesh_2 = smesh.Mesh(piquage, "Mesh_with_imprinting")
 155
 156 NETGEN_2D_1_1 = Mesh_2.Triangle(algo=smeshBuilder.NETGEN_1D2D)
 157 status = Mesh_2.AddHypothesis(NETGEN_2D_Parameters)
 158
 159 Inlet_x_id, Inlet_z_id, Outlet_id = geompy.GetSubShapesIDs(piquage, [Inlet_x, Inlet_z, Outlet])
 160
 161 MG_Hybrid_2 = Mesh_2.Tetrahedron(algo=smeshBuilder.HYBRID)
 162 MG_Hybrid_Parameters_2 = MG_Hybrid_2.Parameters()
 163 MG_Hybrid_Parameters_2.SetLayersOnAllWrap( 0 )
 164 MG_Hybrid_Parameters_2.SetElementGeneration( HYBRIDPluginBuilder.Generation_Tetra_Dominant )
 165 MG_Hybrid_Parameters_2.SetHeightFirstLayer( 0.01 )
 166 MG_Hybrid_Parameters_2.SetBoundaryLayersProgression( 1.1 )
 167 MG_Hybrid_Parameters_2.SetNbOfBoundaryLayers( 3 )
 168 MG_Hybrid_Parameters_2.SetFacesWithLayers( wall_ids )
 169 MG_Hybrid_Parameters_2.SetFacesWithImprinting( [ Inlet_x_id, Inlet_z_id, Outlet_id ] )
 170
 171 isDone = Mesh_2.Compute()
 172
 173 if not isDone:
 174   raise Exception("Error when computing Mesh_with_imprinting")
 175
 176 assert Mesh_2.NbQuadrangles() > 0
 177
 178 d_groups_2 = {}
 179
 180 for geom_group in geom_groups:
 181   name = geom_group.GetName()
 182   gr = Mesh_2.Group(geom_group)
 183   d_groups_2[name] = gr
 184
 185 faces_imprinted = ["Inlet_x", "Inlet_z", "Outlet"]
 186
 187 # Check viscous layers with imprinting
 188 for name in faces_imprinted:
 189   geom_area = geompy.BasicProperties(d_geom_groups[name])[1]
 190   gr_1 = d_groups_1[name]
 191   gr_2 = d_groups_2[name]
 192   #assert gr_1.Size() > 0
 193   assert gr_2.Size() > 0
 194   # Nb of quadrangles is in 7th index of mesh info
 195   assert gr_2.GetMeshInfo()[7] > 0
 196
 197   # Compare mesh group and geom group
 198   area_error_1 = abs(geom_area-smesh.GetArea(gr_2))/geom_area
 199   assert area_error_1 < 0.025
 200   # Compare mesh group with imprinting and mesh group without imprinting
 201   area_error_2 = abs(smesh.GetArea(gr_1)-smesh.GetArea(gr_2))/smesh.GetArea(gr_1)
 202   assert area_error_2 < 1e-08
 203
 204 # Compare whole mesh volume
 205 mesh_2_volume = Mesh_2.GetVolume()
 206 volume_error_2 = abs(shape_volume-mesh_2_volume)/shape_volume
 207
 208 assert volume_error_2 < 0.02
 209
 210 # Viscous layers with imprinting set by groups
 211 # ============================================
 212
 213 Mesh_3 = smesh.Mesh(piquage, "Mesh_with_imprinting_set_by_groups")
 214
 215 NETGEN_2D_1_1 = Mesh_3.Triangle(algo=smeshBuilder.NETGEN_1D2D)
 216 status = Mesh_3.AddHypothesis(NETGEN_2D_Parameters)
 217
 218 MG_Hybrid_3 = Mesh_3.Tetrahedron(algo=smeshBuilder.HYBRID)
 219 MG_Hybrid_Parameters_3 = MG_Hybrid_3.Parameters()
 220 MG_Hybrid_Parameters_3.SetLayersOnAllWrap( 0 )
 221 MG_Hybrid_Parameters_3.SetElementGeneration( HYBRIDPluginBuilder.Generation_Tetra_Dominant )
 222 MG_Hybrid_Parameters_3.SetHeightFirstLayer( 0.01 )
 223 MG_Hybrid_Parameters_3.SetBoundaryLayersProgression( 1.1 )
 224 MG_Hybrid_Parameters_3.SetNbOfBoundaryLayers( 3 )
 225 MG_Hybrid_3.SetFacesWithLayers( Wall )
 226 MG_Hybrid_3.SetFacesWithImprinting( [ Inlet_x, Inlet_z, Outlet ] )
 227
 228 isDone = Mesh_3.Compute()
 229
 230 if not isDone:
 231   raise Exception("Error when computing Mesh_with_imprinting_set_by_groups")
 232
 233 d_groups_3 = {}
 234
 235 for geom_group in geom_groups:
 236   name = geom_group.GetName()
 237   gr = Mesh_3.Group(geom_group)
 238   d_groups_3[name] = gr
 239
 240 # Compare whole mesh volume
 241 mesh_3_volume = Mesh_3.GetVolume()
 242 volume_error_3 = abs(mesh_2_volume-mesh_3_volume)/mesh_2_volume
 243
 244 assert math.isclose(mesh_3_volume,mesh_2_volume,rel_tol=1e-7)
 245
 246 # Check viscous layers with imprinting
 247 for name in faces_imprinted:
 248   gr_2 = d_groups_2[name]
 249   gr_3 = d_groups_3[name]
 250   assert gr_3.Size() > 0
 251   # Nb of quadrangles is in 7th index of mesh info
 252   assert gr_2.GetMeshInfo()[7] == gr_3.GetMeshInfo()[7]
 253
 254   # Compare mesh group with imprinting set by ids and mesh group with imprinting set by geom group
 255   assert math.isclose( smesh.GetArea(gr_2), smesh.GetArea(gr_3))
 256