src/Tools/blocFissure/gmu/calculePointsAxiauxPipe.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 """Préparation maillage du pipe"""
  21
  22 import logging
  23 import math
  24
  25 from .geomsmesh import geompy
  26 from .geomsmesh import smesh
  27
  28 from .putName import putName
  29
  30 def calculePointsAxiauxPipe(edgesFondFiss, edgesIdByOrientation, facesDefaut,
  31                             centreFondFiss, wireFondFiss, wirePipeFiss,
  32                             lenSegPipe, rayonPipe, nbsegCercle, nbsegRad, \
  33                             nro_cas=-1):
  34   """Préparation maillage du pipe :
  35
  36   - détections des points a respecter : jonction des edges/faces constituant
  37     la face de fissure externe au pipe
  38   - points sur les edges de fond de fissure et edges pipe/face fissure,
  39   - vecteurs tangents au fond de fissure (normal au disque maillé)
  40   """
  41
  42   logging.info('start')
  43   logging.info("Pour le cas n°%d", nro_cas)
  44
  45   # --- option de maillage selon le rayon de courbure du fond de fissure
  46   lenEdgeFondExt = 0
  47   for edff in edgesFondFiss:
  48     lenEdgeFondExt += geompy.BasicProperties(edff)[0]
  49
  50   disfond = list()
  51   for filling in facesDefaut:
  52     disfond.append(geompy.MinDistance(centreFondFiss, filling))
  53   disfond.sort()
  54   rcourb = disfond[0]
  55   texte = "rcourb: {}, lenEdgeFondExt: {}, lenSegPipe: {}".format(rcourb, lenEdgeFondExt, lenSegPipe)
  56   logging.info(texte)
  57   nbSegQuart = 5 # on veut 5 segments min sur un quart de cercle
  58   alpha = math.pi/(4*nbSegQuart)
  59   deflexion = rcourb*(1.0 -math.cos(alpha))
  60   lgmin = lenSegPipe*0.25
  61   lgmax = lenSegPipe*1.5
  62   texte = "==> deflexion: {}, lgmin: {}, lgmax: {}".format(deflexion, lgmin, lgmax)
  63   logging.info(texte)
  64
  65   meshFondExt = smesh.Mesh(wireFondFiss)
  66   algo1d = meshFondExt.Segment()
  67   hypo1d = algo1d.Adaptive(lgmin, lgmax, deflexion) # a ajuster selon la profondeur de la fissure
  68   putName(algo1d.GetSubMesh(), "wireFondFiss", i_pref=nro_cas)
  69   putName(algo1d, "algo1d_wireFondFiss", i_pref=nro_cas)
  70   putName(hypo1d, "hypo1d_wireFondFiss", i_pref=nro_cas)
  71
  72   is_done = meshFondExt.Compute()
  73   text = "calculePointsAxiauxPipe meshFondExt.Compute"
  74   if is_done:
  75     logging.info(text)
  76   else:
  77     text = "Erreur au calcul du maillage.\n" + text
  78     logging.info(text)
  79     raise Exception(text)
  80
  81   ptGSdic = dict() # dictionnaire [paramètre sur la courbe] --> point géométrique
  82   allNodeIds = meshFondExt.GetNodesId()
  83   for nodeId in allNodeIds:
  84     xyz = meshFondExt.GetNodeXYZ(nodeId)
  85     #logging.debug("nodeId %s, coords %s", nodeId, str(xyz))
  86     point = geompy.MakeVertex(xyz[0], xyz[1], xyz[2])
  87     parametre, _, EdgeInWireIndex = geompy.MakeProjectionOnWire(point, wireFondFiss) # parametre compris entre 0 et 1
  88     edgeOrder = edgesIdByOrientation[EdgeInWireIndex]
  89     ptGSdic[(edgeOrder, EdgeInWireIndex, parametre)] = point
  90     #logging.debug("nodeId %s, parametre %s", nodeId, str(parametre))
  91   usort = sorted(ptGSdic)
  92   logging.debug("nombre de points obtenus par deflexion %s",len(usort))
  93
  94   centres = list()
  95   origins = list()
  96   normals = list()
  97   for edu in usort:
  98     ied = edu[1]
  99     parametre = edu[2]
 100     vertcx = ptGSdic[edu]
 101     norm = geompy.MakeTangentOnCurve(edgesFondFiss[ied], parametre)
 102     plan = geompy.MakePlane(vertcx, norm, 3.*rayonPipe)
 103     part = geompy.MakePartition([plan], [wirePipeFiss], list(), list(), geompy.ShapeType["VERTEX"], 0, list(), 0)
 104     liste = geompy.ExtractShapes(part, geompy.ShapeType["VERTEX"], True)
 105     if len(liste) == 5: # 4 coins du plan plus intersection recherchée
 106       for point in liste:
 107         if geompy.MinDistance(point, vertcx) < 1.1*rayonPipe: # les quatre coins sont plus loin
 108           vertpx = point
 109           break
 110       centres.append(vertcx)
 111       origins.append(vertpx)
 112       normals.append(norm)
 113 #      name = "vertcx%d"%i
 114 #      geompy.addToStudyInFather(wireFondFiss, vertcx, name)
 115 #      name = "vertpx%d"%i
 116 #      geompy.addToStudyInFather(wireFondFiss, vertpx, name)
 117 #      name = "plan%d"%i
 118 #      geompy.addToStudyInFather(wireFondFiss, plan, name)
 119
 120   # --- maillage du pipe étendu, sans tenir compte de l'intersection avec la face de peau
 121
 122   logging.debug("nbsegCercle %s", nbsegCercle)
 123
 124   # -----------------------------------------------------------------------
 125   # --- points géométriques
 126
 127   gptsdisks = list() # vertices géométrie de tous les disques
 128   raydisks = [list() for _ in range(nbsegCercle)]
 129   for indice, centres_i in enumerate(centres): # boucle sur les disques
 130     gptdsk = list() # vertices géométrie d'un disque
 131     vertcx = centres_i
 132     vertpx = origins[indice]
 133     normal = normals[indice]
 134     vec1 = geompy.MakeVector(vertcx, vertpx)
 135
 136     points = [vertcx] # les points du rayon de référence
 137     dist_0 = rayonPipe/float(nbsegRad)
 138     for j in range(nbsegRad):
 139       point = geompy.MakeTranslationVectorDistance(vertcx, vec1, (j+1)*dist_0)
 140       points.append(point)
 141     gptdsk.append(points)
 142     point = geompy.MakeTranslationVectorDistance(vertcx, vec1, 1.5*rayonPipe)
 143     rayon = geompy.MakeLineTwoPnt(vertcx, point)
 144     raydisks[0].append(rayon)
 145
 146     angle_0 = 2.*math.pi/float(nbsegCercle)
 147     for k in range(nbsegCercle-1):
 148       angle = float(k+1)*angle_0
 149       pts = [vertcx] # les points d'un rayon obtenu par rotation
 150       for j in range(nbsegRad):
 151         point = geompy.MakeRotation(points[j+1], normal, angle)
 152         pts.append(point)
 153       gptdsk.append(pts)
 154       ray = geompy.MakeRotation(rayon, normal, angle)
 155       raydisks[k+1].append(ray)
 156
 157     gptsdisks.append(gptdsk)
 158
 159   return (centres, gptsdisks, raydisks)