src/TEST_PY/cas_2013/cas_2013.py

   1 # !/bin/python
   2 # -*- coding: latin-1 -*-
   3 # Hexa : Creation d'hexaedres
   4
   5 import hexablock
   6 import os
   7
   8 #---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
   9
  10 doc  = hexablock.addDocument ("default")
  11 vx   = doc.addVector (1,0,0)
  12 vy   = doc.addVector (0,1,0)
  13 vz   = doc.addVector (0,0,1)
  14
  15 vxy  = doc.addVector (1,1,0)
  16
  17 nbr_files = 0
  18
  19
  20 # ======================================================= save_vtk
  21 def save_vtk () :
  22
  23     global nbr_files
  24     nom = "lecas%d.vtk" % nbr_files
  25     nbr_files += 1
  26     doc.saveVtk (nom)
  27
  28 # ======================================================= carre
  29 def carre (x) :
  30     return x*x
  31
  32 # ======================================================= get_center
  33 def get_center (quad) :
  34     px = 0
  35     py = 0
  36     pz = 0
  37     for nv in range (4) :
  38         vertex = quad.getVertex (nv)
  39         px += vertex.getX() / 4
  40         py += vertex.getY() / 4
  41         pz += vertex.getZ() / 4
  42     return [ px, py, pz ]
  43 # ======================================================= nearest
  44 def nearest (grid, vertex) :
  45     nbre = grid.countVertex()
  46     dmin   = 1e+6
  47     result = None
  48     px = vertex.getX()
  49     py = vertex.getY()
  50     pz = vertex.getZ()
  51     for nro in range (nbre) :
  52         v1 = grid.getVertex (nro)
  53         d2 = carre(px-v1.getX()) + carre(py-v1.getY()) + carre(pz-v1.getZ())
  54         if (d2 < dmin) :
  55            result = v1
  56            dmin   = d2
  57
  58     print  vertex.getName () , px, py, pz, " -> ", result.getName()
  59     return result
  60
  61 # ======================================================= nearest_quad
  62 def nearest_quad (grid, quad) :
  63     dmin   = 1e+16
  64     result = None
  65     [ox, oy, oz]   = get_center (quad)
  66     nbre   = grid.countQuad ()
  67     for nro in range (nbre) :
  68         q1  = grid.getQuad (nro)
  69         if q1 != None :
  70            [px, py, pz] = get_center (q1)
  71            d2 = carre(px-ox) + carre(py-oy) + carre(pz-oz)
  72            if (d2 < dmin) :
  73               result = q1
  74               dmin   = d2
  75
  76     print  quad.getName () , px, py, pz, " -> ", result.getName()
  77     return result
  78
  79 # ======================================================= insert_cylinder
  80 def insert_cylinder (plaque, nx, ny) :
  81
  82     hexa   = plaque.getHexaIJK (nx,   ny,   0)
  83     xmin =  666  ;   ymin = xmin ;   zmin = xmin
  84     xmax = -666  ;   ymax = xmax ;   zmax = xmax
  85
  86     tabv1 = []
  87     for nv in range (8) :
  88         node = hexa.getVertex (nv)
  89         xmin = min (xmin, node.getX())   ;  xmax = max (xmax, node.getX())
  90         ymin = min (ymin, node.getY())   ;  ymax = max (ymax, node.getY())
  91         zmin = min (zmin, node.getZ())   ;  zmax = max (zmax, node.getZ())
  92         tabv1.append (node)
  93
  94     doc.removeHexa (hexa)
  95     save_vtk ()
  96
  97     dx    = (xmax - xmin)/2
  98     dz    = (zmax - zmin)/2
  99     xorig = (xmin + xmax)/2
 100     yorig = (ymin + ymax)/2
 101     zorig = (zmin + zmax)/2 - dz
 102
 103     orig = doc.addVertex (xorig, yorig, zorig)
 104     nr = 1
 105     na = 4
 106     nh = 1
 107     rext  = dx
 108     rint  = rext/2
 109     haut  = 1
 110     angle = 360
 111     pipe  = doc.makePipeUni (orig, vxy,vz, rint,rext,angle,haut, nr,na,nh)
 112
 113     hexablock.what ()
 114
 115     tabquad  = []
 116     tabv0 = []
 117     for nq in range (4) :
 118         quad = pipe.getQuadJK (1, nq, 0)
 119         tabquad.append (quad)
 120
 121     print  " .. tabquad[0] = ", tabquad[0].getName ()
 122     cible = nearest_quad (plaque, tabquad[0])
 123     tabquad[0]. setColor (5)
 124     cible . setColor (5)
 125     save_vtk ()
 126
 127     va1   = tabquad[0].getVertex (0)
 128     va2   = tabquad[0].getVertex (1)
 129     vb1   = cible.nearestVertex  (va1)
 130     vb2   = cible.nearestVertex  (va2)
 131     doc.setLevel (1)
 132     doc.joinQuadsUni (tabquad, cible, va1, vb1, va2, vb2, 1)
 133     hexablock.what ()
 134     save_vtk ()
 135
 136     return
 137     doc.setLevel (1)
 138     for nv in range (8) :
 139         ier = doc.mergeVertices (tabv0[nv], tabv1[nv])
 140         print "ier = ", ier
 141         save_vtk ()
 142
 143
 144 # ======================================================= test_2013
 145 def test_2013 () :
 146
 147     orig = doc.addVertex (0,0,0)
 148
 149     lx  = 3
 150     ly  = lx
 151     lz  = 1
 152     nx  = 3
 153     ny  = nx
 154     nz  = 1
 155
 156     plaque = doc.makeCartesianUni (orig, vx,vy,vz, lx, ly, lz, nx,ny,nz)
 157     save_vtk ()
 158
 159     insert_cylinder (plaque, 1, 1)
 160     return doc
 161
 162 # ================================================================= Begin
 163
 164 doc = test_2013  ()
 165 doc.addLaws (0.1, True)
 166
 167 mesh_hexas = hexablock.mesh (doc)