src/Tools/blocFissure/gmu/fissureCoude.py

   1 # -*- coding: utf-8 -*-
   2 # Copyright (C) 2014-2020  EDF R&D
   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 import os
  22
  23 from .geomsmesh import geompy, smesh
  24 from .geomsmesh import geomPublish
  25 from .geomsmesh import geomPublishInFather
  26 from . import initLog
  27
  28 import math
  29 import GEOM
  30 import SALOMEDS
  31 import SMESH
  32 import logging
  33
  34 from .fissureGenerique import fissureGenerique
  35
  36 from .triedreBase import triedreBase
  37 from .genereMeshCalculZoneDefaut import genereMeshCalculZoneDefaut
  38 from .creeZoneDefautDansObjetSain import creeZoneDefautDansObjetSain
  39 from .construitFissureGenerale import construitFissureGenerale
  40 from .sortEdges import sortEdges
  41
  42 O, OX, OY, OZ = triedreBase()
  43
  44 class fissureCoude(fissureGenerique):
  45   """
  46   problème de fissure du Coude : version de base
  47   maillage hexa
  48   """
  49
  50   nomProbleme = "fissureCoude"
  51
  52   # ---------------------------------------------------------------------------
  53   def setParamGeometrieSaine(self):
  54     """
  55     Paramètres géométriques du tuyau coudé sain:
  56     angleCoude
  57     r_cintr
  58     l_tube_p1
  59     l_tube_p2
  60     epais
  61     de
  62     """
  63     self.geomParams = dict(angleCoude = 60,
  64                            r_cintr    = 1200,
  65                            l_tube_p1  = 1600,
  66                            l_tube_p2  = 1200,
  67                            epais      = 40,
  68                            de         = 760)
  69
  70   # ---------------------------------------------------------------------------
  71   def genereGeometrieSaine(self, geomParams):
  72     """a écrire"""
  73     logging.info("genereGeometrieSaine %s", self.nomCas)
  74
  75     angleCoude = geomParams['angleCoude']
  76     r_cintr    = geomParams['r_cintr']
  77     l_tube_p1  = geomParams['l_tube_p1']
  78     l_tube_p2  = geomParams['l_tube_p2']
  79     epais      = geomParams['epais']
  80     de         = geomParams['de']
  81
  82     centre = geompy.MakeVertex(0, 0, -l_tube_p1)
  83     diskext = geompy.MakeDiskPntVecR(centre, OZ, de/2.)
  84     diskint = geompy.MakeDiskPntVecR(centre, OZ, de/2. -epais)
  85     couronne = geompy.MakeCut(diskext, diskint)
  86     tube_1 = geompy.MakePrismVecH(couronne, OZ, l_tube_p1)
  87     axe = geompy.MakeTranslation(OY, -r_cintr, 0, -l_tube_p1)
  88     coude = geompy.MakeRevolution(couronne, axe, angleCoude*math.pi/180.0)
  89     Rotation_1 = geompy.MakeRotation(couronne, axe, angleCoude*math.pi/180.0)
  90     Rotation_2 = geompy.MakeRotation(OZ, OY, angleCoude*math.pi/180.0)
  91     tube_2 = geompy.MakePrismVecH(Rotation_1, Rotation_2, -l_tube_p2)
  92     plan_y = geompy.MakePlaneLCS(None, 100000, 3)
  93     geomPublish(initLog.debug,  plan_y, "plan_y" )
  94     geomPublish(initLog.debug,  tube_1, "tube_1" )
  95     geomPublish(initLog.debug,  coude, "coude" )
  96     geomPublish(initLog.debug,  tube_2, "tube_2" )
  97
  98     P1 = O
  99     geompy.addToStudy(P1, "P1" )
 100     op2 = geompy.MakeVertex(0, 0, -l_tube_p1)
 101     P2 = geompy.MakeRotation(op2, axe, angleCoude*math.pi/180.0)
 102     P2 = geompy.MakeTranslationVectorDistance(P2, Rotation_2, -l_tube_p2)
 103     geompy.addToStudy(P2, "P2" )
 104
 105     # --- tube coude sain
 106
 107     geometrieSaine = geompy.MakePartition([tube_1, coude, tube_2, P1, P2], [plan_y], [], [], geompy.ShapeType["SOLID"], 0, [], 1)
 108     geomPublish(initLog.debug,  geometrieSaine, self.nomCas )
 109     [P1, P2] = geompy.RestoreGivenSubShapes(geometrieSaine, [P1, P2], GEOM.FSM_GetInPlaceByHistory, False, True)
 110
 111     xmin = -de -r_cintr -l_tube_p2
 112     zmin = -l_tube_p1 -r_cintr -l_tube_p2 -de
 113     ymax = de +100.
 114     boxypos = geompy.MakeBox(xmin, 0, zmin, ymax, ymax, 100, "boxypos")
 115     boxyneg = geompy.MakeBox(xmin, 0, zmin, ymax, -ymax, 100, "boxyneg")
 116     edgesypos = geompy.GetShapesOnShape(boxypos, geometrieSaine, geompy.ShapeType["EDGE"], GEOM.ST_IN)
 117     edgesyneg = geompy.GetShapesOnShape(boxyneg, geometrieSaine, geompy.ShapeType["EDGE"], GEOM.ST_IN)
 118     circ_g = geompy.CreateGroup(geometrieSaine, geompy.ShapeType["EDGE"])
 119     geompy.UnionList(circ_g, edgesyneg)
 120     circ_d = geompy.CreateGroup(geometrieSaine, geompy.ShapeType["EDGE"])
 121     geompy.UnionList(circ_d, edgesypos)
 122     edgesy0pos = geompy.GetShapesOnShape(boxypos, geometrieSaine, geompy.ShapeType["EDGE"], GEOM.ST_ONIN)
 123     grpedpos = geompy.CreateGroup(geometrieSaine, geompy.ShapeType["EDGE"])
 124     geompy.UnionList(grpedpos, edgesy0pos)
 125     grpedy0 = geompy.CutGroups(grpedpos, circ_d, "edges_y0")
 126     boxtub1 = geompy.MakeBox(-de/2.0 -1, -1, -l_tube_p1, de, de, 0, "boxtub1")
 127     edgestub1 = geompy.GetShapesOnShape(boxtub1, geometrieSaine, geompy.ShapeType["EDGE"], GEOM.ST_IN)
 128     grped = geompy.CreateGroup(geometrieSaine, geompy.ShapeType["EDGE"])
 129     geompy.UnionList(grped, edgestub1)
 130     long_p1 = geompy.IntersectGroups(grped, grpedy0)
 131     boxtub  = geompy.MakeBox(-de/2.0 -1, -1, -l_tube_p1 -l_tube_p2, de, de, -l_tube_p1)
 132     boxtub2 = geompy.MakeRotation(boxtub, axe, angleCoude*math.pi/180.0, "boxttub2")
 133     edgestub2 = geompy.GetShapesOnShape(boxtub2, geometrieSaine, geompy.ShapeType["EDGE"], GEOM.ST_IN)
 134     grped = geompy.CreateGroup(geometrieSaine, geompy.ShapeType["EDGE"])
 135     geompy.UnionList(grped, edgestub2)
 136     long_p2 = geompy.IntersectGroups(grped, grpedy0)
 137     boxtub1t = geompy.MakeTranslationVectorDistance(boxtub1, OZ, -l_tube_p1)
 138     facer = geompy.GetShapesOnShape(boxtub1t, boxtub1, geompy.ShapeType["FACE"], GEOM.ST_ONIN, "facer")
 139     boxcoud = geompy.MakeRevolution(facer[0], axe, angleCoude*math.pi/180.0, "boxcoud")
 140     edgescoud = geompy.GetShapesOnShape(boxcoud, geometrieSaine, geompy.ShapeType["EDGE"], GEOM.ST_IN)
 141     grped = geompy.CreateGroup(geometrieSaine, geompy.ShapeType["EDGE"])
 142     geompy.UnionList(grped, edgescoud)
 143     long_coude = geompy.IntersectGroups(grped, grpedy0)
 144     grped = geompy.CutGroups(grpedy0, long_p1)
 145     grped = geompy.CutGroups(grped, long_p2)
 146     ep = geompy.CutGroups(grped, long_coude)
 147     geomPublishInFather(initLog.debug, geometrieSaine, long_p1, 'long_p1' )
 148     geomPublishInFather(initLog.debug, geometrieSaine, ep, 'ep' )
 149     geomPublishInFather(initLog.debug, geometrieSaine, long_coude, 'long_coude' )
 150     geomPublishInFather(initLog.debug, geometrieSaine, circ_g, 'circ_g' )
 151     geomPublishInFather(initLog.debug, geometrieSaine, circ_d, 'circ_d' )
 152     geomPublishInFather(initLog.debug, geometrieSaine, long_p2, 'long_p2' )
 153
 154     # --- face extremite tube (EXTUBE)
 155
 156     facesIds = geompy.GetShapesOnPlaneIDs(geometrieSaine, geompy.ShapeType["FACE"], OZ, GEOM.ST_ON)
 157     EXTUBE = geompy.CreateGroup(geometrieSaine, geompy.ShapeType["FACE"])
 158     geompy.UnionIDs(EXTUBE, facesIds)
 159     geomPublishInFather(initLog.debug, geometrieSaine, EXTUBE, 'EXTUBE' )
 160
 161     # --- edge bord extremite tube (BORDTU)
 162
 163     edge1Ids = geompy.GetShapesOnPlaneIDs(geometrieSaine, geompy.ShapeType["EDGE"], OZ, GEOM.ST_ON)
 164     edge2Ids = geompy.GetShapesOnCylinderIDs(geometrieSaine, geompy.ShapeType["EDGE"], OZ, de/2. -epais, GEOM.ST_ON)
 165     edgesIds = []
 166     for edge in edge1Ids:
 167       if edge in edge2Ids:
 168         edgesIds.append(edge)
 169     BORDTU = geompy.CreateGroup(geometrieSaine, geompy.ShapeType["EDGE"])
 170     geompy.UnionIDs(BORDTU, edgesIds)
 171     geomPublishInFather(initLog.debug, geometrieSaine, BORDTU, 'BORDTU' )
 172
 173     # --- face origine tube (CLGV)
 174
 175     pp2 = geompy.MakeTranslationVectorDistance(P2, Rotation_2, 10)
 176     vec2 = geompy.MakeVector(P2, pp2)
 177     #geomPublish(initLog.debug, vec2, 'vec2')
 178     facesIds = geompy.GetShapesOnPlaneIDs(geometrieSaine, geompy.ShapeType["FACE"], vec2, GEOM.ST_ON)
 179     CLGV = geompy.CreateGroup(geometrieSaine, geompy.ShapeType["FACE"])
 180     geompy.UnionIDs(CLGV, facesIds)
 181     geomPublishInFather(initLog.debug, geometrieSaine, CLGV, 'CLGV' )
 182
 183     # --- peau tube interieur (PEAUINT)
 184
 185     extru1 = geompy.MakePrismVecH(diskint, OZ, l_tube_p1)
 186     revol1 = geompy.MakeRevolution(diskint, axe, angleCoude*math.pi/180.0)
 187     rot1 = geompy.MakeRotation(diskint, axe, angleCoude*math.pi/180.0)
 188     extru2 = geompy.MakePrismVecH(rot1, Rotation_2, -l_tube_p2)
 189     interne = geompy.MakeFuse(extru1, revol1)
 190     interne = geompy.MakeFuse(extru2, interne)
 191     geomPublish(initLog.debug, interne, 'interne')
 192     facesIds = geompy.GetShapesOnShapeIDs(interne, geometrieSaine, geompy.ShapeType["FACE"], GEOM.ST_ONIN)
 193     PEAUINT = geompy.CreateGroup(geometrieSaine, geompy.ShapeType["FACE"])
 194     geompy.UnionIDs(PEAUINT, facesIds)
 195     geomPublishInFather(initLog.debug, geometrieSaine, PEAUINT, 'PEAUINT' )
 196
 197     # --- peau tube exterieur (PEAUEXT)
 198
 199     Disk_3 = geompy.MakeDiskPntVecR(centre, OZ, de/2. +epais)
 200     extru1 = geompy.MakePrismVecH(Disk_3, OZ, l_tube_p1)
 201     revol1 = geompy.MakeRevolution(Disk_3, axe, angleCoude*math.pi/180.0)
 202     rot1 = geompy.MakeRotation(Disk_3, axe, angleCoude*math.pi/180.0)
 203     extru2 = geompy.MakePrismVecH(rot1, Rotation_2, -l_tube_p2)
 204     externe = geompy.MakeFuse(extru1, revol1)
 205     externe = geompy.MakeFuse(extru2, externe)
 206     geomPublish(initLog.debug, externe, 'externe')
 207     facesIds = geompy.GetShapesOnShapeIDs(externe, geometrieSaine, geompy.ShapeType["FACE"], GEOM.ST_ON)
 208     PEAUEXT = geompy.CreateGroup(geometrieSaine, geompy.ShapeType["FACE"])
 209     geompy.UnionIDs(PEAUEXT, facesIds)
 210     geomPublishInFather(initLog.debug, geometrieSaine, PEAUEXT, 'PEAUEXT' )
 211
 212     # --- solide sain
 213
 214     volIds = geompy.SubShapeAllIDs(geometrieSaine, geompy.ShapeType["SOLID"])
 215     COUDE = geompy.CreateGroup(geometrieSaine, geompy.ShapeType["SOLID"])
 216     geompy.UnionIDs(COUDE, volIds)
 217     geomPublishInFather(initLog.debug, geometrieSaine, COUDE, 'COUDSAIN' )
 218
 219     geometriesSaines = [geometrieSaine, long_p1, ep, long_coude, circ_g, circ_d, long_p2, P1, P2, EXTUBE, BORDTU, CLGV, PEAUINT, PEAUEXT, COUDE]
 220
 221     return geometriesSaines
 222
 223   # ---------------------------------------------------------------------------
 224   def setParamMaillageSain(self):
 225     self.meshParams = dict(n_long_p1    = 16,
 226                            n_ep         = 3,
 227                            n_long_coude = 15,
 228                            n_circ_g     = 20,
 229                            n_circ_d     = 20,
 230                            n_long_p2    = 12)
 231
 232   # ---------------------------------------------------------------------------
 233   def genereMaillageSain(self, geometriesSaines, meshParams):
 234     logging.info("genereMaillageSain %s", self.nomCas)
 235
 236     geometrieSaine = geometriesSaines[0]
 237     long_p1        = geometriesSaines[1]
 238     ep             = geometriesSaines[2]
 239     long_coude     = geometriesSaines[3]
 240     circ_g         = geometriesSaines[4]
 241     circ_d         = geometriesSaines[5]
 242     long_p2        = geometriesSaines[6]
 243     P1             = geometriesSaines[7]
 244     P2             = geometriesSaines[8]
 245     EXTUBE         = geometriesSaines[9]
 246     BORDTU         = geometriesSaines[10]
 247     CLGV           = geometriesSaines[11]
 248     PEAUINT        = geometriesSaines[12]
 249     PEAUEXT        = geometriesSaines[13]
 250     COUDE          = geometriesSaines[14]
 251
 252     n_long_p1    = meshParams['n_long_p1']
 253     n_ep         = meshParams['n_ep']
 254     n_long_coude = meshParams['n_long_coude']
 255     n_circ_g     = meshParams['n_circ_g']
 256     n_circ_d     = meshParams['n_circ_d']
 257     n_long_p2    = meshParams['n_long_p2']
 258
 259     maillageSain = smesh.Mesh(geometrieSaine)
 260
 261     algo3d = maillageSain.Hexahedron()
 262     algo2d = maillageSain.Quadrangle()
 263     smesh.SetName(algo3d, "algo3d_maillageSain")
 264     smesh.SetName(algo2d, "algo2d_maillageSain")
 265
 266     algo1d_long_p1 = maillageSain.Segment(geom=long_p1)
 267     hypo1d_long_p1 = algo1d_long_p1.NumberOfSegments(n_long_p1)
 268     smesh.SetName(algo1d_long_p1, "algo1d_long_p1")
 269     smesh.SetName(hypo1d_long_p1, "hypo1d_long_p1")
 270
 271     algo1d_ep = maillageSain.Segment(geom=ep)
 272     hypo1d_ep = algo1d_ep.NumberOfSegments(n_ep)
 273     smesh.SetName(algo1d_ep, "algo1d_ep")
 274     smesh.SetName(hypo1d_ep, "hypo1d_ep")
 275
 276     algo1d_long_coude = maillageSain.Segment(geom=long_coude)
 277     hypo1d_long_coude = algo1d_long_coude.NumberOfSegments(n_long_coude)
 278     smesh.SetName(algo1d_long_coude, "algo1d_long_coude")
 279     smesh.SetName(hypo1d_long_coude, "hypo1d_long_coude")
 280
 281     algo1d_circ_g = maillageSain.Segment(geom=circ_g)
 282     hypo1d_circ_g = algo1d_circ_g.NumberOfSegments(n_circ_g)
 283     smesh.SetName(algo1d_circ_g, "algo1d_circ_g")
 284     smesh.SetName(hypo1d_circ_g, "hypo1d_circ_g")
 285
 286     algo1d_circ_d = maillageSain.Segment(geom=circ_d)
 287     hypo1d_circ_d = algo1d_circ_d.NumberOfSegments(n_circ_d)
 288     smesh.SetName(algo1d_circ_d, "algo1d_circ_d")
 289     smesh.SetName(hypo1d_circ_d, "hypo1d_circ_d")
 290
 291     algo1d_long_p2 = maillageSain.Segment(geom=long_p2)
 292     hypo1d_long_p2 = algo1d_long_p2.NumberOfSegments(n_long_p2)
 293     smesh.SetName(algo1d_long_p2, "algo1d_long_p2")
 294     smesh.SetName(hypo1d_long_p2, "hypo1d_long_p2")
 295
 296     is_done = maillageSain.Compute()
 297     text = "maillageSain.Compute"
 298     if is_done:
 299       logging.info(text+" OK")
 300     else:
 301       text = "Erreur au calcul du maillage.\n" + text
 302       logging.info(text)
 303       raise Exception(text)
 304
 305     _ = maillageSain.GroupOnGeom(P1,'P1',SMESH.NODE)
 306     _ = maillageSain.GroupOnGeom(P2,'P2',SMESH.NODE)
 307     _ = maillageSain.GroupOnGeom(EXTUBE,'EXTUBE',SMESH.FACE)
 308     _ = maillageSain.GroupOnGeom(BORDTU,'BORDTU',SMESH.EDGE)
 309     _ = maillageSain.GroupOnGeom(CLGV,'CLGV',SMESH.FACE)
 310     _ = maillageSain.GroupOnGeom(PEAUINT,'PEAUINT',SMESH.FACE)
 311     _ = maillageSain.GroupOnGeom(PEAUEXT,'PEAUEXT',SMESH.FACE)
 312     _ = maillageSain.GroupOnGeom(COUDE,'COUDSAIN',SMESH.VOLUME)
 313
 314     return [maillageSain, True] # True : maillage hexa
 315
 316   # ---------------------------------------------------------------------------
 317   def setParamShapeFissure(self):
 318     """
 319     paramètres de la fissure pour le tuyau coude
 320     profondeur  : 0 < profondeur <= épaisseur
 321     rayonPipe   : rayon du pipe correspondant au maillage rayonnant
 322     lenSegPipe  : longueur des mailles rayonnantes le long du fond de fissure (= rayonPipe par défaut)
 323     azimut      : entre 0 et 360°
 324     alpha       : 0 < alpha < angleCoude
 325     longueur    : <=2*profondeur ==> force une fissure elliptique (longueur/profondeur = grand axe/petit axe).
 326     orientation : 0° : longitudinale, 90° : circonférentielle, autre : uniquement fissures elliptiques
 327     lgInfluence : distance autour de la shape de fissure a remailler (si 0, pris égal à profondeur. A ajuster selon le maillage)
 328     elliptique  : True : fissure elliptique (longueur/profondeur = grand axe/petit axe); False : fissure longue (fond de fissure de profondeur constante, demi-cercles aux extrémites)
 329     pointIn_x   : optionnel coordonnées x d'un point dans le solide, pas trop loin du centre du fond de fissure (idem y,z)
 330     externe     : True : fissure face externe, False : fissure face interne
 331     """
 332     logging.info("setParamShapeFissure %s", self.nomCas)
 333     self.shapeFissureParams = dict(profondeur  = 10,
 334                                    rayonPipe   = 2.5,
 335                                    lenSegPipe  = 2.5,
 336                                    azimut      = 160,
 337                                    alpha       = 20,
 338                                    longueur    = 400,
 339                                    orientation = 90,
 340                                    lgInfluence = 50,
 341                                    elliptique  = False,
 342                                    externe     = True)
 343
 344   # ---------------------------------------------------------------------------
 345   def genereShapeFissure( self, geometriesSaines, geomParams, shapeFissureParams, \
 346                                 mailleur="MeshGems"):
 347     logging.info("genereShapeFissure %s", self.nomCas)
 348     logging.info("shapeFissureParams %s", shapeFissureParams)
 349
 350     angleCoude = geomParams['angleCoude']
 351     r_cintr    = geomParams['r_cintr']
 352     l_tube_p1  = geomParams['l_tube_p1']
 353     l_tube_p2  = geomParams['l_tube_p2']
 354     epais      = geomParams['epais']
 355     de         = geomParams['de']
 356
 357     profondeur  = shapeFissureParams['profondeur']
 358     azimut      = shapeFissureParams['azimut']
 359     alpha       = shapeFissureParams['alpha']
 360     longueur    = shapeFissureParams['longueur']
 361     orientation = shapeFissureParams['orientation']
 362     externe     = shapeFissureParams['externe']
 363     lgInfluence = shapeFissureParams['lgInfluence']
 364     self.elliptique  = False
 365     if 'elliptique' in shapeFissureParams:
 366       self.elliptique = shapeFissureParams['elliptique']
 367
 368
 369
 370     azimut = -azimut # axe inverse / ASCOUF
 371     axe = geompy.MakeTranslation(OY, -r_cintr, 0, -l_tube_p1)
 372     geomPublish(initLog.debug, axe,"axe")
 373
 374     if not lgInfluence:
 375       lgInfluence = profondeur
 376
 377     if longueur > 2*profondeur:
 378       self.fissureLongue=True
 379     else:
 380       self.fissureLongue=False
 381       self.elliptique = True
 382
 383     self.circonferentielle = False
 384     self.longitudinale = False
 385     if self.fissureLongue and not self.elliptique:
 386       if abs(orientation) < 45 :
 387         self.longitudinale = True
 388       else:
 389         self.circonferentielle = True
 390
 391     nbp1 = 10
 392     if self.circonferentielle:
 393       if externe:
 394         dp = -1.0
 395         raybor = de/2.
 396         rayint = raybor - profondeur
 397         rayext = raybor + profondeur/5.0
 398       else:
 399         dp = 1.0
 400         raybor = de/2. - epais
 401         rayint = raybor + profondeur
 402         rayext = raybor - profondeur/5.0
 403       lgfond = longueur -2*profondeur
 404       angle = lgfond/(2*raybor)
 405       pb = geompy.MakeVertex(raybor, 0, 0)
 406       pi = geompy.MakeVertex(rayint, 0, 0)
 407       pbl = geompy.MakeRotation(pb, OZ, angle)
 408       pbr = geompy.MakeRotation(pb, OZ, -angle)
 409       geomPublish(initLog.debug, pbl,"pbl")
 410       geomPublish(initLog.debug, pbr,"pbr")
 411       pal = geompy.MakeTranslationVector(pbl, OZ)
 412       par = geompy.MakeTranslationVector(pbr, OZ)
 413       axl = geompy.MakeVector(pbl,pal)
 414       axr = geompy.MakeVector(pbr,par)
 415       pil = geompy.MakeRotation(pi, OZ, angle)
 416       pir = geompy.MakeRotation(pi, OZ, -angle)
 417       points = []
 418       nbp = 3*nbp1
 419       for i in range(nbp):
 420         angi = dp*(nbp -i)*(2.0*math.pi/3.0)/nbp
 421         pt = geompy.MakeRotation(pil, axl, angi)
 422         points.append(pt)
 423       for i in range(nbp):
 424         angi = angle -2.0*i*angle/nbp
 425         pt = geompy.MakeRotation(pi, OZ, angi)
 426         points.append(pt)
 427       for i in range(nbp+1):
 428         angi = -dp*i*(2.0*math.pi/3.0)/nbp
 429         pt = geompy.MakeRotation(pir, axr, angi)
 430         points.append(pt)
 431       for i, pt in enumerate(points):
 432         pt = geompy.MakeRotation(pt, OZ, azimut*math.pi/180.)
 433         pt = geompy.MakeTranslation(pt, 0, 0, -l_tube_p1)
 434         pt = geompy.MakeRotation(pt, axe, alpha*math.pi/180.)
 435         points[i] = pt
 436       wire0 = geompy.MakeInterpol(points[0:nbp+1])
 437       wire1 = geompy.MakeInterpol(points[nbp:2*nbp+1])
 438       wire2 = geompy.MakeInterpol(points[2*nbp:3*nbp+1])
 439       #wiretube = geompy.MakeInterpol(points)
 440       wiretube=geompy.MakeWire([wire0,wire1,wire2])
 441       geomPublish(initLog.debug, wiretube,"wiretube")
 442
 443       pe = geompy.MakeVertex(rayext, 0, 0)
 444       pe = geompy.MakeRotation(pe, OZ, azimut*math.pi/180.)
 445       pe = geompy.MakeTranslation(pe, 0, 0, -l_tube_p1)
 446       pe = geompy.MakeRotation(pe, axe, alpha*math.pi/180.)
 447
 448       arce = geompy.MakeArc(points[0], pe, points[-1])
 449       geomPublish(initLog.debug, arce,"arce")
 450
 451       facefiss = geompy.MakeFaceWires([arce, wiretube], 1)
 452       geomPublish(initLog.debug,  facefiss, 'facefissPlace' )
 453
 454       pc = geompy.MakeVertex((raybor + rayint)/2.0, 0, 0)
 455       centre = geompy.MakeRotation(pc, OZ, azimut*math.pi/180.)
 456       centre = geompy.MakeTranslation(centre, 0, 0, -l_tube_p1)
 457       centre = geompy.MakeRotation(centre, axe, alpha*math.pi/180.)
 458       geomPublish(initLog.debug,  centre, 'centrefissPlace' )
 459
 460       wiretube = geompy.GetInPlace(facefiss, wiretube)
 461       geomPublish(initLog.debug, wiretube, 'wiretubePlace' )
 462       try:
 463         edgetube = geompy.MakeEdgeWire(wiretube)
 464         geomPublish(initLog.debug, edgetube,"edgetube")
 465       except:
 466         logging.debug("erreur MakeEdgeWire sur fond de fissure, on fait sans")
 467         edgetube = None
 468
 469     # ---------------------------------------------------------
 470
 471     elif self.longitudinale:
 472       if externe:
 473         raybor = de/2.
 474         dp = -1.0
 475       else:
 476         raybor = de/2. - epais
 477         dp = +1.0
 478       prof = dp * profondeur
 479       lgfond = longueur -2*profondeur
 480       cosaz = math.cos(azimut*math.pi/180.)
 481       sinaz = math.sin(azimut*math.pi/180.)
 482       alfrd = alpha*math.pi/180.
 483       rayxy = r_cintr + raybor*cosaz
 484       angle = lgfond/(2.*rayxy)
 485       logging.debug("longueur: %s, angle: %s, rayon: %s",lgfond, angle, rayxy)
 486       pb = geompy.MakeVertex(raybor*cosaz, raybor*sinaz, -l_tube_p1, "pb")
 487       pi = geompy.MakeTranslation(pb, prof*cosaz, prof*sinaz, 0., "pi")
 488       pbv = geompy.MakeTranslation(pb, -sinaz, cosaz, 0., "pbv")
 489       axb  = geompy.MakeVector(pb,pbv, "axb")
 490       pbl = geompy.MakeRotation(pb, axe, alfrd -angle, "pbl")
 491       pbr = geompy.MakeRotation(pb, axe, alfrd +angle, "pbr")
 492       axl = geompy.MakeRotation(axb, axe, alfrd -angle, "axl")
 493       axr = geompy.MakeRotation(axb, axe, alfrd +angle, "axr")
 494       pil = geompy.MakeRotation(pi, axe, alfrd -angle, "pil")
 495       pir = geompy.MakeRotation(pi, axe, alfrd +angle, "pir")
 496
 497       curves = []
 498
 499       points = []
 500       nbp = 3*nbp1
 501       xs = []
 502       totx = 0
 503       for i in range(nbp+2):
 504         x = math.sin(i*math.pi/(nbp+1)) # fonction de répartition des points : distance relative
 505         x2 = x*x
 506         totx += x2
 507         xs.append(totx)
 508         logging.debug("x2: %s, totx: %s", x2, totx)
 509       for i in range(nbp+1):
 510         #posi = nbp -i             # répartition équidistante des points sur la courbe
 511         posi = nbp*(1 -xs[i]/totx) # points plus resserrés aux extrémités de la courbe
 512         angi = -dp*posi*(5.0*math.pi/8.0)/nbp
 513         pt = geompy.MakeRotation(pil, axl, angi)
 514         points.append(pt)
 515       curves.append(geompy.MakeInterpol(points))
 516       point0 = points[0]
 517       geomPublish(initLog.debug, curves[-1],"curve0")
 518 #      for i, pt in enumerate(points):
 519 #        name = "point%d"%i
 520 #        geomPublishInFather(initLog.debug,curves[-1], pt, name)
 521
 522       points = []
 523       nbp = 3*nbp1
 524       xs =[]
 525       totx = 0
 526       for i in range(nbp+1):
 527         x = math.sin(i*math.pi/nbp)
 528         #x = 1.0 # répartition équidistante des points sur la courbe
 529         x2 = x*x # points plus resserrés aux extrémités de la courbe
 530         totx += x2
 531         xs.append(totx)
 532         logging.debug("x2: %s, totx: %s", x2, totx)
 533       for i in range(nbp):
 534         angi = alfrd -angle +2.0*angle*xs[i]/totx
 535         pt = geompy.MakeRotation(pi, axe, angi)
 536         points.append(pt)
 537       curves.append(geompy.MakeInterpol(points))
 538       geomPublish(initLog.debug, curves[-1],"curve1")
 539 #      for i, pt in enumerate(points):
 540 #        name = "point%d"%i
 541 #        geomPublishInFather(initLog.debug,curves[-1], pt, name)
 542
 543       points = []
 544       nbp = 3*nbp1
 545       xs = []
 546       totx = 0
 547       for i in range(nbp+2):
 548         x = math.sin(i*math.pi/(nbp+1))
 549         x2 = x*x
 550         totx += x2
 551         xs.append(totx)
 552         logging.debug("x2: %s, totx: %s", x2, totx)
 553       for i in range(nbp+1):
 554         #posi = nbp -i        # répartition équidistante des points sur la courbe
 555         posi = nbp*xs[i]/totx # points plus resserrés aux extrémités de la courbe
 556         angi = dp*posi*(5.0*math.pi/8.0)/nbp
 557         pt = geompy.MakeRotation(pir, axr, angi)
 558         points.append(pt)
 559       curves.append(geompy.MakeInterpol(points))
 560       point1 = points[-1]
 561       geomPublish(initLog.debug, curves[-1],"curve2")
 562 #      for i, pt in enumerate(points):
 563 #        name = "point%d"%i
 564 #        geomPublishInFather(initLog.debug,curves[-1], pt, name)
 565
 566       wiretube = geompy.MakeWire(curves)
 567       geomPublish(initLog.debug, wiretube,"wiretube")
 568       try:
 569         edgetube = geompy.MakeEdgeWire(wiretube)
 570         geomPublish(initLog.debug, edgetube,"edgetube")
 571       except:
 572         logging.debug("erreur MakeEdgeWire sur fond de fissure, on fait sans")
 573         edgetube = None
 574
 575       pts = []
 576       pts.append(point0)
 577       dpr = prof*math.cos(5.0*math.pi/8.0)
 578       pe = geompy.MakeTranslation(pb, dpr*cosaz, dpr*sinaz, 0., "pe")
 579       for i in range(nbp):
 580         angi = alfrd -angle +2.0*i*angle/nbp
 581         pt = geompy.MakeRotation(pe, axe, angi)
 582         pts.append(pt)
 583       pts.append(point1)
 584       arce = geompy.MakeInterpol(pts)
 585       geomPublish(initLog.debug, arce,"arce")
 586
 587       facefiss = geompy.MakeFaceWires([arce, wiretube], 0)
 588       geomPublish(initLog.debug,  facefiss, 'facefissPlace' )
 589
 590       pc = geompy.MakeTranslation(pb, 0.5*prof*cosaz, 0.5*prof*sinaz, 0.)
 591       centre = geompy.MakeRotation(pc, axe, alfrd)
 592       geomPublish(initLog.debug,  centre, 'centrefissPlace' )
 593
 594       edges = geompy.ExtractShapes(facefiss, geompy.ShapeType["EDGE"], True)
 595       edgesTriees, minl, maxl = sortEdges(edges)
 596       edges = edgesTriees[:-1] # la plus grande correspond à arce, on l'elimine
 597       wiretube = geompy.MakeWire(edges)
 598       #wiretube = edgesTriees[-1]
 599       geomPublish(initLog.debug, wiretube, 'wiretubePlace' )
 600
 601     # ---------------------------------------------------------
 602
 603     else: # fissure elliptique, longue ou courte
 604       if externe:
 605         raybor = de/2.
 606         dp = -1.0
 607       else:
 608         raybor = de/2. - epais
 609         dp = +1.0
 610       prof = dp * profondeur
 611       cosaz = math.cos(azimut*math.pi/180.)
 612       sinaz = math.sin(azimut*math.pi/180.)
 613       alfrd = alpha*math.pi/180.
 614       pb = geompy.MakeVertex(raybor*cosaz, raybor*sinaz, -l_tube_p1, "pb")
 615       pi = geompy.MakeTranslation(pb, prof*cosaz, prof*sinaz, 0., "pi")
 616       pbv = geompy.MakeTranslation(pb, -profondeur*sinaz, profondeur*cosaz, 0., "pbv")
 617       ayb  = geompy.MakeVector(pb,pbv, "ayb")
 618       pb0 = geompy.MakeRotation(pb, axe, alfrd, "pb0")
 619       ay0 = geompy.MakeRotation(ayb, axe, alfrd, "ay0")
 620       pi0 = geompy.MakeRotation(pi, axe, alfrd, "pi0")
 621       az_ = geompy.MakeVector(pi0, pb0, "az_")
 622       az0 = geompy.MakeTranslationVector(az_, az_, "az0") #normale sortante
 623       ax0 = geompy.MakeRotation(ay0, az0, -math.pi/2.0, "ax0")
 624       ax1 = geompy.MakeRotation(ax0, az0, orientation*math.pi/180., "ax1")
 625       ay1 = geompy.MakeRotation(ay0, az0, orientation*math.pi/180., "ay1")
 626       originLCS = geompy.MakeMarker(0, 0, 0, 1, 0, 0, 0, 1, 0, "originLCS")
 627       coo = geompy.PointCoordinates(pb0)
 628       cox = geompy.VectorCoordinates(ax1)
 629       coy = geompy.VectorCoordinates(ay1)
 630       localLCS = geompy.MakeMarker(coo[0], coo[1], coo[2], cox[0], cox[1], cox[2], coy[0], coy[1], coy[2], "localLCS")
 631
 632       pco = geompy.MakeVertex(0, 0, -profondeur, "pco")
 633       pao = geompy.MakeRotation(pco, OY, 0.6*math.pi, "pao")
 634       pbo = geompy.MakeRotation(pco, OY, -0.6*math.pi, "pbo")
 635       pce = geompy.MakeVertex(0, 0, 0.1*profondeur,"pce")
 636       arcoo = geompy.MakeArc(pao, pco, pbo, "arcoo")
 637       linoo = geompy.MakeArc(pao, pce, pbo, "linoo")
 638       scalex = longueur/profondeur
 639       arco =geompy.MakeScaleAlongAxes(arcoo, O, scalex, 1., 1., "arco")
 640       lino =geompy.MakeScaleAlongAxes(linoo, O, scalex, 1., 1., "lino")
 641       arci = geompy.MakePosition(arco, originLCS, localLCS, "arci")
 642       arce = geompy.MakePosition(lino, originLCS, localLCS, "arce")
 643       facefiss = geompy.MakeFaceWires([arce, arci], 0)
 644       geomPublish(initLog.debug,  facefiss, 'facefissPlace' )
 645       edges = geompy.ExtractShapes(facefiss, geompy.ShapeType["EDGE"], True)
 646       edgesTriees, minl, maxl = sortEdges(edges)
 647       edgetube = edgesTriees[-1] # la plus grande correspond à arci
 648       wiretube = edgetube
 649
 650       pc = geompy.MakeTranslation(pb, 0.5*prof*cosaz, 0.5*prof*sinaz, 0.)
 651       centre = geompy.MakeRotation(pc, axe, alfrd)
 652       geomPublish(initLog.debug,  centre, 'centrefissPlace' )
 653
 654     coordsNoeudsFissure = genereMeshCalculZoneDefaut(facefiss, profondeur/2. ,profondeur, \
 655                                                      mailleur)
 656
 657     return [facefiss, centre, lgInfluence, coordsNoeudsFissure, wiretube, edgetube]
 658
 659   # ---------------------------------------------------------------------------
 660   def setParamMaillageFissure(self):
 661     """
 662     Paramètres du maillage de la fissure pour le tuyau coudé
 663     Voir également setParamShapeFissure, paramètres rayonPipe et lenSegPipe.
 664     nbSegRad = nombre de couronnes
 665     nbSegCercle = nombre de secteurs
 666     areteFaceFissure = taille cible de l'arête des triangles en face de fissure.
 667     """
 668     self.maillageFissureParams = dict(nomRep        = os.curdir,
 669                                       nomFicSain    = self.nomCas,
 670                                       nomFicFissure = 'fissure_' + self.nomCas,
 671                                       nbsegRad      = 5,
 672                                       nbsegCercle   = 6,
 673                                       areteFaceFissure = 5)
 674
 675   # ---------------------------------------------------------------------------
 676   def genereZoneDefaut(self, geometriesSaines, maillagesSains, shapesFissure, shapeFissureParams, maillageFissureParams):
 677     elementsDefaut = creeZoneDefautDansObjetSain(geometriesSaines, maillagesSains, shapesFissure, shapeFissureParams, maillageFissureParams)
 678     return elementsDefaut
 679
 680   # ---------------------------------------------------------------------------
 681   def genereMaillageFissure(self, geometriesSaines, maillagesSains, \
 682                             shapesFissure, shapeFissureParams, \
 683                             maillageFissureParams, elementsDefaut, step, \
 684                             mailleur="MeshGems"):
 685     maillageFissure = construitFissureGenerale(shapesFissure, shapeFissureParams, \
 686                                                maillageFissureParams, elementsDefaut, \
 687                                                step, mailleur)
 688     return maillageFissure
 689
 690   # ---------------------------------------------------------------------------
 691   def setReferencesMaillageFissure(self):
 692     self.referencesMaillageFissure = dict(Entity_Node            = 77917,
 693                                           Entity_Quad_Edge       = 975,
 694                                           Entity_Quad_Triangle   = 2182,
 695                                           Entity_Quad_Quadrangle = 6842,
 696                                           Entity_Quad_Tetra      = 20135,
 697                                           Entity_Quad_Hexa       = 8994,
 698                                           Entity_Quad_Penta      = 972,
 699                                           Entity_Quad_Pyramid    = 1038)