src/Tools/blocFissure/gmu/fissureCoude.py

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