src/Tools/blocFissure/gmu/fissureCoude.py

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