src/Tools/blocFissure/gmu/fissureCoude.py

   1 # -*- coding: utf-8 -*-
   2 # Copyright (C) 2014-2021  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 os
  23
  24 import logging
  25 import math
  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
 263     algo3d = maillageSain.Hexahedron()
 264     algo2d = maillageSain.Quadrangle()
 265     putName(algo3d, "{}_3d_maillageSain".format(self.mailleur2d3d()), i_pref=self.numeroCas)
 266     putName(algo2d, "{}_2d_maillageSain".format(self.mailleur2d3d()), i_pref=self.numeroCas)
 267
 268     algo1d_long_p1 = maillageSain.Segment(geom=long_p1)
 269     hypo1d_long_p1 = algo1d_long_p1.NumberOfSegments(n_long_p1)
 270     putName(algo1d_long_p1, "algo1d_long_p1", i_pref=self.numeroCas)
 271     putName(hypo1d_long_p1, "hypo1d_long_p1", i_pref=self.numeroCas)
 272
 273     algo1d_ep = maillageSain.Segment(geom=ep)
 274     hypo1d_ep = algo1d_ep.NumberOfSegments(n_ep)
 275     putName(algo1d_ep, "algo1d_ep", i_pref=self.numeroCas)
 276     putName(hypo1d_ep, "hypo1d_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(algo1d_long_coude, "algo1d_long_coude", i_pref=self.numeroCas)
 281     putName(hypo1d_long_coude, "hypo1d_long_coude", i_pref=self.numeroCas)
 282
 283     algo1d_circ_g = maillageSain.Segment(geom=circ_g)
 284     hypo1d_circ_g = algo1d_circ_g.NumberOfSegments(n_circ_g)
 285     putName(algo1d_circ_g, "algo1d_circ_g", i_pref=self.numeroCas)
 286     putName(hypo1d_circ_g, "hypo1d_circ_g", i_pref=self.numeroCas)
 287
 288     algo1d_circ_d = maillageSain.Segment(geom=circ_d)
 289     hypo1d_circ_d = algo1d_circ_d.NumberOfSegments(n_circ_d)
 290     putName(algo1d_circ_d, "algo1d_circ_d", i_pref=self.numeroCas)
 291     putName(hypo1d_circ_d, "hypo1d_circ_d", i_pref=self.numeroCas)
 292
 293     algo1d_long_p2 = maillageSain.Segment(geom=long_p2)
 294     hypo1d_long_p2 = algo1d_long_p2.NumberOfSegments(n_long_p2)
 295     putName(algo1d_long_p2, "algo1d_long_p2", i_pref=self.numeroCas)
 296     putName(hypo1d_long_p2, "hypo1d_long_p2", i_pref=self.numeroCas)
 297
 298     is_done = maillageSain.Compute()
 299     text = "maillageSain.Compute"
 300     if is_done:
 301       logging.info(text+" OK")
 302     else:
 303       text = "Erreur au calcul du maillage.\n" + text
 304       logging.info(text)
 305       raise Exception(text)
 306
 307     _ = maillageSain.GroupOnGeom(P1,'P1',SMESH.NODE)
 308     _ = maillageSain.GroupOnGeom(P2,'P2',SMESH.NODE)
 309     _ = maillageSain.GroupOnGeom(EXTUBE,'EXTUBE',SMESH.FACE)
 310     _ = maillageSain.GroupOnGeom(BORDTU,'BORDTU',SMESH.EDGE)
 311     _ = maillageSain.GroupOnGeom(CLGV,'CLGV',SMESH.FACE)
 312     _ = maillageSain.GroupOnGeom(PEAUINT,'PEAUINT',SMESH.FACE)
 313     _ = maillageSain.GroupOnGeom(PEAUEXT,'PEAUEXT',SMESH.FACE)
 314     _ = maillageSain.GroupOnGeom(COUDE,'COUDSAIN',SMESH.VOLUME)
 315
 316     return [maillageSain, True] # True : maillage hexa
 317
 318   # ---------------------------------------------------------------------------
 319   def setParamShapeFissure(self):
 320     """
 321     paramètres de la fissure pour le tuyau coude
 322     profondeur  : 0 < profondeur <= épaisseur
 323     rayonPipe   : rayon du pipe correspondant au maillage rayonnant
 324     lenSegPipe  : longueur des mailles rayonnantes le long du fond de fissure (= rayonPipe par défaut)
 325     azimut      : entre 0 et 360°
 326     alpha       : 0 < alpha < angleCoude
 327     longueur    : <=2*profondeur ==> force une fissure elliptique (longueur/profondeur = grand axe/petit axe).
 328     orientation : 0° : longitudinale, 90° : circonférentielle, autre : uniquement fissures elliptiques
 329     lgInfluence : distance autour de la shape de fissure a remailler (si 0, pris égal à profondeur. A ajuster selon le maillage)
 330     elliptique  : True : fissure elliptique (longueur/profondeur = grand axe/petit axe); False : fissure longue (fond de fissure de profondeur constante, demi-cercles aux extrémites)
 331     pointIn_x   : optionnel coordonnées x d'un point dans le solide, pas trop loin du centre du fond de fissure (idem y,z)
 332     externe     : True : fissure face externe, False : fissure face interne
 333     """
 334     logging.info("setParamShapeFissure %s", self.nomCas)
 335     self.shapeFissureParams = dict(profondeur  = 10,
 336                                    rayonPipe   = 2.5,
 337                                    lenSegPipe  = 2.5,
 338                                    azimut      = 160,
 339                                    alpha       = 20,
 340                                    longueur    = 400,
 341                                    orientation = 90,
 342                                    lgInfluence = 50,
 343                                    elliptique  = False,
 344                                    externe     = True)
 345
 346   # ---------------------------------------------------------------------------
 347   def genereShapeFissure( self, geometriesSaines, geomParams, shapeFissureParams, \
 348                                 mailleur="MeshGems"):
 349     logging.info("genereShapeFissure %s", self.nomCas)
 350     logging.info("shapeFissureParams %s", shapeFissureParams)
 351
 352     r_cintr    = geomParams['r_cintr']
 353     l_tube_p1  = geomParams['l_tube_p1']
 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     azimut = -azimut # axe inverse / ASCOUF
 369     axe = geompy.MakeTranslation(OY, -r_cintr, 0, -l_tube_p1)
 370     geomPublish(initLog.debug, axe,"axe", self.numeroCas)
 371
 372     if not lgInfluence:
 373       lgInfluence = profondeur
 374
 375     if longueur > 2*profondeur:
 376       self.fissureLongue=True
 377     else:
 378       self.fissureLongue=False
 379       self.elliptique = True
 380
 381     self.circonferentielle = False
 382     self.longitudinale = False
 383     if self.fissureLongue and not self.elliptique:
 384       self.longitudinale = bool(abs(orientation) < 45)
 385       self.circonferentielle = not bool(abs(orientation) < 45)
 386
 387     nbp1 = 10
 388     if self.circonferentielle:
 389       if externe:
 390         dp = -1.0
 391         raybor = de/2.
 392         rayint = raybor - profondeur
 393         rayext = raybor + profondeur/5.0
 394       else:
 395         dp = 1.0
 396         raybor = de/2. - epais
 397         rayint = raybor + profondeur
 398         rayext = raybor - profondeur/5.0
 399       lgfond = longueur -2.*profondeur
 400       angle = lgfond/(2.*raybor)
 401       pb = geompy.MakeVertex(raybor, 0, 0)
 402       pi = geompy.MakeVertex(rayint, 0, 0)
 403       pbl = geompy.MakeRotation(pb, OZ, angle)
 404       pbr = geompy.MakeRotation(pb, OZ, -angle)
 405       geomPublish(initLog.debug, pbl,"pbl", self.numeroCas)
 406       geomPublish(initLog.debug, pbr,"pbr", self.numeroCas)
 407       pal = geompy.MakeTranslationVector(pbl, OZ)
 408       par = geompy.MakeTranslationVector(pbr, OZ)
 409       axl = geompy.MakeVector(pbl,pal)
 410       axr = geompy.MakeVector(pbr,par)
 411       pil = geompy.MakeRotation(pi, OZ, angle)
 412       pir = geompy.MakeRotation(pi, OZ, -angle)
 413       points = list()
 414       nbp = 3*nbp1
 415       for i_aux in range(nbp):
 416         angi = dp*(nbp - i_aux)*(2.0*math.pi/3.0)/nbp
 417         pt = geompy.MakeRotation(pil, axl, angi)
 418         points.append(pt)
 419       for i_aux in range(nbp):
 420         angi = angle -2.0*float(i_aux)*angle/nbp
 421         pt = geompy.MakeRotation(pi, OZ, angi)
 422         points.append(pt)
 423       for i_aux in range(nbp+1):
 424         angi = -dp*float(i_aux)*(2.0*math.pi/3.0)/nbp
 425         pt = geompy.MakeRotation(pir, axr, angi)
 426         points.append(pt)
 427       for i_aux, pt in enumerate(points):
 428         pt = geompy.MakeRotation(pt, OZ, azimut*math.pi/180.)
 429         pt = geompy.MakeTranslation(pt, 0, 0, -l_tube_p1)
 430         pt = geompy.MakeRotation(pt, axe, alpha*math.pi/180.)
 431         points[i_aux] = pt
 432       wire0 = geompy.MakeInterpol(points[0:nbp+1])
 433       wire1 = geompy.MakeInterpol(points[nbp:2*nbp+1])
 434       wire2 = geompy.MakeInterpol(points[2*nbp:3*nbp+1])
 435       #wiretube = geompy.MakeInterpol(points)
 436       wiretube=geompy.MakeWire([wire0,wire1,wire2])
 437       geomPublish(initLog.debug, wiretube,"wiretube", self.numeroCas)
 438
 439       pe = geompy.MakeVertex(rayext, 0, 0)
 440       pe = geompy.MakeRotation(pe, OZ, azimut*math.pi/180.)
 441       pe = geompy.MakeTranslation(pe, 0, 0, -l_tube_p1)
 442       pe = geompy.MakeRotation(pe, axe, alpha*math.pi/180.)
 443
 444       arce = geompy.MakeArc(points[0], pe, points[-1])
 445       geomPublish(initLog.debug, arce,"arce", self.numeroCas)
 446
 447       facefiss = geompy.MakeFaceWires([arce, wiretube], 1)
 448       geomPublish(initLog.debug, facefiss, 'facefissPlace', self.numeroCas )
 449
 450       pc = geompy.MakeVertex((raybor + rayint)/2.0, 0, 0)
 451       centre = geompy.MakeRotation(pc, OZ, azimut*math.pi/180.)
 452       centre = geompy.MakeTranslation(centre, 0, 0, -l_tube_p1)
 453       centre = geompy.MakeRotation(centre, axe, alpha*math.pi/180.)
 454       geomPublish(initLog.debug, centre, 'centrefissPlace', self.numeroCas )
 455
 456       wiretube = geompy.GetInPlace(facefiss, wiretube)
 457       geomPublish(initLog.debug, wiretube, 'wiretubePlace', self.numeroCas )
 458       try:
 459         edgetube = geompy.MakeEdgeWire(wiretube)
 460         geomPublish(initLog.debug, edgetube,"edgetube", self.numeroCas)
 461       except:
 462         logging.debug("erreur MakeEdgeWire sur fond de fissure, on fait sans")
 463         edgetube = None
 464
 465     # ---------------------------------------------------------
 466
 467     elif self.longitudinale:
 468       if externe:
 469         raybor = de/2.
 470         dp = -1.0
 471       else:
 472         raybor = de/2. - epais
 473         dp = +1.0
 474       prof = dp * profondeur
 475       lgfond = longueur -2*profondeur
 476       cosaz = math.cos(azimut*math.pi/180.)
 477       sinaz = math.sin(azimut*math.pi/180.)
 478       alfrd = alpha*math.pi/180.
 479       rayxy = r_cintr + raybor*cosaz
 480       angle = lgfond/(2.*rayxy)
 481       logging.debug("longueur: %s, angle: %s, rayon: %s",lgfond, angle, rayxy)
 482       pb = geompy.MakeVertex(raybor*cosaz, raybor*sinaz, -l_tube_p1, "pb")
 483       pi = geompy.MakeTranslation(pb, prof*cosaz, prof*sinaz, 0., "pi")
 484       pbv = geompy.MakeTranslation(pb, -sinaz, cosaz, 0., "pbv")
 485       axb  = geompy.MakeVector(pb,pbv, "axb")
 486       pbl = geompy.MakeRotation(pb, axe, alfrd -angle, "pbl")
 487       pbr = geompy.MakeRotation(pb, axe, alfrd +angle, "pbr")
 488       axl = geompy.MakeRotation(axb, axe, alfrd -angle, "axl")
 489       axr = geompy.MakeRotation(axb, axe, alfrd +angle, "axr")
 490       pil = geompy.MakeRotation(pi, axe, alfrd -angle, "pil")
 491       pir = geompy.MakeRotation(pi, axe, alfrd +angle, "pir")
 492
 493       curves = list()
 494
 495       points = list()
 496       nbp = 3*nbp1
 497       xs = list()
 498       totx = 0
 499       for i_aux in range(nbp+2):
 500         x = math.sin(float(i_aux)*math.pi/(nbp+1)) # fonction de répartition des points : distance relative
 501         x2 = x*x
 502         totx += x2
 503         xs.append(totx)
 504         logging.debug("x2: %s, totx: %s", x2, totx)
 505       for i_aux in range(nbp+1):
 506         #posi = nbp - i_aux             # répartition équidistante des points sur la courbe
 507         posi = nbp*(1 -xs[i_aux]/totx) # points plus resserrés aux extrémités de la courbe
 508         angi = -dp*posi*(5.0*math.pi/8.0)/nbp
 509         pt = geompy.MakeRotation(pil, axl, angi)
 510         points.append(pt)
 511       curves.append(geompy.MakeInterpol(points))
 512       point0 = points[0]
 513       geomPublish(initLog.debug, curves[-1],"curve0", self.numeroCas)
 514 #      for i_aux, pt in enumerate(points):
 515 #        name = "point%d"%i_aux
 516 #        geomPublishInFather(initLog.debug,curves[-1], pt, name)
 517
 518       points = list()
 519       nbp = 3*nbp1
 520       xs = list()
 521       totx = 0
 522       for i_aux in range(nbp+1):
 523         x = math.sin(float(i_aux)*math.pi/nbp)
 524         #x = 1.0 # répartition équidistante des points sur la courbe
 525         x2 = x*x # points plus resserrés aux extrémités de la courbe
 526         totx += x2
 527         xs.append(totx)
 528         logging.debug("x2: %s, totx: %s", x2, totx)
 529       for i_aux in range(nbp):
 530         angi = alfrd -angle +2.0*angle*xs[i_aux]/totx
 531         pt = geompy.MakeRotation(pi, axe, angi)
 532         points.append(pt)
 533       curves.append(geompy.MakeInterpol(points))
 534       geomPublish(initLog.debug, curves[-1],"curve1", self.numeroCas)
 535 #      for i_aux, pt in enumerate(points):
 536 #        name = "point%d"%i_aux
 537 #        geomPublishInFather(initLog.debug,curves[-1], pt, name)
 538
 539       points = list()
 540       nbp = 3*nbp1
 541       xs = list()
 542       totx = 0
 543       for i_aux in range(nbp+2):
 544         x = math.sin(float(i_aux)*math.pi/(nbp+1))
 545         x2 = x*x
 546         totx += x2
 547         xs.append(totx)
 548         logging.debug("x2: %s, totx: %s", x2, totx)
 549       for i_aux in range(nbp+1):
 550         #posi = nbp - i_aux        # répartition équidistante des points sur la courbe
 551         posi = nbp*xs[i_aux]/totx # points plus resserrés aux extrémités de la courbe
 552         angi = dp*posi*(5.0*math.pi/8.0)/nbp
 553         pt = geompy.MakeRotation(pir, axr, angi)
 554         points.append(pt)
 555       curves.append(geompy.MakeInterpol(points))
 556       point1 = points[-1]
 557       geomPublish(initLog.debug, curves[-1],"curve2", self.numeroCas)
 558 #      for i_aux, pt in enumerate(points):
 559 #        name = "point%d"%i_aux
 560 #        geomPublishInFather(initLog.debug,curves[-1], pt, name)
 561
 562       wiretube = geompy.MakeWire(curves)
 563       geomPublish(initLog.debug, wiretube,"wiretube", self.numeroCas)
 564       try:
 565         edgetube = geompy.MakeEdgeWire(wiretube)
 566         geomPublish(initLog.debug, edgetube,"edgetube", self.numeroCas)
 567       except:
 568         logging.debug("erreur MakeEdgeWire sur fond de fissure, on fait sans")
 569         edgetube = None
 570
 571       pts = list()
 572       pts.append(point0)
 573       dpr = prof*math.cos(5.0*math.pi/8.0)
 574       pe = geompy.MakeTranslation(pb, dpr*cosaz, dpr*sinaz, 0., "pe")
 575       for i_aux in range(nbp):
 576         angi = alfrd -angle +2.0*float(i_aux)*angle/nbp
 577         pt = geompy.MakeRotation(pe, axe, angi)
 578         pts.append(pt)
 579       pts.append(point1)
 580       arce = geompy.MakeInterpol(pts)
 581       geomPublish(initLog.debug, arce,"arce", self.numeroCas)
 582
 583       facefiss = geompy.MakeFaceWires([arce, wiretube], 0)
 584       geomPublish(initLog.debug, facefiss, 'facefissPlace', self.numeroCas)
 585
 586       pc = geompy.MakeTranslation(pb, 0.5*prof*cosaz, 0.5*prof*sinaz, 0.)
 587       centre = geompy.MakeRotation(pc, axe, alfrd)
 588       geomPublish(initLog.debug, centre, 'centrefissPlace', self.numeroCas)
 589
 590       edges = geompy.ExtractShapes(facefiss, geompy.ShapeType["EDGE"], True)
 591       edgesTriees, _, _ = sortEdges(edges)
 592       edges = edgesTriees[:-1] # la plus grande correspond à arce, on l'elimine
 593       wiretube = geompy.MakeWire(edges)
 594       #wiretube = edgesTriees[-1]
 595       geomPublish(initLog.debug, wiretube, 'wiretubePlace', self.numeroCas)
 596
 597     # ---------------------------------------------------------
 598
 599     else: # fissure elliptique, longue ou courte
 600       if externe:
 601         raybor = de/2.
 602         dp = -1.0
 603       else:
 604         raybor = de/2. - epais
 605         dp = +1.0
 606       prof = dp * profondeur
 607       cosaz = math.cos(azimut*math.pi/180.)
 608       sinaz = math.sin(azimut*math.pi/180.)
 609       alfrd = alpha*math.pi/180.
 610       pb = geompy.MakeVertex(raybor*cosaz, raybor*sinaz, -l_tube_p1, "pb")
 611       pi = geompy.MakeTranslation(pb, prof*cosaz, prof*sinaz, 0., "pi")
 612       pbv = geompy.MakeTranslation(pb, -profondeur*sinaz, profondeur*cosaz, 0., "pbv")
 613       ayb  = geompy.MakeVector(pb,pbv, "ayb")
 614       pb0 = geompy.MakeRotation(pb, axe, alfrd, "pb0")
 615       ay0 = geompy.MakeRotation(ayb, axe, alfrd, "ay0")
 616       pi0 = geompy.MakeRotation(pi, axe, alfrd, "pi0")
 617       az_ = geompy.MakeVector(pi0, pb0, "az_")
 618       az0 = geompy.MakeTranslationVector(az_, az_, "az0") #normale sortante
 619       ax0 = geompy.MakeRotation(ay0, az0, -math.pi/2.0, "ax0")
 620       ax1 = geompy.MakeRotation(ax0, az0, orientation*math.pi/180., "ax1")
 621       ay1 = geompy.MakeRotation(ay0, az0, orientation*math.pi/180., "ay1")
 622       originLCS = geompy.MakeMarker(0, 0, 0, 1, 0, 0, 0, 1, 0, "originLCS")
 623       coo = geompy.PointCoordinates(pb0)
 624       cox = geompy.VectorCoordinates(ax1)
 625       coy = geompy.VectorCoordinates(ay1)
 626       localLCS = geompy.MakeMarker(coo[0], coo[1], coo[2], cox[0], cox[1], cox[2], coy[0], coy[1], coy[2], "localLCS")
 627
 628       pco = geompy.MakeVertex(0, 0, -profondeur, "pco")
 629       pao = geompy.MakeRotation(pco, OY, 0.6*math.pi, "pao")
 630       pbo = geompy.MakeRotation(pco, OY, -0.6*math.pi, "pbo")
 631       pce = geompy.MakeVertex(0, 0, 0.1*profondeur,"pce")
 632       arcoo = geompy.MakeArc(pao, pco, pbo, "arcoo")
 633       linoo = geompy.MakeArc(pao, pce, pbo, "linoo")
 634       scalex = longueur/profondeur
 635       arco =geompy.MakeScaleAlongAxes(arcoo, O, scalex, 1., 1., "arco")
 636       lino =geompy.MakeScaleAlongAxes(linoo, O, scalex, 1., 1., "lino")
 637       arci = geompy.MakePosition(arco, originLCS, localLCS, "arci")
 638       arce = geompy.MakePosition(lino, originLCS, localLCS, "arce")
 639       facefiss = geompy.MakeFaceWires([arce, arci], 0)
 640       geomPublish(initLog.debug, facefiss, 'facefissPlace', self.numeroCas )
 641       edges = geompy.ExtractShapes(facefiss, geompy.ShapeType["EDGE"], True)
 642       edgesTriees, _, _ = sortEdges(edges)
 643       edgetube = edgesTriees[-1] # la plus grande correspond à arci
 644       wiretube = edgetube
 645
 646       pc = geompy.MakeTranslation(pb, 0.5*prof*cosaz, 0.5*prof*sinaz, 0.)
 647       centre = geompy.MakeRotation(pc, axe, alfrd)
 648       geomPublish(initLog.debug, centre, 'centrefissPlace', self.numeroCas )
 649
 650     coordsNoeudsFissure = genereMeshCalculZoneDefaut(facefiss, profondeur/2. ,profondeur, \
 651                                                      mailleur, self.numeroCas)
 652
 653     return [facefiss, centre, lgInfluence, coordsNoeudsFissure, wiretube, edgetube]
 654
 655   # ---------------------------------------------------------------------------
 656   def setParamMaillageFissure(self):
 657     """
 658     Paramètres du maillage de la fissure pour le tuyau coudé
 659     Voir également setParamShapeFissure, paramètres rayonPipe et lenSegPipe.
 660     nbSegRad = nombre de couronnes
 661     nbSegCercle = nombre de secteurs
 662     areteFaceFissure = taille cible de l'arête des triangles en face de fissure.
 663     """
 664     self.maillageFissureParams = dict(nomRep        = os.curdir,
 665                                       nomFicSain    = self.nomCas,
 666                                       nomFicFissure = 'fissure_' + self.nomCas,
 667                                       nbsegRad      = 5,
 668                                       nbsegCercle   = 6,
 669                                       areteFaceFissure = 5)
 670
 671   # ---------------------------------------------------------------------------
 672   def genereZoneDefaut(self, geometriesSaines, maillagesSains, shapesFissure, shapeFissureParams, maillageFissureParams):
 673     elementsDefaut = creeZoneDefautDansObjetSain(geometriesSaines, maillagesSains, shapesFissure, shapeFissureParams, maillageFissureParams, \
 674                           self.numeroCas)
 675     return elementsDefaut
 676
 677   # ---------------------------------------------------------------------------
 678   def genereMaillageFissure(self, geometriesSaines, maillagesSains, \
 679                             shapesFissure, shapeFissureParams, \
 680                             maillageFissureParams, elementsDefaut, step, \
 681                             mailleur="MeshGems"):
 682     maillageFissure = construitFissureGenerale(shapesFissure, shapeFissureParams, \
 683                                                maillageFissureParams, elementsDefaut, \
 684                                                mailleur, self.numeroCas)
 685     return maillageFissure
 686
 687   # ---------------------------------------------------------------------------
 688   def setReferencesMaillageFissure(self):
 689     self.referencesMaillageFissure = dict( \
 690                                           Entity_Quad_Edge = 975, \
 691                                           Entity_Quad_Quadrangle = 6842, \
 692                                           Entity_Quad_Hexa = 8994, \
 693                                           Entity_Node = 77917, \
 694                                           Entity_Quad_Triangle = 2182, \
 695                                           Entity_Quad_Tetra = 20135, \
 696                                           Entity_Quad_Pyramid = 1038, \
 697                                           Entity_Quad_Penta = 972 \
 698                                          )