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