src/TEST_PY/cas_2013/cas_2013.py

   1 # !/bin/python
   2 # -*- coding: latin-1 -*-
   3 # Copyright (C) 2009-2023  CEA, EDF
   4 #
   5 # This library is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU Lesser General Public
   7 # License as published by the Free Software Foundation; either
   8 # version 2.1 of the License, or (at your option) any later version.
   9 #
  10 # This library is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 # Lesser General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with this library; if not, write to the Free Software
  17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21
  22 # Hexa : Creation d'hexaedres
  23
  24 import hexablock
  25 import os
  26
  27 #---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
  28
  29 doc  = hexablock.addDocument ("default")
  30 vx   = doc.addVector (1,0,0)
  31 vy   = doc.addVector (0,1,0)
  32 vz   = doc.addVector (0,0,1)
  33
  34 vxy  = doc.addVector (1,1,0)
  35
  36 nbr_files = 0
  37
  38
  39 # ======================================================= save_vtk
  40 def save_vtk () :
  41
  42     global nbr_files
  43     nom = "lecas%d.vtk" % nbr_files
  44     nbr_files += 1
  45     doc.saveVtk (nom)
  46
  47 # ======================================================= carre
  48 def carre (x) :
  49     return x*x
  50
  51 # ======================================================= get_center
  52 def get_center (quad) :
  53     px = 0
  54     py = 0
  55     pz = 0
  56     for nv in range (4) :
  57         vertex = quad.getVertex (nv)
  58         px += vertex.getX() / 4
  59         py += vertex.getY() / 4
  60         pz += vertex.getZ() / 4
  61     return [ px, py, pz ]
  62 # ======================================================= nearest
  63 def nearest (grid, vertex) :
  64     nbre = grid.countVertex()
  65     dmin   = 1e+6
  66     result = None
  67     px = vertex.getX()
  68     py = vertex.getY()
  69     pz = vertex.getZ()
  70     for nro in range (nbre) :
  71         v1 = grid.getVertex (nro)
  72         d2 = carre(px-v1.getX()) + carre(py-v1.getY()) + carre(pz-v1.getZ())
  73         if (d2 < dmin) :
  74            result = v1
  75            dmin   = d2
  76
  77     print(vertex.getName () , px, py, pz, " -> ", result.getName())
  78     return result
  79
  80 # ======================================================= nearest_quad
  81 def nearest_quad (grid, quad) :
  82     dmin   = 1e+16
  83     result = None
  84     [ox, oy, oz]   = get_center (quad)
  85     nbre   = grid.countQuad ()
  86     for nro in range (nbre) :
  87         q1  = grid.getQuad (nro)
  88         if q1 != None :
  89            [px, py, pz] = get_center (q1)
  90            d2 = carre(px-ox) + carre(py-oy) + carre(pz-oz)
  91            if (d2 < dmin) :
  92               result = q1
  93               dmin   = d2
  94
  95     print(quad.getName () , px, py, pz, " -> ", result.getName())
  96     return result
  97
  98 # ======================================================= insert_cylinder
  99 def insert_cylinder (plaque, nx, ny) :
 100
 101     hexa   = plaque.getHexaIJK (nx,   ny,   0)
 102     xmin =  666  ;   ymin = xmin ;   zmin = xmin
 103     xmax = -666  ;   ymax = xmax ;   zmax = xmax
 104
 105     tabv1 = []
 106     for nv in range (8) :
 107         node = hexa.getVertex (nv)
 108         xmin = min (xmin, node.getX())   ;  xmax = max (xmax, node.getX())
 109         ymin = min (ymin, node.getY())   ;  ymax = max (ymax, node.getY())
 110         zmin = min (zmin, node.getZ())   ;  zmax = max (zmax, node.getZ())
 111         tabv1.append (node)
 112
 113     doc.removeHexa (hexa)
 114     save_vtk ()
 115
 116     dx    = (xmax - xmin)/2
 117     dz    = (zmax - zmin)/2
 118     xorig = (xmin + xmax)/2
 119     yorig = (ymin + ymax)/2
 120     zorig = (zmin + zmax)/2 - dz
 121
 122     orig = doc.addVertex (xorig, yorig, zorig)
 123     nr = 1
 124     na = 4
 125     nh = 1
 126     rext  = dx
 127     rint  = rext/2
 128     haut  = 1
 129     angle = 360
 130     pipe  = doc.makePipeUni (orig, vxy,vz, rint,rext,angle,haut, nr,na,nh)
 131
 132     hexablock.what ()
 133
 134     tabquad  = []
 135     tabv0 = []
 136     for nq in range (4) :
 137         quad = pipe.getQuadJK (1, nq, 0)
 138         tabquad.append (quad)
 139
 140     print(" .. tabquad[0] = ", tabquad[0].getName ())
 141     cible = nearest_quad (plaque, tabquad[0])
 142     tabquad[0]. setColor (5)
 143     cible . setColor (5)
 144     save_vtk ()
 145
 146     va1   = tabquad[0].getVertex (0)
 147     va2   = tabquad[0].getVertex (1)
 148     vb1   = cible.nearestVertex  (va1)
 149     vb2   = cible.nearestVertex  (va2)
 150     doc.setLevel (1)
 151     doc.joinQuadsUni (tabquad, cible, va1, vb1, va2, vb2, 1)
 152     hexablock.what ()
 153     save_vtk ()
 154
 155     return
 156     doc.setLevel (1)
 157     for nv in range (8) :
 158         ier = doc.mergeVertices (tabv0[nv], tabv1[nv])
 159         print("ier = ", ier)
 160         save_vtk ()
 161
 162
 163 # ======================================================= test_2013
 164 def test_2013 () :
 165
 166     orig = doc.addVertex (0,0,0)
 167
 168     lx  = 3
 169     ly  = lx
 170     lz  = 1
 171     nx  = 3
 172     ny  = nx
 173     nz  = 1
 174
 175     plaque = doc.makeCartesianUni (orig, vx,vy,vz, lx, ly, lz, nx,ny,nz)
 176     save_vtk ()
 177
 178     insert_cylinder (plaque, 1, 1)
 179     return doc
 180
 181 # ================================================================= Begin
 182
 183 doc = test_2013  ()
 184 doc.addLaws (0.1, True)
 185
 186 mesh_hexas = hexablock.mesh (doc)