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