src/Tools/blocFissure/gmu/fissureCoude.py

   1 # -*- coding: utf-8 -*-
   2 # Copyright (C) 2014-2020  EDF R&D
   3 #
   4 # This library is free software; you can redistribute it and/or
   5 # modify it under the terms of the GNU Lesser General Public
   6 # License as published by the Free Software Foundation; either
   7 # version 2.1 of the License, or (at your option) any later version.
   8 #
   9 # This library is distributed in the hope that it will be useful,
  10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 # Lesser General Public License for more details.
  13 #
  14 # You should have received a copy of the GNU Lesser General Public
  15 # License along with this library; if not, write to the Free Software
  16 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 #
  18 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 #
  20
  21 import os
  22
  23 from .geomsmesh import geompy, smesh
  24 from .geomsmesh import geomPublish
  25 from .geomsmesh import geomPublishInFather
  26 from . import initLog
  27
  28 import math
  29 import GEOM
  30 import SALOMEDS
  31 import SMESH
  32 #import StdMeshers
  33 #import GHS3DPlugin
  34 #import NETGENPlugin
  35 import logging
  36
  37 from .fissureGenerique import fissureGenerique
  38
  39 from .triedreBase import triedreBase
  40 from .genereMeshCalculZoneDefaut import genereMeshCalculZoneDefaut
  41 from .creeZoneDefautDansObjetSain import creeZoneDefautDansObjetSain
  42 from .construitFissureGenerale import construitFissureGenerale
  43 from .sortEdges import sortEdges
  44
  45 O, OX, OY, OZ = triedreBase()
  46
  47 class fissureCoude(fissureGenerique):
  48   """
  49   problème de fissure du Coude : version de base
  50   maillage hexa
  51   """
  52
  53   nomProbleme = "fissureCoude"
  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     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" )
  96     geomPublish(initLog.debug,  tube_1, "tube_1" )
  97     geomPublish(initLog.debug,  coude, "coude" )
  98     geomPublish(initLog.debug,  tube_2, "tube_2" )
  99
 100     P1 = O
 101     geompy.addToStudy(P1, "P1" )
 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     geompy.addToStudy(P2, "P2" )
 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 )
 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 = []
 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')
 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')
 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')
 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     smesh.SetName(algo3d, "algo3d_maillageSain")
 266     smesh.SetName(algo2d, "algo2d_maillageSain")
 267
 268     algo1d_long_p1 = maillageSain.Segment(geom=long_p1)
 269     hypo1d_long_p1 = algo1d_long_p1.NumberOfSegments(n_long_p1)
 270     smesh.SetName(algo1d_long_p1, "algo1d_long_p1")
 271     smesh.SetName(hypo1d_long_p1, "hypo1d_long_p1")
 272
 273     algo1d_ep = maillageSain.Segment(geom=ep)
 274     hypo1d_ep = algo1d_ep.NumberOfSegments(n_ep)
 275     smesh.SetName(algo1d_ep, "algo1d_ep")
 276     smesh.SetName(hypo1d_ep, "hypo1d_ep")
 277
 278     algo1d_long_coude = maillageSain.Segment(geom=long_coude)
 279     hypo1d_long_coude = algo1d_long_coude.NumberOfSegments(n_long_coude)
 280     smesh.SetName(algo1d_long_coude, "algo1d_long_coude")
 281     smesh.SetName(hypo1d_long_coude, "hypo1d_long_coude")
 282
 283     algo1d_circ_g = maillageSain.Segment(geom=circ_g)
 284     hypo1d_circ_g = algo1d_circ_g.NumberOfSegments(n_circ_g)
 285     smesh.SetName(algo1d_circ_g, "algo1d_circ_g")
 286     smesh.SetName(hypo1d_circ_g, "hypo1d_circ_g")
 287
 288     algo1d_circ_d = maillageSain.Segment(geom=circ_d)
 289     hypo1d_circ_d = algo1d_circ_d.NumberOfSegments(n_circ_d)
 290     smesh.SetName(algo1d_circ_d, "algo1d_circ_d")
 291     smesh.SetName(hypo1d_circ_d, "hypo1d_circ_d")
 292
 293     algo1d_long_p2 = maillageSain.Segment(geom=long_p2)
 294     hypo1d_long_p2 = algo1d_long_p2.NumberOfSegments(n_long_p2)
 295     smesh.SetName(algo1d_long_p2, "algo1d_long_p2")
 296     smesh.SetName(hypo1d_long_p2, "hypo1d_long_p2")
 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     mp1 = maillageSain.GroupOnGeom(P1,'P1',SMESH.NODE)
 308     mp2 = maillageSain.GroupOnGeom(P2,'P2',SMESH.NODE)
 309     ext = maillageSain.GroupOnGeom(EXTUBE,'EXTUBE',SMESH.FACE)
 310     btu = maillageSain.GroupOnGeom(BORDTU,'BORDTU',SMESH.EDGE)
 311     clg = maillageSain.GroupOnGeom(CLGV,'CLGV',SMESH.FACE)
 312     pei = maillageSain.GroupOnGeom(PEAUINT,'PEAUINT',SMESH.FACE)
 313     pex = maillageSain.GroupOnGeom(PEAUEXT,'PEAUEXT',SMESH.FACE)
 314     cou = 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     logging.info("genereShapeFissure %s", self.nomCas)
 349     logging.info("shapeFissureParams %s", shapeFissureParams)
 350
 351     angleCoude = geomParams['angleCoude']
 352     r_cintr    = geomParams['r_cintr']
 353     l_tube_p1  = geomParams['l_tube_p1']
 354     l_tube_p2  = geomParams['l_tube_p2']
 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
 370
 371     azimut = -azimut # axe inverse / ASCOUF
 372     axe = geompy.MakeTranslation(OY, -r_cintr, 0, -l_tube_p1)
 373     geomPublish(initLog.debug, axe,"axe")
 374
 375     if not lgInfluence:
 376       lgInfluence = profondeur
 377
 378     if longueur > 2*profondeur:
 379       self.fissureLongue=True
 380     else:
 381       self.fissureLongue=False
 382       self.elliptique = True
 383
 384     self.circonferentielle = False
 385     self.longitudinale = False
 386     if self.fissureLongue and not self.elliptique:
 387       if abs(orientation) < 45 :
 388         self.longitudinale = True
 389       else:
 390         self.circonferentielle = True
 391
 392     nbp1 = 10
 393     if self.circonferentielle:
 394       if externe:
 395         dp = -1.0
 396         raybor = de/2.
 397         rayint = raybor - profondeur
 398         rayext = raybor + profondeur/5.0
 399       else:
 400         dp = 1.0
 401         raybor = de/2. - epais
 402         rayint = raybor + profondeur
 403         rayext = raybor - profondeur/5.0
 404       lgfond = longueur -2*profondeur
 405       angle = lgfond/(2*raybor)
 406       pb = geompy.MakeVertex(raybor, 0, 0)
 407       pi = geompy.MakeVertex(rayint, 0, 0)
 408       pbl = geompy.MakeRotation(pb, OZ, angle)
 409       pbr = geompy.MakeRotation(pb, OZ, -angle)
 410       geomPublish(initLog.debug, pbl,"pbl")
 411       geomPublish(initLog.debug, pbr,"pbr")
 412       pal = geompy.MakeTranslationVector(pbl, OZ)
 413       par = geompy.MakeTranslationVector(pbr, OZ)
 414       axl = geompy.MakeVector(pbl,pal)
 415       axr = geompy.MakeVector(pbr,par)
 416       pil = geompy.MakeRotation(pi, OZ, angle)
 417       pir = geompy.MakeRotation(pi, OZ, -angle)
 418       points = []
 419       nbp = 3*nbp1
 420       for i in range(nbp):
 421         angi = dp*(nbp -i)*(2.0*math.pi/3.0)/nbp
 422         pt = geompy.MakeRotation(pil, axl, angi)
 423         points.append(pt)
 424       for i in range(nbp):
 425         angi = angle -2.0*i*angle/nbp
 426         pt = geompy.MakeRotation(pi, OZ, angi)
 427         points.append(pt)
 428       for i in range(nbp+1):
 429         angi = -dp*i*(2.0*math.pi/3.0)/nbp
 430         pt = geompy.MakeRotation(pir, axr, angi)
 431         points.append(pt)
 432       for i, pt in enumerate(points):
 433         pt = geompy.MakeRotation(pt, OZ, azimut*math.pi/180.)
 434         pt = geompy.MakeTranslation(pt, 0, 0, -l_tube_p1)
 435         pt = geompy.MakeRotation(pt, axe, alpha*math.pi/180.)
 436         points[i] = pt
 437       wire0 = geompy.MakeInterpol(points[0:nbp+1])
 438       wire1 = geompy.MakeInterpol(points[nbp:2*nbp+1])
 439       wire2 = geompy.MakeInterpol(points[2*nbp:3*nbp+1])
 440       #wiretube = geompy.MakeInterpol(points)
 441       wiretube=geompy.MakeWire([wire0,wire1,wire2])
 442       geomPublish(initLog.debug, wiretube,"wiretube")
 443
 444       pe = geompy.MakeVertex(rayext, 0, 0)
 445       pe = geompy.MakeRotation(pe, OZ, azimut*math.pi/180.)
 446       pe = geompy.MakeTranslation(pe, 0, 0, -l_tube_p1)
 447       pe = geompy.MakeRotation(pe, axe, alpha*math.pi/180.)
 448
 449       arce = geompy.MakeArc(points[0], pe, points[-1])
 450       geomPublish(initLog.debug, arce,"arce")
 451
 452       facefiss = geompy.MakeFaceWires([arce, wiretube], 1)
 453       geomPublish(initLog.debug,  facefiss, 'facefissPlace' )
 454
 455       pc = geompy.MakeVertex((raybor + rayint)/2.0, 0, 0)
 456       centre = geompy.MakeRotation(pc, OZ, azimut*math.pi/180.)
 457       centre = geompy.MakeTranslation(centre, 0, 0, -l_tube_p1)
 458       centre = geompy.MakeRotation(centre, axe, alpha*math.pi/180.)
 459       geomPublish(initLog.debug,  centre, 'centrefissPlace' )
 460
 461       wiretube = geompy.GetInPlace(facefiss, wiretube)
 462       geomPublish(initLog.debug, wiretube, 'wiretubePlace' )
 463       try:
 464         edgetube = geompy.MakeEdgeWire(wiretube)
 465         geomPublish(initLog.debug, edgetube,"edgetube")
 466       except:
 467         logging.debug("erreur MakeEdgeWire sur fond de fissure, on fait sans")
 468         edgetube = None
 469
 470     # ---------------------------------------------------------
 471
 472     elif self.longitudinale:
 473       if externe:
 474         raybor = de/2.
 475         dp = -1.0
 476       else:
 477         raybor = de/2. - epais
 478         dp = +1.0
 479       prof = dp * profondeur
 480       lgfond = longueur -2*profondeur
 481       cosaz = math.cos(azimut*math.pi/180.)
 482       sinaz = math.sin(azimut*math.pi/180.)
 483       alfrd = alpha*math.pi/180.
 484       rayxy = r_cintr + raybor*cosaz
 485       angle = lgfond/(2.*rayxy)
 486       logging.debug("longueur: %s, angle: %s, rayon: %s",lgfond, angle, rayxy)
 487       pb = geompy.MakeVertex(raybor*cosaz, raybor*sinaz, -l_tube_p1, "pb")
 488       pi = geompy.MakeTranslation(pb, prof*cosaz, prof*sinaz, 0., "pi")
 489       pbv = geompy.MakeTranslation(pb, -sinaz, cosaz, 0., "pbv")
 490       axb  = geompy.MakeVector(pb,pbv, "axb")
 491       pbl = geompy.MakeRotation(pb, axe, alfrd -angle, "pbl")
 492       pbr = geompy.MakeRotation(pb, axe, alfrd +angle, "pbr")
 493       axl = geompy.MakeRotation(axb, axe, alfrd -angle, "axl")
 494       axr = geompy.MakeRotation(axb, axe, alfrd +angle, "axr")
 495       pil = geompy.MakeRotation(pi, axe, alfrd -angle, "pil")
 496       pir = geompy.MakeRotation(pi, axe, alfrd +angle, "pir")
 497
 498       curves = []
 499
 500       points = []
 501       nbp = 3*nbp1
 502       xs = []
 503       totx = 0
 504       for i in range(nbp+2):
 505         x = math.sin(i*math.pi/(nbp+1)) # fonction de répartition des points : distance relative
 506         x2 = x*x
 507         totx += x2
 508         xs.append(totx)
 509         logging.debug("x2: %s, totx: %s", x2, totx)
 510       for i in range(nbp+1):
 511         #posi = nbp -i             # répartition équidistante des points sur la courbe
 512         posi = nbp*(1 -xs[i]/totx) # points plus resserrés aux extrémités de la courbe
 513         angi = -dp*posi*(5.0*math.pi/8.0)/nbp
 514         pt = geompy.MakeRotation(pil, axl, angi)
 515         points.append(pt)
 516       curves.append(geompy.MakeInterpol(points))
 517       point0 = points[0]
 518       geomPublish(initLog.debug, curves[-1],"curve0")
 519 #      for i, pt in enumerate(points):
 520 #        name = "point%d"%i
 521 #        geomPublishInFather(initLog.debug,curves[-1], pt, name)
 522
 523       points = []
 524       nbp = 3*nbp1
 525       xs =[]
 526       totx = 0
 527       for i in range(nbp+1):
 528         x = math.sin(i*math.pi/nbp)
 529         #x = 1.0 # répartition équidistante des points sur la courbe
 530         x2 = x*x # points plus resserrés aux extrémités de la courbe
 531         totx += x2
 532         xs.append(totx)
 533         logging.debug("x2: %s, totx: %s", x2, totx)
 534       for i in range(nbp):
 535         angi = alfrd -angle +2.0*angle*xs[i]/totx
 536         pt = geompy.MakeRotation(pi, axe, angi)
 537         points.append(pt)
 538       curves.append(geompy.MakeInterpol(points))
 539       geomPublish(initLog.debug, curves[-1],"curve1")
 540 #      for i, pt in enumerate(points):
 541 #        name = "point%d"%i
 542 #        geomPublishInFather(initLog.debug,curves[-1], pt, name)
 543
 544       points = []
 545       nbp = 3*nbp1
 546       xs = []
 547       totx = 0
 548       for i in range(nbp+2):
 549         x = math.sin(i*math.pi/(nbp+1))
 550         x2 = x*x
 551         totx += x2
 552         xs.append(totx)
 553         logging.debug("x2: %s, totx: %s", x2, totx)
 554       for i in range(nbp+1):
 555         #posi = nbp -i        # répartition équidistante des points sur la courbe
 556         posi = nbp*xs[i]/totx # points plus resserrés aux extrémités de la courbe
 557         angi = dp*posi*(5.0*math.pi/8.0)/nbp
 558         pt = geompy.MakeRotation(pir, axr, angi)
 559         points.append(pt)
 560       curves.append(geompy.MakeInterpol(points))
 561       point1 = points[-1]
 562       geomPublish(initLog.debug, curves[-1],"curve2")
 563 #      for i, pt in enumerate(points):
 564 #        name = "point%d"%i
 565 #        geomPublishInFather(initLog.debug,curves[-1], pt, name)
 566
 567       wiretube = geompy.MakeWire(curves)
 568       geomPublish(initLog.debug, wiretube,"wiretube")
 569       try:
 570         edgetube = geompy.MakeEdgeWire(wiretube)
 571         geomPublish(initLog.debug, edgetube,"edgetube")
 572       except:
 573         logging.debug("erreur MakeEdgeWire sur fond de fissure, on fait sans")
 574         edgetube = None
 575
 576       pts = []
 577       pts.append(point0)
 578       dpr = prof*math.cos(5.0*math.pi/8.0)
 579       pe = geompy.MakeTranslation(pb, dpr*cosaz, dpr*sinaz, 0., "pe")
 580       for i in range(nbp):
 581         angi = alfrd -angle +2.0*i*angle/nbp
 582         pt = geompy.MakeRotation(pe, axe, angi)
 583         pts.append(pt)
 584       pts.append(point1)
 585       arce = geompy.MakeInterpol(pts)
 586       geomPublish(initLog.debug, arce,"arce")
 587
 588       facefiss = geompy.MakeFaceWires([arce, wiretube], 0)
 589       geomPublish(initLog.debug,  facefiss, 'facefissPlace' )
 590
 591       pc = geompy.MakeTranslation(pb, 0.5*prof*cosaz, 0.5*prof*sinaz, 0.)
 592       centre = geompy.MakeRotation(pc, axe, alfrd)
 593       geomPublish(initLog.debug,  centre, 'centrefissPlace' )
 594
 595       edges = geompy.ExtractShapes(facefiss, geompy.ShapeType["EDGE"], True)
 596       edgesTriees, minl, maxl = sortEdges(edges)
 597       edges = edgesTriees[:-1] # la plus grande correspond à arce, on l'elimine
 598       wiretube = geompy.MakeWire(edges)
 599       #wiretube = edgesTriees[-1]
 600       geomPublish(initLog.debug, wiretube, 'wiretubePlace' )
 601
 602     # ---------------------------------------------------------
 603
 604     else: # fissure elliptique, longue ou courte
 605       if externe:
 606         raybor = de/2.
 607         dp = -1.0
 608       else:
 609         raybor = de/2. - epais
 610         dp = +1.0
 611       prof = dp * profondeur
 612       cosaz = math.cos(azimut*math.pi/180.)
 613       sinaz = math.sin(azimut*math.pi/180.)
 614       alfrd = alpha*math.pi/180.
 615       pb = geompy.MakeVertex(raybor*cosaz, raybor*sinaz, -l_tube_p1, "pb")
 616       pi = geompy.MakeTranslation(pb, prof*cosaz, prof*sinaz, 0., "pi")
 617       pbv = geompy.MakeTranslation(pb, -profondeur*sinaz, profondeur*cosaz, 0., "pbv")
 618       ayb  = geompy.MakeVector(pb,pbv, "ayb")
 619       pb0 = geompy.MakeRotation(pb, axe, alfrd, "pb0")
 620       ay0 = geompy.MakeRotation(ayb, axe, alfrd, "ay0")
 621       pi0 = geompy.MakeRotation(pi, axe, alfrd, "pi0")
 622       az_ = geompy.MakeVector(pi0, pb0, "az_")
 623       az0 = geompy.MakeTranslationVector(az_, az_, "az0") #normale sortante
 624       ax0 = geompy.MakeRotation(ay0, az0, -math.pi/2.0, "ax0")
 625       ax1 = geompy.MakeRotation(ax0, az0, orientation*math.pi/180., "ax1")
 626       ay1 = geompy.MakeRotation(ay0, az0, orientation*math.pi/180., "ay1")
 627       originLCS = geompy.MakeMarker(0, 0, 0, 1, 0, 0, 0, 1, 0, "originLCS")
 628       coo = geompy.PointCoordinates(pb0)
 629       cox = geompy.VectorCoordinates(ax1)
 630       coy = geompy.VectorCoordinates(ay1)
 631       localLCS = geompy.MakeMarker(coo[0], coo[1], coo[2], cox[0], cox[1], cox[2], coy[0], coy[1], coy[2], "localLCS")
 632
 633       pco = geompy.MakeVertex(0, 0, -profondeur, "pco")
 634       pao = geompy.MakeRotation(pco, OY, 0.6*math.pi, "pao")
 635       pbo = geompy.MakeRotation(pco, OY, -0.6*math.pi, "pbo")
 636       pce = geompy.MakeVertex(0, 0, 0.1*profondeur,"pce")
 637       arcoo = geompy.MakeArc(pao, pco, pbo, "arcoo")
 638       linoo = geompy.MakeArc(pao, pce, pbo, "linoo")
 639       scalex = longueur/profondeur
 640       arco =geompy.MakeScaleAlongAxes(arcoo, O, scalex, 1., 1., "arco")
 641       lino =geompy.MakeScaleAlongAxes(linoo, O, scalex, 1., 1., "lino")
 642       arci = geompy.MakePosition(arco, originLCS, localLCS, "arci")
 643       arce = geompy.MakePosition(lino, originLCS, localLCS, "arce")
 644       facefiss = geompy.MakeFaceWires([arce, arci], 0)
 645       geomPublish(initLog.debug,  facefiss, 'facefissPlace' )
 646       edges = geompy.ExtractShapes(facefiss, geompy.ShapeType["EDGE"], True)
 647       edgesTriees, minl, maxl = sortEdges(edges)
 648       edgetube = edgesTriees[-1] # la plus grande correspond à arci
 649       wiretube = edgetube
 650
 651       pc = geompy.MakeTranslation(pb, 0.5*prof*cosaz, 0.5*prof*sinaz, 0.)
 652       centre = geompy.MakeRotation(pc, axe, alfrd)
 653       geomPublish(initLog.debug,  centre, 'centrefissPlace' )
 654
 655     coordsNoeudsFissure = genereMeshCalculZoneDefaut(facefiss, profondeur/2. ,profondeur)
 656
 657     return [facefiss, centre, lgInfluence, coordsNoeudsFissure, wiretube, edgetube]
 658
 659   # ---------------------------------------------------------------------------
 660   def setParamMaillageFissure(self):
 661     """
 662     Paramètres du maillage de la fissure pour le tuyau coudé
 663     Voir également setParamShapeFissure, paramètres rayonPipe et lenSegPipe.
 664     nbSegRad = nombre de couronnes
 665     nbSegCercle = nombre de secteurs
 666     areteFaceFissure = taille cible de l'arête des triangles en face de fissure.
 667     """
 668     self.maillageFissureParams = dict(nomRep        = os.curdir,
 669                                       nomFicSain    = self.nomCas,
 670                                       nomFicFissure = 'fissure_' + self.nomCas,
 671                                       nbsegRad      = 5,
 672                                       nbsegCercle   = 6,
 673                                       areteFaceFissure = 5)
 674
 675   # ---------------------------------------------------------------------------
 676   def genereZoneDefaut(self, geometriesSaines, maillagesSains, shapesFissure, shapeFissureParams, maillageFissureParams):
 677     elementsDefaut = creeZoneDefautDansObjetSain(geometriesSaines, maillagesSains, shapesFissure, shapeFissureParams, maillageFissureParams)
 678     return elementsDefaut
 679
 680   # ---------------------------------------------------------------------------
 681   def genereMaillageFissure(self, geometriesSaines, maillagesSains,
 682                             shapesFissure, shapeFissureParams,
 683                             maillageFissureParams, elementsDefaut, step):
 684     maillageFissure = construitFissureGenerale(maillagesSains,
 685                                                shapesFissure, shapeFissureParams,
 686                                                maillageFissureParams, elementsDefaut, step)
 687     return maillageFissure
 688
 689   # ---------------------------------------------------------------------------
 690   def setReferencesMaillageFissure(self):
 691     self.referencesMaillageFissure = dict(Entity_Node            = 77917,
 692                                           Entity_Quad_Edge       = 975,
 693                                           Entity_Quad_Triangle   = 2182,
 694                                           Entity_Quad_Quadrangle = 6842,
 695                                           Entity_Quad_Tetra      = 20135,
 696                                           Entity_Quad_Hexa       = 8994,
 697                                           Entity_Quad_Penta      = 972,
 698                                           Entity_Quad_Pyramid    = 1038)
 699