src/SMESH_SWIG/ex29_refine.py

   1 # CEA/LGLS-LETR 2008, Francis KLOSS (OCC)
   2 # =======================================
   3
   4 # Procedure that take a triangulation and split all triangles in 4 others triangles
   5
   6 import geompy
   7 import smesh
   8
   9 import os
  10
  11 # Values
  12 # ------
  13
  14 # Path for ".med" files
  15 path = "/tmp/ex29_%s_" % os.getenv('USER','unknown')
  16
  17 # Name of the shape and the mesh
  18 name = "Carre"
  19
  20 # Add a node and needed edges
  21 # ---------------------------
  22
  23 def node(m, f, n1, n2, lnv):
  24     x1, y1, z1 = m.GetNodeXYZ(n1)
  25     x2, y2, z2 = m.GetNodeXYZ(n2)
  26
  27     x = (x1 + x2) / 2.0
  28     y = (y1 + y2) / 2.0
  29     z = (z1 + z2) / 2.0
  30
  31     i = m.AddNode(x, y, z)
  32
  33     in1 = m.GetShapeID(n1)
  34     in2 = m.GetShapeID(n2)
  35
  36     if (in1==f) or (in2==f):
  37         m.SetNodeOnFace(i, f, 0, 0)
  38
  39     else:
  40         e1 = m.AddEdge([ n1, i  ])
  41         e2 = m.AddEdge([ i , n2 ])
  42
  43         if n1 in lnv:
  44             e = in2
  45         else:
  46             e = in1
  47
  48         m.SetMeshElementOnShape(e1, e)
  49         m.SetMeshElementOnShape(e2, e)
  50         m.SetNodeOnEdge(i, e, 0)
  51
  52     return i
  53
  54 # Add a triangle and associate to the CAD face
  55 # --------------------------------------------
  56
  57 def triangle(m, f, n1, n2, n3):
  58     i = m.AddFace([ n1, n2, n3 ])
  59     m.SetMeshElementOnShape(i, f)
  60
  61 # Split all triangles in 4 triangles
  62 # ----------------------------------
  63
  64 def SplitTrianglesIn4(m):
  65     # Get all triangles
  66     triangles = m.GetElementsByType(smesh.FACE)
  67
  68     # Remove all edges
  69     m.RemoveElements(m.GetElementsByType(smesh.EDGE))
  70
  71     # Get the list of nodes (ids) associated with the CAD vertices
  72     shape = m.GetShape()
  73     lnv = []
  74     for v in geompy.SubShapeAll(shape, geompy.ShapeType["VERTEX"]):
  75         lnv = lnv + m.GetSubMeshNodesId(v, True)
  76
  77     # Split every triangle
  78     for t in triangles:
  79         noeud_1, noeud_2, noeud_3 = m.GetElemNodes(t)
  80
  81         face = m.GetShapeIDForElem(t)
  82
  83         noeud_12 = node(m, face, noeud_1, noeud_2, lnv)
  84         noeud_23 = node(m, face, noeud_2, noeud_3, lnv)
  85         noeud_13 = node(m, face, noeud_1, noeud_3, lnv)
  86
  87         triangle(m, face, noeud_1 , noeud_12, noeud_13)
  88         triangle(m, face, noeud_2 , noeud_23, noeud_12)
  89         triangle(m, face, noeud_3 , noeud_13, noeud_23)
  90         triangle(m, face, noeud_12, noeud_23, noeud_13)
  91
  92     # Remove all initial triangles
  93     m.RemoveElements(triangles)
  94
  95     # Merge all identical nodes
  96     m.MergeNodes(m.FindCoincidentNodes(0.0001))
  97
  98 # Build a CAD square
  99 # ------------------
 100
 101 x0 = 0.0 ; y0 = 0.0 ; z0 = 0.0
 102 x1 = 1.0 ; y1 = 0.0 ; z1 = 0.0
 103 x2 = 1.0 ; y2 = 1.0 ; z2 = 0.0
 104 x3 = 0.0 ; y3 = 1.0 ; z3 = 0.0
 105
 106 P0 = geompy.MakeVertex(x0, y0, z0)
 107 P1 = geompy.MakeVertex(x1, y1, z1)
 108 P2 = geompy.MakeVertex(x2, y2, z2)
 109 P3 = geompy.MakeVertex(x3, y3, z3)
 110
 111 square = geompy.MakeQuad4Vertices(P0, P1, P2, P3)
 112 geompy.addToStudy(square, name)
 113
 114 # Refine edges and create group of mesh
 115 # -------------------------------------
 116
 117 def refine(m, p1, p2, n, k, name):
 118     s = m.GetShape()
 119
 120     g = geompy.CreateGroup(s, geompy.ShapeType["EDGE"])
 121     e = geompy.GetEdge(s, p1, p2)
 122     i = geompy.GetSubShapeID(s, e)
 123     geompy.AddObject(g, i)
 124     m.Group(g, name)
 125
 126     a = m.Segment(e)
 127     a.NumberOfSegments(n, k)
 128
 129 # Mesh the square
 130 # ---------------
 131
 132 MyMesh = smesh.Mesh(square)
 133
 134 refine(MyMesh, P1, P2,  8,  7, "Droite")
 135 refine(MyMesh, P3, P0,  9, 10, "Gauche")
 136 refine(MyMesh, P0, P1,  7,  9, "Bas"   )
 137 refine(MyMesh, P2, P3, 12, 14, "Haut"  )
 138
 139 algo2D = MyMesh.Triangle()
 140 algo2D.MaxElementArea(0.07)
 141
 142 MyMesh.Compute()
 143
 144 MyMesh.ExportMED(path+"110_triangles.med", 0)
 145
 146 # Disturb the mesh
 147 # ----------------
 148
 149 MyMesh.MoveNode( 37, 0.05    , 0.368967 , 0 )
 150 MyMesh.MoveNode( 38, 0.34    , 0.0762294, 0 )
 151 MyMesh.MoveNode( 40, 0.8     , 0.42     , 0 )
 152 MyMesh.MoveNode( 42, 0.702662, 0.74     , 0 )
 153 MyMesh.MoveNode( 46, 0.4     , 0.374656 , 0 )
 154 MyMesh.MoveNode( 47, 0.13    , 0.63     , 0 )
 155 MyMesh.MoveNode( 49, 0.222187, 0.3      , 0 )
 156 MyMesh.MoveNode( 54, 0.557791, 0.05     , 0 )
 157 MyMesh.MoveNode( 55, 0.7     , 0.2      , 0 )
 158 MyMesh.MoveNode( 56, 0.73    , 0.52     , 0 )
 159 MyMesh.MoveNode( 58, 0.313071, 0.31     , 0 )
 160 MyMesh.MoveNode( 59, 0.8     , 0.56     , 0 )
 161 MyMesh.MoveNode( 62, 0.592703, 0.95     , 0 )
 162 MyMesh.MoveNode( 63, 0.28    , 0.5      , 0 )
 163 MyMesh.MoveNode( 65, 0.49    , 0.93     , 0 )
 164 MyMesh.MoveNode( 68, 0.501038, 0.65     , 0 )
 165 MyMesh.MoveNode( 69, 0.37    , 0.63     , 0 )
 166 MyMesh.MoveNode( 70, 0.597025, 0.52     , 0 )
 167 MyMesh.MoveNode( 72, 0.899   , 0.878589 , 0 )
 168 MyMesh.MoveNode( 73, 0.92    , 0.85     , 0 )
 169 MyMesh.MoveNode( 74, 0.820851, 0.75     , 0 )
 170
 171 NbCells1 = 110
 172 MyMesh.ExportMED(path+"110_triangles_2.med", 0)
 173
 174 # First mesh refining
 175 # -------------------
 176
 177 SplitTrianglesIn4(MyMesh)
 178
 179 NbCells2 = NbCells1*4
 180 print("Mesh with "+str(NbCells2)+" cells computed.")
 181
 182 MyMesh.ExportMED(path+str(NbCells2)+"_triangles.med", 0)
 183
 184 # Second mesh refining
 185 # --------------------
 186
 187 SplitTrianglesIn4(MyMesh)
 188
 189 NbCells3 = NbCells2*4
 190 print("Mesh with "+str(NbCells3)+" cells computed.")
 191
 192 MyMesh.ExportMED(path+str(NbCells3)+"_triangles.med",0)
 193
 194 # Third mesh refining
 195 # -------------------
 196
 197 SplitTrianglesIn4(MyMesh)
 198
 199 NbCells4 = NbCells3*4
 200 print("Mesh with "+str(NbCells4)+" cells computed.")
 201
 202 MyMesh.ExportMED(path+str(NbCells4)+"_triangles.med", 0)
 203
 204 # Update the object browser
 205 # -------------------------
 206
 207 geompy.salome.sg.updateObjBrowser(1)