test/ex29_refine.py

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