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=None):
  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°%s", 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   putName(meshFondExt, "wireFondFiss", i_pref=nro_cas)
  67   algo1d = meshFondExt.Segment()
  68   hypo1d = algo1d.Adaptive(lgmin, lgmax, deflexion) # a ajuster selon la profondeur de la fissure
  69   putName(algo1d.GetSubMesh(), "wireFondFiss", i_pref=nro_cas)
  70   putName(algo1d, "algo1d_wireFondFiss", i_pref=nro_cas)
  71   putName(hypo1d, "hypo1d_wireFondFiss", i_pref=nro_cas)
  72
  73   is_done = meshFondExt.Compute()
  74   text = "calculePointsAxiauxPipe meshFondExt.Compute"
  75   if is_done:
  76     logging.info(text)
  77   else:
  78     text = "Erreur au calcul du maillage.\n" + text
  79     logging.info(text)
  80     raise Exception(text)
  81
  82   ptGSdic = dict() # dictionnaire [paramètre sur la courbe] --> point géométrique
  83   allNodeIds = meshFondExt.GetNodesId()
  84   for nodeId in allNodeIds:
  85     xyz = meshFondExt.GetNodeXYZ(nodeId)
  86     #logging.debug("nodeId %s, coords %s", nodeId, str(xyz))
  87     point = geompy.MakeVertex(xyz[0], xyz[1], xyz[2])
  88     parametre, _, EdgeInWireIndex = geompy.MakeProjectionOnWire(point, wireFondFiss) # parametre compris entre 0 et 1
  89     edgeOrder = edgesIdByOrientation[EdgeInWireIndex]
  90     ptGSdic[(edgeOrder, EdgeInWireIndex, parametre)] = point
  91     #logging.debug("nodeId %s, parametre %s", nodeId, str(parametre))
  92   usort = sorted(ptGSdic)
  93   logging.debug("nombre de points obtenus par deflexion %s",len(usort))
  94
  95   centres = list()
  96   origins = list()
  97   normals = list()
  98   for edu in usort:
  99     vertcx = ptGSdic[edu]
 100     norm = geompy.MakeTangentOnCurve(edgesFondFiss[edu[1]], edu[2])
 101     plan = geompy.MakePlane(vertcx, norm, 3.*rayonPipe)
 102     part = geompy.MakePartition([plan], [wirePipeFiss], list(), list(), geompy.ShapeType["VERTEX"], 0, list(), 0)
 103     liste = geompy.ExtractShapes(part, geompy.ShapeType["VERTEX"], True)
 104     if ( len(liste) == 5 ): # 4 coins du plan plus intersection recherchée
 105       for point in liste:
 106         if geompy.MinDistance(point, vertcx) < 1.1*rayonPipe: # les quatre coins sont plus loin
 107           vertpx = point
 108           break
 109       centres.append(vertcx)
 110       origins.append(vertpx)
 111       normals.append(norm)
 112 #      name = "vertcx%d"%i
 113 #      geompy.addToStudyInFather(wireFondFiss, vertcx, name)
 114 #      name = "vertpx%d"%i
 115 #      geompy.addToStudyInFather(wireFondFiss, vertpx, name)
 116 #      name = "plan%d"%i
 117 #      geompy.addToStudyInFather(wireFondFiss, plan, name)
 118
 119   # --- maillage du pipe étendu, sans tenir compte de l'intersection avec la face de peau
 120
 121   logging.debug("nbsegCercle %s", nbsegCercle)
 122
 123   # -----------------------------------------------------------------------
 124   # --- points géométriques
 125
 126   gptsdisks = list() # vertices géométrie de tous les disques
 127   raydisks = [list() for _ in range(nbsegCercle)]
 128   for indice, centres_i in enumerate(centres): # boucle sur les disques
 129     gptdsk = list() # vertices géométrie d'un disque
 130     vertcx = centres_i
 131     vertpx = origins[indice]
 132     normal = normals[indice]
 133     vec1 = geompy.MakeVector(vertcx, vertpx)
 134
 135     points = [vertcx] # les points du rayon de référence
 136     dist_0 = rayonPipe/float(nbsegRad)
 137     for j_aux in range(nbsegRad):
 138       point = geompy.MakeTranslationVectorDistance(vertcx, vec1, float(j_aux+1)*dist_0)
 139       points.append(point)
 140     gptdsk.append(points)
 141     point = geompy.MakeTranslationVectorDistance(vertcx, vec1, 1.5*rayonPipe)
 142     rayon = geompy.MakeLineTwoPnt(vertcx, point)
 143     raydisks[0].append(rayon)
 144
 145     angle_0 = 2.*math.pi/float(nbsegCercle)
 146     for k_aux in range(nbsegCercle-1):
 147       angle = float(k_aux+1)*angle_0
 148       pts = [vertcx] # les points d'un rayon obtenu par rotation
 149       for j_aux in range(nbsegRad):
 150         point = geompy.MakeRotation(points[j_aux+1], normal, angle)
 151         pts.append(point)
 152       gptdsk.append(pts)
 153       ray = geompy.MakeRotation(rayon, normal, angle)
 154       raydisks[k_aux+1].append(ray)
 155
 156     gptsdisks.append(gptdsk)
 157
 158   return (centres, gptsdisks, raydisks)