Salome HOME
Corrections cas-tests 14 et 27
[modules/smesh.git] / src / Tools / blocFissure / gmu / calculePointsAxiauxPipe.py
index c3f4d5ff29f1fc6e1f07f268c76e4ea24165af64..f67262afca3cc67f6e56c7597f092bdc90efebc4 100644 (file)
 # -*- coding: utf-8 -*-
+# Copyright (C) 2014-2021  EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+"""Préparation maillage du pipe"""
 
 import logging
 import math
 
 from .geomsmesh import geompy
 from .geomsmesh import smesh
-  
-def calculePointsAxiauxPipe(edgesFondFiss, edgesIdByOrientation, facesDefaut, 
+
+from .putName import putName
+
+def calculePointsAxiauxPipe(edgesFondFiss, edgesIdByOrientation, facesDefaut,
                             centreFondFiss, wireFondFiss, wirePipeFiss,
-                            lenSegPipe, rayonPipe, nbsegCercle, nbsegRad):
-  """
-  preparation maillage du pipe :
+                            lenSegPipe, rayonPipe, nbsegCercle, nbsegRad, \
+                            nro_cas=None):
+  """Préparation maillage du pipe :
+
   - détections des points a respecter : jonction des edges/faces constituant
     la face de fissure externe au pipe
   - points sur les edges de fond de fissure et edges pipe/face fissure,
-  - vecteurs tangents au fond de fissure (normal au disque maillé)  
+  - vecteurs tangents au fond de fissure (normal au disque maillé)
   """
-  
+
   logging.info('start')
+  logging.info("Pour le cas n°%s", nro_cas)
 
-  # --- option de maillage selon le rayon de courbure du fond de fissure 
+  # --- option de maillage selon le rayon de courbure du fond de fissure
   lenEdgeFondExt = 0
   for edff in edgesFondFiss:
     lenEdgeFondExt += geompy.BasicProperties(edff)[0]
-  
-  disfond = []
+
+  disfond = list()
   for filling in facesDefaut:
     disfond.append(geompy.MinDistance(centreFondFiss, filling))
   disfond.sort()
   rcourb = disfond[0]
+  texte = "rcourb: {}, lenEdgeFondExt: {}, lenSegPipe: {}".format(rcourb, lenEdgeFondExt, lenSegPipe)
+  logging.info(texte)
   nbSegQuart = 5 # on veut 5 segments min sur un quart de cercle
   alpha = math.pi/(4*nbSegQuart)
   deflexion = rcourb*(1.0 -math.cos(alpha))
   lgmin = lenSegPipe*0.25
-  lgmax = lenSegPipe*1.5               
-  logging.debug("rcourb: %s, lenFond:%s, deflexion: %s, lgmin: %s, lgmax: %s", rcourb, lenEdgeFondExt, deflexion, lgmin, lgmax)  
+  lgmax = lenSegPipe*1.5
+  texte = "==> deflexion: {}, lgmin: {}, lgmax: {}".format(deflexion, lgmin, lgmax)
+  logging.info(texte)
 
   meshFondExt = smesh.Mesh(wireFondFiss)
+  putName(meshFondExt, "wireFondFiss", i_pref=nro_cas)
   algo1d = meshFondExt.Segment()
   hypo1d = algo1d.Adaptive(lgmin, lgmax, deflexion) # a ajuster selon la profondeur de la fissure
-  isDone = meshFondExt.Compute()
-  
-  ptGSdic = {} # dictionnaire [paramètre sur la courbe] --> point géométrique
+  putName(algo1d.GetSubMesh(), "wireFondFiss", i_pref=nro_cas)
+  putName(algo1d, "algo1d_wireFondFiss", i_pref=nro_cas)
+  putName(hypo1d, "hypo1d_wireFondFiss", i_pref=nro_cas)
+
+  is_done = meshFondExt.Compute()
+  text = "calculePointsAxiauxPipe meshFondExt.Compute"
+  if is_done:
+    logging.info(text)
+  else:
+    text = "Erreur au calcul du maillage.\n" + text
+    logging.info(text)
+    raise Exception(text)
+
+  ptGSdic = dict() # dictionnaire [paramètre sur la courbe] --> point géométrique
   allNodeIds = meshFondExt.GetNodesId()
   for nodeId in allNodeIds:
     xyz = meshFondExt.GetNodeXYZ(nodeId)
     #logging.debug("nodeId %s, coords %s", nodeId, str(xyz))
-    pt = geompy.MakeVertex(xyz[0], xyz[1], xyz[2])
-    u, PointOnEdge, EdgeInWireIndex = geompy.MakeProjectionOnWire(pt, wireFondFiss) # u compris entre 0 et 1
+    point = geompy.MakeVertex(xyz[0], xyz[1], xyz[2])
+    parametre, _, EdgeInWireIndex = geompy.MakeProjectionOnWire(point, wireFondFiss) # parametre compris entre 0 et 1
     edgeOrder = edgesIdByOrientation[EdgeInWireIndex]
-    ptGSdic[(edgeOrder, EdgeInWireIndex, u)] = pt
-    #logging.debug("nodeId %s, u %s", nodeId, str(u))
-  usort = sorted(ptGSdic)  
+    ptGSdic[(edgeOrder, EdgeInWireIndex, parametre)] = point
+    #logging.debug("nodeId %s, parametre %s", nodeId, str(parametre))
+  usort = sorted(ptGSdic)
   logging.debug("nombre de points obtenus par deflexion %s",len(usort))
-     
-  centres = []
-  origins = []
-  normals = []      
+
+  centres = list()
+  origins = list()
+  normals = list()
   for edu in usort:
-    ied = edu[1]
-    u = edu[2]
     vertcx = ptGSdic[edu]
-    norm = geompy.MakeTangentOnCurve(edgesFondFiss[ied], u)
-    plan = geompy.MakePlane(vertcx, norm, 3*rayonPipe)
-    part = geompy.MakePartition([plan], [wirePipeFiss], [], [], geompy.ShapeType["VERTEX"], 0, [], 0)
+    norm = geompy.MakeTangentOnCurve(edgesFondFiss[edu[1]], edu[2])
+    plan = geompy.MakePlane(vertcx, norm, 3.*rayonPipe)
+    part = geompy.MakePartition([plan], [wirePipeFiss], list(), list(), geompy.ShapeType["VERTEX"], 0, list(), 0)
     liste = geompy.ExtractShapes(part, geompy.ShapeType["VERTEX"], True)
-    if len(liste) == 5: # 4 coins du plan plus intersection recherchée
+    if ( len(liste) == 5 ): # 4 coins du plan plus intersection recherchée
       for point in liste:
         if geompy.MinDistance(point, vertcx) < 1.1*rayonPipe: # les quatre coins sont plus loin
           vertpx = point
@@ -81,40 +117,42 @@ def calculePointsAxiauxPipe(edgesFondFiss, edgesIdByOrientation, facesDefaut,
 #      geompy.addToStudyInFather(wireFondFiss, plan, name)
 
   # --- maillage du pipe étendu, sans tenir compte de l'intersection avec la face de peau
-      
+
   logging.debug("nbsegCercle %s", nbsegCercle)
-  
+
   # -----------------------------------------------------------------------
   # --- points géométriques
-  
-  gptsdisks = [] # vertices géométrie de tous les disques
-  raydisks = [[] for i in range(nbsegCercle)]
-  for i in range(len(centres)): # boucle sur les disques
-    gptdsk = [] # vertices géométrie d'un disque
-    vertcx = centres[i]
-    vertpx = origins[i]
-    normal = normals[i]
+
+  gptsdisks = list() # vertices géométrie de tous les disques
+  raydisks = [list() for _ in range(nbsegCercle)]
+  for indice, centres_i in enumerate(centres): # boucle sur les disques
+    gptdsk = list() # vertices géométrie d'un disque
+    vertcx = centres_i
+    vertpx = origins[indice]
+    normal = normals[indice]
     vec1 = geompy.MakeVector(vertcx, vertpx)
-    
+
     points = [vertcx] # les points du rayon de référence
-    for j in range(nbsegRad):
-      pt = geompy.MakeTranslationVectorDistance(vertcx, vec1, (j+1)*float(rayonPipe)/nbsegRad)
-      points.append(pt)
+    dist_0 = rayonPipe/float(nbsegRad)
+    for j_aux in range(nbsegRad):
+      point = geompy.MakeTranslationVectorDistance(vertcx, vec1, float(j_aux+1)*dist_0)
+      points.append(point)
     gptdsk.append(points)
-    pt = geompy.MakeTranslationVectorDistance(vertcx, vec1, 1.5*rayonPipe)
-    rayon = geompy.MakeLineTwoPnt(vertcx, pt)
+    point = geompy.MakeTranslationVectorDistance(vertcx, vec1, 1.5*rayonPipe)
+    rayon = geompy.MakeLineTwoPnt(vertcx, point)
     raydisks[0].append(rayon)
-    
-    for k in range(nbsegCercle-1):
-      angle = (k+1)*2*math.pi/nbsegCercle
+
+    angle_0 = 2.*math.pi/float(nbsegCercle)
+    for k_aux in range(nbsegCercle-1):
+      angle = float(k_aux+1)*angle_0
       pts = [vertcx] # les points d'un rayon obtenu par rotation
-      for j in range(nbsegRad):
-        pt = geompy.MakeRotation(points[j+1], normal, angle)
-        pts.append(pt)
+      for j_aux in range(nbsegRad):
+        point = geompy.MakeRotation(points[j_aux+1], normal, angle)
+        pts.append(point)
       gptdsk.append(pts)
       ray = geompy.MakeRotation(rayon, normal, angle)
-      raydisks[k+1].append(ray)
-      
+      raydisks[k_aux+1].append(ray)
+
     gptsdisks.append(gptdsk)
-    
+
   return (centres, gptsdisks, raydisks)