test/test_import_1D2D_with_tetras_and_pyras.py

   1 #!/usr/bin/env python
   2
   3 # test fix bos #33557 Bad pyramids generated
   4
   5 import sys
   6 import salome
   7
   8
   9 #nb_segs_right = 30
  10 #nb_segs_right = 100
  11 nb_segs_right = 150
  12 #nb_segs_right = 200
  13
  14 algo = "MG-Tetra"
  15 #algo = "Netgen"
  16 #algo = "GMSH"
  17
  18 algo_gmsh = "Delaunay"
  19 #algo_gmsh = "HXT"
  20
  21 salome.salome_init()
  22
  23
  24
  25 ###
  26 ### GEOM component
  27 ###
  28
  29 import GEOM
  30 from salome.geom import geomBuilder
  31 import math
  32 import SALOMEDS
  33
  34
  35 geompy = geomBuilder.New()
  36
  37 O = geompy.MakeVertex(0, 0, 0)
  38 OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
  39 OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
  40 OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
  41 Vertex_1 = geompy.MakeVertex(-0.62375, 0.0575, 0.02)
  42 Vertex_2 = geompy.MakeVertex(-0.62375, -0.0575, 0.02)
  43 Vertex_3 = geompy.MakeVertex(-0.365, -0.0575, 0.02)
  44 Vertex_4 = geompy.MakeVertex(-0.365, 0.0575, 0.02)
  45 geomObj_1 = geompy.MakeMarker(0, 0, 0, 1, 0, 0, 0, 1, 0)
  46 Curve_1 = geompy.MakePolyline([Vertex_1, Vertex_4, Vertex_3, Vertex_2], True)
  47 Face_1 = geompy.MakeFaceWires([Curve_1], 1)
  48 [right_side] = geompy.SubShapes(Face_1, [6])
  49 [right_side] = geompy.GetExistingSubObjects(Face_1, False)
  50 path = geompy.MakePrismDXDYDZ(Vertex_3, 0, 0, 0.1)
  51 Translation_1 = geompy.MakeTranslation(Vertex_3, 0, -0.02, 0)
  52 Translation_2 = geompy.MakeTranslation(Vertex_4, 0, 0.02, 0)
  53 Translation_3 = geompy.MakeTranslation(Vertex_4, 0, 0.02, 0)
  54 Line_1 = geompy.MakeLineTwoPnt(Translation_3, Translation_1)
  55 Extrusion_1 = geompy.MakePrismDXDYDZ(Line_1, 0.2, 0, 0)
  56 Partition_1 = geompy.MakePartition([Extrusion_1], [Face_1], [], [], geompy.ShapeType["FACE"], 0, [], 0)
  57 Extrusion_2 = geompy.MakePrismDXDYDZ(Partition_1, 0, 0, 0.1)
  58 [Face_to_enforce] = geompy.SubShapes(Extrusion_2, [30])
  59 geompy.addToStudy( O, 'O' )
  60 geompy.addToStudy( OX, 'OX' )
  61 geompy.addToStudy( OY, 'OY' )
  62 geompy.addToStudy( OZ, 'OZ' )
  63 geompy.addToStudy( Vertex_1, 'Vertex_1' )
  64 geompy.addToStudy( Vertex_2, 'Vertex_2' )
  65 geompy.addToStudy( Vertex_3, 'Vertex_3' )
  66 geompy.addToStudy( Vertex_4, 'Vertex_4' )
  67 geompy.addToStudy( Curve_1, 'Curve_1' )
  68 geompy.addToStudy( Face_1, 'Face_1' )
  69 geompy.addToStudyInFather( Face_1, right_side, 'right_side' )
  70 geompy.addToStudy( path, 'path' )
  71 geompy.addToStudy( Translation_2, 'Translation_2' )
  72 geompy.addToStudy( Translation_1, 'Translation_1' )
  73 geompy.addToStudy( Translation_3, 'Translation_3' )
  74 geompy.addToStudy( Line_1, 'Line_1' )
  75 geompy.addToStudy( Extrusion_1, 'Extrusion_1' )
  76 geompy.addToStudy( Partition_1, 'Partition_1' )
  77 geompy.addToStudy( Extrusion_2, 'Extrusion_2' )
  78 geompy.addToStudyInFather( Extrusion_2, Face_to_enforce, 'Face_to_enforce' )
  79
  80
  81 # Create a group with other faces than Face_to_enforce
  82 all_faces = geompy.SubShapeAll(Extrusion_2, geompy.ShapeType["FACE"])
  83 gr_other_faces = geompy.CreateGroup(Extrusion_2, geompy.ShapeType["FACE"])
  84 geompy.UnionList(gr_other_faces, all_faces)
  85 geompy.DifferenceList(gr_other_faces, [Face_to_enforce])
  86 geompy.addToStudyInFather( Extrusion_2, gr_other_faces, 'other_faces' )
  87
  88 geom_volume = geompy.BasicProperties(Extrusion_2)[2]
  89
  90 ###
  91 ### SMESH component
  92 ###
  93
  94 import  SMESH, SALOMEDS
  95 from salome.smesh import smeshBuilder
  96
  97 smesh = smeshBuilder.New()
  98
  99 Mesh_1 = smesh.Mesh(Face_1,'Mesh_1')
 100
 101 Regular_1D = Mesh_1.Segment()
 102 Number_of_Segments_1 = Regular_1D.NumberOfSegments(10)
 103
 104 Quadrangle_2D = Mesh_1.Quadrangle(algo=smeshBuilder.QUADRANGLE)
 105
 106 right_side_1 = Mesh_1.GroupOnGeom(right_side,'right_side',SMESH.EDGE)
 107 Regular_1D_1 = Mesh_1.Segment(geom=right_side)
 108
 109
 110 Number_of_Segments_2 = Regular_1D_1.NumberOfSegments(nb_segs_right)
 111 Propagation_of_1D_Hyp = Regular_1D_1.Propagation()
 112
 113 isDone = Mesh_1.Compute()
 114
 115 Mesh_path = smesh.Mesh(path,'Mesh_path')
 116
 117 Regular_1D_2 = Mesh_path.Segment()
 118 Local_Length_1 = Regular_1D_2.LocalLength(0.01,None,1e-07)
 119
 120 isDone = Mesh_path.Compute()
 121
 122 ([ right_side_extruded, right_side_top ], error) = Mesh_1.ExtrusionAlongPathObjects( [], [ Mesh_1 ], [ Mesh_1 ], Mesh_path, None, 1, 0, [  ], 0, 0, [ 0, 0, 0 ], 1, [  ], 0 )
 123
 124 Mesh_2 = smesh.Mesh(Extrusion_2,'Mesh_2')
 125
 126 # Enforce elements from Mesh_1
 127 Import_1D2D = Mesh_2.UseExisting2DElements(geom=Face_to_enforce)
 128 Source_Faces_1 = Import_1D2D.SourceFaces([ right_side_extruded ],0,0)
 129
 130 if algo != "GMSH":
 131   # CADSurf is global mesh
 132   MG_CADSurf = Mesh_2.Triangle(algo=smeshBuilder.MG_CADSurf)
 133 else:
 134   # CADSurf is a submesh on other faces
 135   MG_CADSurf = Mesh_2.Triangle(algo=smeshBuilder.MG_CADSurf, geom=gr_other_faces)
 136 MG_CADSurf_Parameters_1 = MG_CADSurf.Parameters()
 137 MG_CADSurf_Parameters_1.SetGeometricMesh( 0 )
 138 MG_CADSurf_Parameters_1.SetPhySize( 0.005 )
 139 MG_CADSurf_Parameters_1.SetMinSize( 0.005 )
 140 MG_CADSurf_Parameters_1.SetMaxSize( 0.005 )
 141 MG_CADSurf_Parameters_1.SetCorrectSurfaceIntersection( False )
 142 MG_CADSurf_Parameters_1.SetElementType( 1 ) # Quadrangle preference
 143 if algo == "MG-Tetra":
 144   Mesh_2.Tetrahedron(algo=smeshBuilder.MG_Tetra)
 145 elif algo == "Netgen":
 146   Mesh_2.Tetrahedron()
 147 elif algo == "GMSH":
 148   GMSH = Mesh_2.Tetrahedron(algo=smeshBuilder.GMSH)
 149   Gmsh_Parameters = GMSH.Parameters()
 150   Gmsh_Parameters.Set2DAlgo( 0 )
 151   Gmsh_Parameters.SetMinSize( 0 )
 152   Gmsh_Parameters.SetMaxSize(  0.005 )
 153   if algo_gmsh == "HXT":
 154     Gmsh_Parameters.Set3DAlgo( 4 )
 155   else:
 156     Gmsh_Parameters.Set3DAlgo( 0 )
 157   Gmsh_Parameters.SetIs2d( 0 )
 158
 159 Face_to_enforce_1 = Mesh_2.GroupOnGeom(Face_to_enforce,'Face_to_enforce',SMESH.FACE)
 160
 161
 162 isDone = Mesh_2.Compute()
 163
 164 Mesh_2.MakeGroup("pyramids", SMESH.VOLUME, CritType=SMESH.FT_ElemGeomType, Threshold=SMESH.Geom_PYRAMID)
 165
 166 if salome.sg.hasDesktop():
 167   salome.sg.updateObjBrowser()
 168
 169 if not isDone:
 170   raise Exception("Error on mesh compute")
 171
 172 mesh_volume = Mesh_2.GetVolume()
 173 assert abs(mesh_volume-geom_volume)/geom_volume < 1e-7, "Wrong mesh volume"
 174
 175 min_aspect_ratio, max_aspect_ratio = Mesh_2.GetMinMax(SMESH.FT_AspectRatio3D)
 176 assert max_aspect_ratio < 200, "Bad aspect ratio 3D: %.1f"%max_aspect_ratio
 177
 178 min_volume, max_volume = Mesh_2.GetMinMax(SMESH.FT_Volume3D)
 179 assert min_volume > 1e-10, "Bad min volume: %s"%min_volume