Back to index

salome-smesh  6.5.0
ex29_refine.py
Go to the documentation of this file.
00001 #  -*- coding: iso-8859-1 -*-
00002 # Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00003 #
00004 # This library is free software; you can redistribute it and/or
00005 # modify it under the terms of the GNU Lesser General Public
00006 # License as published by the Free Software Foundation; either
00007 # version 2.1 of the License.
00008 #
00009 # This library is distributed in the hope that it will be useful,
00010 # but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012 # Lesser General Public License for more details.
00013 #
00014 # You should have received a copy of the GNU Lesser General Public
00015 # License along with this library; if not, write to the Free Software
00016 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00017 #
00018 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00019 #
00020 
00021 # =======================================
00022 # Procedure that take a triangulation and split all triangles in 4 others triangles
00023 #
00024 import geompy
00025 import smesh
00026 
00027 import os
00028 
00029 # Values
00030 # ------
00031 
00032 # Path for ".med" files
00033 path = "/tmp/ex29_%s_" % os.getenv('USER','unknown')
00034 
00035 # Name of the shape and the mesh
00036 name = "Carre"
00037 
00038 # Add a node and needed edges
00039 # ---------------------------
00040 
00041 def node(m, f, n1, n2, lnv):
00042     x1, y1, z1 = m.GetNodeXYZ(n1)
00043     x2, y2, z2 = m.GetNodeXYZ(n2)
00044 
00045     x = (x1 + x2) / 2.0
00046     y = (y1 + y2) / 2.0
00047     z = (z1 + z2) / 2.0
00048 
00049     i = m.AddNode(x, y, z)
00050 
00051     in1 = m.GetShapeID(n1)
00052     in2 = m.GetShapeID(n2)
00053 
00054     if (in1==f) or (in2==f):
00055         m.SetNodeOnFace(i, f, 0, 0)
00056 
00057     else:
00058         e1 = m.AddEdge([ n1, i  ])
00059         e2 = m.AddEdge([ i , n2 ])
00060 
00061         if n1 in lnv:
00062             e = in2
00063         else:
00064             e = in1
00065 
00066         m.SetMeshElementOnShape(e1, e)
00067         m.SetMeshElementOnShape(e2, e)
00068         m.SetNodeOnEdge(i, e, 0)
00069 
00070     return i
00071 
00072 # Add a triangle and associate to the CAD face
00073 # --------------------------------------------
00074 
00075 def triangle(m, f, n1, n2, n3):
00076     i = m.AddFace([ n1, n2, n3 ])
00077     m.SetMeshElementOnShape(i, f)
00078 
00079 # Split all triangles in 4 triangles
00080 # ----------------------------------
00081 
00082 def SplitTrianglesIn4(m):
00083     # Get all triangles
00084     triangles = m.GetElementsByType(smesh.FACE)
00085 
00086     # Remove all edges
00087     m.RemoveElements(m.GetElementsByType(smesh.EDGE))
00088 
00089     # Get the list of nodes (ids) associated with the CAD vertices
00090     shape = m.GetShape()
00091     lnv = []
00092     for v in geompy.SubShapeAll(shape, geompy.ShapeType["VERTEX"]):
00093         lnv = lnv + m.GetSubMeshNodesId(v, True)
00094 
00095     # Split every triangle
00096     for t in triangles:
00097         noeud_1, noeud_2, noeud_3 = m.GetElemNodes(t)
00098 
00099         face = m.GetShapeIDForElem(t)
00100 
00101         noeud_12 = node(m, face, noeud_1, noeud_2, lnv)
00102         noeud_23 = node(m, face, noeud_2, noeud_3, lnv)
00103         noeud_13 = node(m, face, noeud_1, noeud_3, lnv)
00104 
00105         triangle(m, face, noeud_1 , noeud_12, noeud_13)
00106         triangle(m, face, noeud_2 , noeud_23, noeud_12)
00107         triangle(m, face, noeud_3 , noeud_13, noeud_23)
00108         triangle(m, face, noeud_12, noeud_23, noeud_13)
00109 
00110     # Remove all initial triangles
00111     m.RemoveElements(triangles)
00112 
00113     # Merge all identical nodes
00114     m.MergeNodes(m.FindCoincidentNodes(0.0001))
00115 
00116 # Build a CAD square
00117 # ------------------
00118 
00119 x0 = 0.0 ; y0 = 0.0 ; z0 = 0.0
00120 x1 = 1.0 ; y1 = 0.0 ; z1 = 0.0
00121 x2 = 1.0 ; y2 = 1.0 ; z2 = 0.0
00122 x3 = 0.0 ; y3 = 1.0 ; z3 = 0.0
00123 
00124 P0 = geompy.MakeVertex(x0, y0, z0)
00125 P1 = geompy.MakeVertex(x1, y1, z1)
00126 P2 = geompy.MakeVertex(x2, y2, z2)
00127 P3 = geompy.MakeVertex(x3, y3, z3)
00128 
00129 square = geompy.MakeQuad4Vertices(P0, P1, P2, P3)
00130 geompy.addToStudy(square, name)
00131 
00132 # Refine edges and create group of mesh
00133 # -------------------------------------
00134 
00135 def refine(m, p1, p2, n, k, name):
00136     s = m.GetShape()
00137 
00138     g = geompy.CreateGroup(s, geompy.ShapeType["EDGE"])
00139     e = geompy.GetEdge(s, p1, p2)
00140     i = geompy.GetSubShapeID(s, e)
00141     geompy.AddObject(g, i)
00142     m.Group(g, name)
00143 
00144     a = m.Segment(e)
00145     a.NumberOfSegments(n, k)
00146 
00147 # Mesh the square
00148 # ---------------
00149 
00150 MyMesh = smesh.Mesh(square)
00151 
00152 refine(MyMesh, P1, P2,  8,  7, "Droite")
00153 refine(MyMesh, P3, P0,  9, 10, "Gauche")
00154 refine(MyMesh, P0, P1,  7,  9, "Bas"   )
00155 refine(MyMesh, P2, P3, 12, 14, "Haut"  )
00156 
00157 algo2D = MyMesh.Triangle()
00158 algo2D.MaxElementArea(0.07)
00159 
00160 MyMesh.Compute()
00161 
00162 MyMesh.ExportMED(path+"110_triangles.med", 0)
00163 
00164 # Disturb the mesh
00165 # ----------------
00166 
00167 MyMesh.MoveNode( 37, 0.05    , 0.368967 , 0 )
00168 MyMesh.MoveNode( 38, 0.34    , 0.0762294, 0 )
00169 MyMesh.MoveNode( 40, 0.8     , 0.42     , 0 )
00170 MyMesh.MoveNode( 42, 0.702662, 0.74     , 0 )
00171 MyMesh.MoveNode( 46, 0.4     , 0.374656 , 0 )
00172 MyMesh.MoveNode( 47, 0.13    , 0.63     , 0 )
00173 MyMesh.MoveNode( 49, 0.222187, 0.3      , 0 )
00174 MyMesh.MoveNode( 54, 0.557791, 0.05     , 0 )
00175 MyMesh.MoveNode( 55, 0.7     , 0.2      , 0 )
00176 MyMesh.MoveNode( 56, 0.73    , 0.52     , 0 )
00177 MyMesh.MoveNode( 58, 0.313071, 0.31     , 0 )
00178 MyMesh.MoveNode( 59, 0.8     , 0.56     , 0 )
00179 MyMesh.MoveNode( 62, 0.592703, 0.95     , 0 )
00180 MyMesh.MoveNode( 63, 0.28    , 0.5      , 0 )
00181 MyMesh.MoveNode( 65, 0.49    , 0.93     , 0 )
00182 MyMesh.MoveNode( 68, 0.501038, 0.65     , 0 )
00183 MyMesh.MoveNode( 69, 0.37    , 0.63     , 0 )
00184 MyMesh.MoveNode( 70, 0.597025, 0.52     , 0 )
00185 MyMesh.MoveNode( 72, 0.899   , 0.878589 , 0 )
00186 MyMesh.MoveNode( 73, 0.92    , 0.85     , 0 )
00187 MyMesh.MoveNode( 74, 0.820851, 0.75     , 0 )
00188 
00189 NbCells1 = 110
00190 MyMesh.ExportMED(path+"110_triangles_2.med", 0)
00191 
00192 # First mesh refining
00193 # -------------------
00194 
00195 SplitTrianglesIn4(MyMesh)
00196 
00197 NbCells2 = NbCells1*4
00198 print("Mesh with "+str(NbCells2)+" cells computed.")
00199 
00200 MyMesh.ExportMED(path+str(NbCells2)+"_triangles.med", 0)
00201 
00202 # Second mesh refining
00203 # --------------------
00204 
00205 SplitTrianglesIn4(MyMesh)
00206 
00207 NbCells3 = NbCells2*4
00208 print("Mesh with "+str(NbCells3)+" cells computed.")
00209 
00210 MyMesh.ExportMED(path+str(NbCells3)+"_triangles.med",0)
00211 
00212 # Third mesh refining
00213 # -------------------
00214 
00215 SplitTrianglesIn4(MyMesh)
00216 
00217 NbCells4 = NbCells3*4
00218 print("Mesh with "+str(NbCells4)+" cells computed.")
00219 
00220 MyMesh.ExportMED(path+str(NbCells4)+"_triangles.med", 0)
00221 
00222 # Update the object browser
00223 # -------------------------
00224 
00225 geompy.salome.sg.updateObjBrowser(1)