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