src/TEST_PY/recettes/tuyauterie.py

   1 # -*- coding: latin-1 -*-
   2 # Copyright (C) 2009-2024  CEA, EDF
   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 # Francis KLOSS : 2012 : CEA-Saclay, DEN, DM2S, SFME, LGLS, F-91191 Gif-sur-Yvette, France
  22 # ========================================================================================
  23
  24 import math
  25 import hexablock
  26 ### import geompy
  27 geompy = hexablock.geompy
  28
  29 # Charger la géométrie
  30 # ====================
  31
  32 nom = "tuyauterie"
  33
  34 geometrie = geompy.ImportBREP(nom+".brep")
  35
  36 # Sélectionner des sous-parties de la géométrie
  37 # ---------------------------------------------
  38
  39 sommets = geompy.SubShapeAllSortedCentres(geometrie, geompy.ShapeType["VERTEX"])
  40 aretes  = geompy.SubShapeAllSortedCentres(geometrie, geompy.ShapeType["EDGE"  ])
  41 faces   = geompy.SubShapeAllSortedCentres(geometrie, geompy.ShapeType["FACE"  ])
  42
  43 petit_arri_v = sommets[0]
  44 petit_arri_x, petit_arri_y, petit_arri_z = geompy.PointCoordinates(petit_arri_v)
  45 petit_avan_x, petit_avan_y, petit_avan_z = geompy.PointCoordinates(sommets[1])
  46
  47 moyen_arri_v = sommets[3]
  48 moyen_arri_x, moyen_arri_y, moyen_arri_z = geompy.PointCoordinates(moyen_arri_v)
  49 moyen_avan_x, moyen_avan_y, moyen_avan_z = geompy.PointCoordinates(sommets[4])
  50
  51 petit0, petit_x, petit_y, petit_z,  petit1, petit2, petit3,  petit_rayon = geompy.KindOfShape(aretes[0])
  52
  53 moyen0, moyen_x, moyen_y, moyen_z,  moyen1, moyen2, moyen3,  moyen_rayon = geompy.KindOfShape(aretes[4])
  54
  55 grand0, grand_xb, grand_yb, grand_zb,  grand1, grand2, grand3,  grand_rayon = geompy.KindOfShape(aretes[7])
  56 grand4, grand_xh, grand_yh, grand_zh,  grand5, grand6, grand7,  grand8      = geompy.KindOfShape(aretes[8])
  57
  58 grand_generat = aretes[11]
  59
  60 arc0,  arc_x, arc_y, arc_z,  arc_dx, arc_dy, arc_dz,  arc_rayon, arc_x1, arc_y1, arc_z1,  arc_x2, arc_y2, arc_z2 = geompy.KindOfShape(aretes[12])
  61
  62 face_cylindre = faces[5]
  63
  64 # Construire des géométries intermédiaires
  65 # ----------------------------------------
  66
  67 arc_sc = geompy.MakeVertex(arc_x , arc_y , arc_z )
  68 arc_s1 = geompy.MakeVertex(arc_x1, arc_y1, arc_z1)
  69 arc_s2 = geompy.MakeVertex(arc_x2, arc_y2, arc_z2)
  70
  71 arc_a1 = geompy.MakeEdge(arc_sc, arc_s1)
  72 arc_a2 = geompy.MakeEdge(arc_sc, arc_s2)
  73
  74 # Construire le modèle de bloc
  75 # ============================
  76
  77 doc = hexablock.addDocument(nom)
  78
  79 # Construire le grand cylindre en 2 parties
  80 # -----------------------------------------
  81
  82 grand_base = doc.addVertex(grand_xb, grand_yb, grand_zb)
  83 grand_oppo = doc.addVertex(grand_xh, grand_yh, grand_zh)
  84
  85 grand_dir_a = doc.addVectorVertices(grand_base, grand_oppo)
  86 grand_dir_b = doc.addVectorVertices(grand_oppo, grand_base)
  87
  88 grand_hauteur_a = grand_dir_a.getNorm()*0.20
  89 grand_hauteur_b = grand_dir_b.getNorm()*0.60
  90
  91 ####   grand_cylindre_a = doc.addCylinder(grand_base, grand_dir_a, grand_rayon, grand_hauteur_a)
  92 ####   grand_cylindre_b = doc.addCylinder(grand_oppo, grand_dir_b, grand_rayon, grand_hauteur_b)
  93
  94 # Construire le moyen cylindre
  95 # ----------------------------
  96
  97 moyen_base = doc.addVertex(moyen_x, moyen_y, moyen_z)
  98
  99 moyen_arri = doc.addVertex(moyen_arri_x, moyen_arri_y, moyen_arri_z)
 100 moyen_avan = doc.addVertex(moyen_avan_x, moyen_avan_y, moyen_avan_z)
 101
 102 moyen_dir = doc.addVectorVertices(moyen_arri, moyen_avan)
 103
 104 moyen_hauteur = geompy.MinDistance(moyen_arri_v, grand_generat) - grand_rayon
 105
 106
 107 # Construire le petit cylindre
 108 # ----------------------------
 109
 110 petit_base = doc.addVertex(petit_x, petit_y, petit_z)
 111
 112 petit_arri = doc.addVertex(petit_arri_x, petit_arri_y, petit_arri_z)
 113 petit_avan = doc.addVertex(petit_avan_x, petit_avan_y, petit_avan_z)
 114
 115 petit_dir = doc.addVectorVertices(petit_arri, petit_avan)
 116
 117 petit_hauteur = geompy.MinDistance(petit_arri_v, grand_generat)
 118
 119
 120 # Construire les 2 cylindres qui s'intersectent
 121 # ---------------------------------------------
 122
 123 moyen_rayon  = moyen_rayon*0.7
 124 cylindres_gm = doc.makeCylinders (grand_base, grand_dir_a, grand_rayon, grand_hauteur_a,
 125                                   moyen_base, moyen_dir, moyen_rayon, moyen_hauteur)
 126 cylindres_gp = doc.makeCylinders (grand_oppo, grand_dir_b, grand_rayon, grand_hauteur_b,
 127                                   petit_base, petit_dir, petit_rayon, petit_hauteur)
 128
 129 doc.saveVtk ("tuyauterie0.vtk")
 130
 131 print("grand_rayon = ", grand_rayon)
 132 print("moyen_rayon = ", moyen_rayon)
 133 print("petit_rayon = ", petit_rayon)
 134
 135
 136 # Joindre les 2 croix
 137 # -------------------
 138
 139 CYL_BIG  = 1
 140 CYL_KMAX = 3
 141 CYL_JMAX = 4
 142
 143 qstart = cylindres_gm.getQuadIJ (CYL_BIG, 0,0,CYL_KMAX)
 144
 145 qdest  = cylindres_gp.getQuadIJ (CYL_BIG, 0,0,CYL_KMAX)
 146
 147 va1 = qstart.getVertex (0)
 148 va2 = qstart.getVertex (1)
 149 vb1 = qdest .nearestVertex (va1)
 150 vb2 = qdest .nearestVertex (va2)
 151
 152 gm_quads = [qstart]
 153
 154 for nj in range (CYL_JMAX) :
 155     quad = cylindres_gm.getQuadIJ (CYL_BIG, 1,nj, CYL_KMAX)
 156     gm_quads.append (quad)
 157
 158 prisme = doc.joinQuadsUni (gm_quads, qdest, va1,vb1,va2,vb2, 1)
 159
 160
 161 # Ajouter le coude au grand cylindre
 162 # ----------------------------------
 163
 164 coude_quads = [ cylindres_gp.getQuadIJ(CYL_BIG, 0, 0, 0) ]
 165
 166 for nj in range (CYL_JMAX) :
 167         quad = cylindres_gp.getQuadIJ(CYL_BIG, 1, nj, 0)
 168         coude_quads.append(quad)
 169
 170 coude_centre = doc.addVertex(arc_x , arc_y , arc_z )
 171 coude_dir    = doc.addVector(arc_dx, arc_dy, arc_dz)
 172 coude_angle  = geompy.GetAngle(arc_a1, arc_a2)
 173
 174 coude = doc.revolutionQuads (coude_quads, coude_centre, coude_dir, [coude_angle])
 175 ##coude = doc.revolutionQuadsUni (coude_quads, coude_centre, coude_dir, coude_angle, 8)
 176
 177 doc.saveVtk ("tuyauterie1.vtk")
 178 doc.addLaws (0.02, True)
 179 blocs = hexablock.mesh(doc)
 180 muv, mue, muq, muh = hexablock.dump(doc, blocs)
 181
 182 arret_procedure ()
 183
 184 # Associer la géométrie au modèle de bloc
 185 # =======================================
 186
 187 doc.addShape(geometrie, nom)
 188
 189 # Associer les arêtes du modèle issues du joinQuads
 190 # -------------------------------------------------
 191
 192 for i in range(hexablock.CV_MAXI_EXT):
 193     asso_h  = prisme.getHexa(i)
 194     asso_a  = [ hexablock.E_CE, hexablock.E_CF, hexablock.E_CE, hexablock.E_CF, hexablock.E_CE, hexablock.E_CF, hexablock.E_CE, hexablock.E_CF ][i]
 195     asso_e  = asso_h.getEdge(asso_a)
 196     asso_s1 = asso_e.getVertex(0)
 197     asso_s2 = asso_e.getVertex(1)
 198     asso_v1 = asso_s1.getAssociation()
 199     asso_v2 = asso_s2.getAssociation()
 200     x1, y1, z1 = geompy.PointCoordinates(asso_v1)
 201     x2, y2, z2 = geompy.PointCoordinates(asso_v2)
 202     asso_pts = [ asso_v1 ]
 203     n = 9
 204     pas = math.sqrt( (x2-x1)**2 + (y2-y1)**2 + (z2-z1)**2 ) / n
 205     for i in range(1, n):
 206         x = x1 + (x2 - x1) * i * pas
 207         y = y1 + (y2 - y1) * i * pas
 208         z = z1 + (z2 - z1) * i * pas
 209         asso_v = geompy.MakeVertexOnSurfaceByCoord(face_cylindre, x, y, z)
 210         asso_pts.append(asso_v)
 211     asso_pts.append(asso_v2)
 212     asso_curve = geompy.MakeInterpol(asso_pts)
 213     asso_e.addAssociation(asso_curve, 0, 1)
 214
 215 # Associer les quadrangles du modèle issues du joinQuads
 216 # ------------------------------------------------------
 217
 218 for i in range(hexablock.CV_MAXI_EXT):
 219     asso_h = prisme.getHexa(i)
 220     asso_f = [ hexablock.Q_C, hexablock.Q_F, hexablock.Q_E, hexablock.Q_F, hexablock.Q_E, hexablock.Q_F, hexablock.Q_E, hexablock.Q_F ][i]
 221     asso_q = asso_h.getQuad(asso_f)
 222     asso_q.addAssociation(face_cylindre)
 223
 224 # Mailler le modèle de bloc
 225 # =========================
 226
 227 # Définir 4 groupes de faces
 228 # --------------------------
 229
 230 groupe_d_grand = doc.addQuadGroup("grand:disque")
 231 groupe_d_moyen = doc.addQuadGroup("moyen:disque")
 232 groupe_d_petit = doc.addQuadGroup("petit:disque")
 233 groupe_d_coude = doc.addQuadGroup("coude:disque")
 234
 235 for i in range(2):
 236     for j in range( [hexablock.CV_MAXI_INT, hexablock.CV_MAXI_EXT][i] ):
 237         quad = cylindres_gm.getQuadIJ(hexablock.CYL_BIG  , i, j, 0)
 238         groupe_d_grand.addElement(quad)
 239
 240         quad = cylindres_gm.getQuadIJ(hexablock.CYL_SMALL, i, j, 0)
 241         groupe_d_moyen.addElement(quad)
 242
 243         quad = cylindres_gp.getQuadIJ(hexablock.CYL_SMALL, i, j, 0)
 244         groupe_d_petit.addElement(quad)
 245
 246 for i in range( hexablock.CV_MAXI_INT + hexablock.CV_MAXI_EXT ):
 247     h = coude.getHexa(i)
 248     quad = h.getQuad(hexablock.Q_B)
 249     groupe_d_coude.addElement(quad)
 250
 251 # Définir 3 groupes de volumes
 252 # ----------------------------
 253
 254 groupe_grand = doc.addHexaGroup("grand")
 255 groupe_moyen = doc.addHexaGroup("moyen")
 256 groupe_petit = doc.addHexaGroup("petit")
 257
 258 for i in range( doc.countUsedHexa() ):
 259     h = doc.getUsedHexa(i)
 260     groupe_grand.addElement(h)
 261
 262 for i in range(2):
 263     for j in range( [hexablock.CV_MAXI_INT, hexablock.CV_MAXI_EXT][i] ):
 264         h = cylindres_gm.getHexaIJK(hexablock.CYL_SMALL, i, j, 0)
 265         groupe_moyen.addElement(h)
 266         groupe_grand.removeElement(h)
 267
 268         h = cylindres_gp.getHexaIJK(hexablock.CYL_SMALL, i, j, 0)
 269         groupe_petit.addElement(h)
 270         groupe_grand.removeElement(h)
 271
 272 # Mailler le modèle de bloc avec ses associations
 273 # -----------------------------------------------
 274
 275 hexablock.addLaws(doc, 0.02, True)
 276
 277 blocs = hexablock.mesh(doc)
 278
 279 muv, mue, muq, muh = hexablock.dump(doc, blocs)