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