test.models/coronavirus.py

   1 #!/usr/bin/env python
   2 #=============================================================================
   3 # This script was written in march and april 2020 during the beginning in France
   4 # of the coronavirus pandemia. It's generating a parametric geometric model of the virus
   5 # in Shaper.
   6 # Author : Raphaël MARC, EDF R&D France, with the precious help of Artem ZHIDKOV from OpenCascade company.
   7 #=============================================================================
   8 import numpy as np
   9 from salome.shaper import model
  10 from SketchAPI import *
  11 from GeomAPI import *
  12
  13 model.begin()
  14 partSet = model.moduleDocument()
  15 Part_1 = model.addPart(partSet)
  16 Part_1_doc = Part_1.document()
  17 #=============================================================================
  18 # Paramètres du modèle SHAPER :
  19 #=============================================================================
  20 Rext = model.addParameter(Part_1_doc, "Rext", "10")
  21 #print(Rext.value())
  22 model.addParameter(Part_1_doc, "eps", "0.5")
  23 Rint = model.addParameter(Part_1_doc, "Rint", "Rext-eps")
  24 model.addParameter(Part_1_doc, "marge", "0.3")
  25 Rint_m = model.addParameter(Part_1_doc, "Rint_m", "Rint-marge")
  26 longueur = model.addParameter(Part_1_doc, "long_tube", "4")
  27 Rtube = model.addParameter(Part_1_doc, "Rtube", "0.8")
  28 ep_tube = model.addParameter(Part_1_doc, "ep_tube", "0.3")
  29 # ratio diamètre cercle milieu tube / diamètre cercle de base (bottom) :
  30 coef1 = model.addParameter(Part_1_doc, "coef1", "0.5")
  31 # ratio diamètre cercle haut (top) du tube / diamètre cercle de base (bottom) :
  32 coef2 = model.addParameter(Part_1_doc, "coef2", "1.2")
  33 fillet_radius_top = model.addParameter(Part_1_doc, "fillet_radius_top", "0.12")
  34 fillet_radius_bottom = model.addParameter(Part_1_doc, "fillet_radius_bottom", "0.8")
  35 # méthodes dispo sur les paramètres : .value() et .setValue()
  36 bruit = 0.012 # bruit dans la position des tubes (qui est régulière si pas de bruit)
  37
  38 #=============================================================================
  39 # Début réel de la création du modèle géométrique SHAPER :
  40 #=============================================================================
  41 Sphere_ext = model.addSphere(Part_1_doc, model.selection("VERTEX", "PartSet/Origin"), "Rext")
  42 Sphere_ext.setName("Sphere_ext")
  43 # Création de num_pts points uniformément répartis sur la sphère, qui serviront à positionner les tubes :
  44 num_pts = 50
  45 indices = np.arange(0, num_pts, dtype=float) + 0.5
  46 phi0 = np.arccos(1 - 2*indices/num_pts)
  47 theta0 = np.pi * (1 + 5**0.5) * indices
  48 # ajout de bruit Gaussien pour rendre légèrement aléatoire les positions des tubes :
  49 bruit=np.pi*np.random.normal(0, bruit, num_pts)  # std = np.pi*0.05
  50 phi = phi0+bruit
  51 theta = theta0+bruit
  52 #print(max(np.abs((phi-phi0)/phi0)))
  53 #print(max(np.abs((theta-theta0)/theta0)))
  54
  55 listeX, listeY, listeZ = Rint_m.value()*np.cos(theta) * np.sin(phi), \
  56                         Rint_m.value()*np.sin(theta) * np.sin(phi),  \
  57                         Rint_m.value()*np.cos(phi);
  58 x,y,z = (listeX[0], listeY[0], listeZ[0])
  59 #=============================================================================
  60 # On construit le premier tube que l'on copiera puis rotations ensuite :
  61 # 3 sketchs successifs composés de 2 cercles // puis un filling (remplissage)
  62 #=============================================================================
  63 Point_init = model.addPoint(Part_1_doc, x,y,z)
  64 Axis = model.addAxis(Part_1_doc, model.selection("VERTEX", "PartSet/Origin"), Point_init.result())
  65 Plane_1 = model.addPlane(Part_1_doc, Axis.result(), Point_init.result(), True)
  66 Sketch_1 = model.addSketch(Part_1_doc, Plane_1.result())
  67 SketchProjection_1 = Sketch_1.addProjection(Point_init.result(), False)
  68 SketchPoint_1 = SketchProjection_1.createdFeature()
  69 SketchCircle_1_int = Sketch_1.addCircle(0, 0, 1.1)
  70 Sketch_1.setCoincident(SketchAPI_Point(SketchPoint_1).coordinates(), SketchCircle_1_int.center())
  71 SketchCircle_1_ext = Sketch_1.addCircle(0, 0, 1.)
  72 Sketch_1.setCoincident(SketchAPI_Point(SketchPoint_1).coordinates(), SketchCircle_1_ext.center())
  73 Sketch_1.setRadius(SketchCircle_1_int.results()[1], "Rtube")
  74 Sketch_1.setRadius(SketchCircle_1_ext.results()[1], "Rtube+ep_tube")
  75 model.do()
  76
  77 Plane_2 = model.addPlane(Part_1_doc, Plane_1.result(), "long_tube", False)
  78 Sketch_2 = model.addSketch(Part_1_doc, Plane_2.result())
  79 SketchProjection_2 = Sketch_2.addProjection(Point_init.result(), False)
  80 SketchPoint_2 = SketchProjection_2.createdFeature()
  81 SketchCircle_2_int = Sketch_2.addCircle(0, 0, 0.4)
  82 Sketch_2.setCoincident(SketchAPI_Point(SketchPoint_2).coordinates(), SketchCircle_2_int.center())
  83 SketchCircle_2_ext = Sketch_2.addCircle(0, 0, 0.7)
  84 Sketch_2.setCoincident(SketchAPI_Point(SketchPoint_2).coordinates(), SketchCircle_2_ext.center())
  85 Sketch_2.setRadius(SketchCircle_2_int.results()[1], "Rtube*coef1")
  86 Sketch_2.setRadius(SketchCircle_2_ext.results()[1], "(Rtube+ep_tube)*coef1")
  87 model.do()
  88
  89 Plane_3 = model.addPlane(Part_1_doc, Plane_2.result(), "long_tube", False)
  90 Sketch_3 = model.addSketch(Part_1_doc, Plane_3.result())
  91 SketchProjection_4 = Sketch_3.addProjection(Point_init.result(), False)
  92 SketchPoint_3 = SketchProjection_4.createdFeature()
  93 SketchCircle_3_int = Sketch_3.addCircle(0, 0, 0.8)
  94 SketchCircle_3_ext = Sketch_3.addCircle(0, 0, 1.1)
  95 Sketch_3.setRadius(SketchCircle_3_int.results()[1], "Rtube*coef2")
  96 Sketch_3.setRadius(SketchCircle_3_ext.results()[1], "(Rtube+ep_tube)*coef2")
  97 Sketch_3.setCoincident(SketchCircle_3_ext.center(), SketchAPI_Point(SketchPoint_3).coordinates())
  98 Sketch_3.setCoincident(SketchCircle_3_int.center(), SketchAPI_Point(SketchPoint_3).coordinates())
  99 model.do()
 100
 101 Filling_1_objects = \
 102     [model.selection(Sketch_1.defaultResult(), SketchCircle_1_int.defaultResult().shape()), \
 103      model.selection(Sketch_2.defaultResult(), SketchCircle_2_int.defaultResult().shape()), \
 104      model.selection(Sketch_3.defaultResult(), SketchCircle_3_int.defaultResult().shape())]
 105 Filling_1 = model.addFilling(Part_1_doc, Filling_1_objects, "curve_info", 2, 5, 0, 0.0001, 0.0001, False)
 106 Copy_filling_1 = model.addCopy(Part_1_doc, [Filling_1.result()], 1)
 107 Copy_filling_1.result().setName("Filling_1_copie")
 108
 109 Filling_2_objects = \
 110     [model.selection(Sketch_1.defaultResult(), SketchCircle_1_ext.defaultResult().shape()), \
 111      model.selection(Sketch_2.defaultResult(), SketchCircle_2_ext.defaultResult().shape()), \
 112      model.selection(Sketch_3.defaultResult(), SketchCircle_3_ext.defaultResult().shape())]
 113 Filling_2 = model.addFilling(Part_1_doc, Filling_2_objects, "curve_info", 2, 5, 0, 0.0001, 0.0001, False)
 114
 115 Filling_3 = model.addFilling(Part_1_doc, \
 116         [model.selection(Sketch_1.defaultResult(), SketchCircle_1_int.defaultResult().shape()), \
 117         model.selection(Sketch_1.defaultResult(), SketchCircle_1_ext.defaultResult().shape())], \
 118         "curve_info", 2, 5, 0, 0.0001, 0.0001, False)
 119
 120 Filling_4 = model.addFilling(Part_1_doc, \
 121         [model.selection(Sketch_3.defaultResult(), SketchCircle_3_int.defaultResult().shape()), \
 122         model.selection(Sketch_3.defaultResult(), SketchCircle_3_ext.defaultResult().shape())], \
 123         "curve_info", 2, 5, 0, 0.0001, 0.0001, False)
 124
 125 Face_inf = model.addFace(Part_1_doc, [model.selection(Sketch_1.defaultResult(), SketchCircle_1_int.defaultResult().shape())])
 126
 127 Face_sup = model.addFace(Part_1_doc, [model.selection(Sketch_3.defaultResult(), SketchCircle_3_int.defaultResult().shape())])
 128
 129 liste_tube_ext = [Filling_1.result(), Filling_2.result(), Filling_3.result(), Filling_4.result()]
 130 Solid_tube_ext = model.addSolid(Part_1_doc, liste_tube_ext)
 131
 132 # Recherche de la face du tube sur laquelle faire un congé de raccordement (fillet) :
 133 face = Filling_4.defaultResult().shape().face()
 134 exp = GeomAPI_ShapeExplorer(Solid_tube_ext.defaultResult().shape(), GeomAPI_Shape.FACE)
 135 while exp.more():
 136    cur = exp.current().face()
 137    if face.isEqual(cur) : # and face.isSameGeometry(cur):
 138       res = cur
 139       break
 140    exp.next()
 141 #print(type(res))
 142 Fillet_1 = model.addFillet(Part_1_doc, [model.selection(Solid_tube_ext.defaultResult(), res)], "fillet_radius_top", keepSubResults = False)
 143 Solid_tube_ext = Fillet_1
 144
 145 #=============================================================================
 146 # Création du tube intérieur avec lequel on va couper la sphère extérieure :
 147 #=============================================================================
 148 liste_tube_int = [Copy_filling_1.result(), Face_inf.result(), Face_sup.result()]
 149 Solid_tube_int = model.addSolid(Part_1_doc, liste_tube_int)
 150 model.do()
 151
 152 # Copie des 2 tubes intérieur et extérieur avant qu'ils ne soient avalés par les
 153 # différentes opérations ultérieures (cut, fuse) :
 154 Sphere_ext.result().setName("Sphere_ext")
 155 Solid_tube_ext.result().setName("Solid_tube_ext")
 156 Solid_tube_int.result().setName("Solid_tube_int")
 157 Copie_tube_ext = model.addCopy(Part_1_doc, [Solid_tube_ext.result()], 1)
 158 Copie_tube_ext.result().setName("Tube_ext_1")
 159 Copie_tube_int = model.addCopy(Part_1_doc, [Solid_tube_int.result()], 1)
 160 Copie_tube_int.result().setName("Tube_int_1")
 161
 162 #=============================================================================
 163 # Création d'un congé de raccordement à la base du premier tube :
 164 #=============================================================================
 165 # Calcul de l'intersection entre tube ext. et sphère de manière à identifier l'arête sur # laquelle on va faire un congé (fillet) :
 166 Intersection_1 = model.addIntersection(Part_1_doc, \
 167                 [model.selection("SOLID", "Sphere_ext"), model.selection("SOLID", "Solid_tube_ext")], keepSubResults = True)
 168 Recover_1 = model.addRecover(Part_1_doc, Intersection_1, [Sphere_ext.result(), Solid_tube_ext.result()])
 169 Solid_tube_ext = Recover_1.results()[0] #.setName("Tube_ext_1")
 170 Sphere_ext = Recover_1.results()[1] #.setName("Sphere_ext")
 171
 172 # Coupe de la sphère par tube int. et fusion avec tube ext. :
 173 Cut_1 = model.addCut(Part_1_doc, [Sphere_ext], [Solid_tube_int.result()], keepSubResults = False)
 174 Fuse_1 = model.addFuse(Part_1_doc, [Solid_tube_ext, Cut_1.result()], [], keepSubResults = False)
 175 model.do()
 176 #=======================================
 177 # Fonction de calcul d'un rayon sur un cercle ou arc de cercle à partir de 3 points :
 178 # 2 points extrêmes et point milieu. Elle nous sert à sélectionner le cercle extérieur
 179 # de l'intersection tube-sphère, qui correspond au cercle sur lequel on veut faire
 180 # un congé de raccordement (fillet).
 181 def radius(theEdge):
 182 #=======================================
 183     p1 = theEdge.firstPoint()
 184     p2 = theEdge.middlePoint()
 185     p3 = theEdge.lastPoint()
 186     dir12 = GeomAPI_Dir(p1.xyz().decreased(p2.xyz()))
 187     dir32 = GeomAPI_Dir(p3.xyz().decreased(p2.xyz()))
 188     ang = 0.5 * dir12.angle(dir32)
 189     return 0.5 * p1.distance(p2) / cos(ang)
 190 #=======================================
 191
 192 maxRad=-999
 193 for i in range(Intersection_1.result().numberOfSubs()):
 194         # The shapes are not circles, but rather NURBS :
 195         rad = radius(GeomAPI_Edge(Intersection_1.result().subResult(i).resultSubShapePair()[1]))
 196         if rad > maxRad:
 197                 imax=i
 198                 maxRad = rad
 199 #print(Intersection_1.result().subResult(imax).name())
 200 edge = GeomAPI_Edge(Intersection_1.result().subResult(imax).resultSubShapePair()[1])
 201 exp = GeomAPI_ShapeExplorer(Fuse_1.defaultResult().shape(), GeomAPI_Shape.EDGE)
 202 while exp.more():
 203    cur = exp.current().edge()
 204    if edge.isEqual(cur) : # and edge.isSameGeometry(cur):
 205       resEdge = cur
 206       break
 207    exp.next()
 208 #print(type(res))
 209 Sphere_ext = model.addFillet(Part_1_doc, [model.selection(Fuse_1.defaultResult(), resEdge)], "fillet_radius_bottom", keepSubResults = False)
 210 #Sphere_ext = Fillet_2
 211 Sphere_ext.result().setName("Sphere_ext")
 212 model.do()
 213 #=============================================================================
 214 # Fin création congé de racc. à la base du premier tube.
 215 #=============================================================================
 216 f2=Sphere_ext # mémorisation pour rangement dans folder à la fin du script
 217
 218 RemoveResults_1 = model.addRemoveResults(Part_1_doc, [Fillet_1.result()])
 219 RemoveResults_2 = model.addRemoveResults(Part_1_doc, [Intersection_1.result()])
 220 model.do()
 221 #=============================================================================
 222 # Boucle de répétition des tubes autour de la sphère :
 223 #=============================================================================
 224 i=1
 225 CutTools = [] # liste pour accumulation des tubes int. et cut de la sphère ext. par tous ces tubes int.
 226 IntTools = [Sphere_ext.result()] # liste pour accumulation de la sphère ext. et de tous les tubes ext. et calcul d'intersection entre tous ces objets (pour déterminer les arêtes sur lesquelles on mettra des fillets).
 227 for x,y,z in zip(listeX[1:], listeY[1:], listeZ[1:]):
 228         i += 1
 229 # on parcourt les listes à partir du 2ème élément car le premier point (Point_init) a déjà été créé.
 230         CurPoint = model.addPoint(Part_1_doc, x,y,z)
 231         if i==2:
 232                 f3 = CurPoint # mémorisation de la feature pour insertion ultérieure dans les dossiers (folders)
 233 #       print("Boucle",i,":") #," : (",round(x,1), round(y,1), round(z,1),")")
 234         Copie_tube_ext = model.addCopy(Part_1_doc, [Solid_tube_ext], 1)
 235         Copie_tube_ext.result().setName("Tube_ext_"+str(i))
 236         Copie_tube_int = model.addCopy(Part_1_doc, [Solid_tube_int.result()], 1)
 237         Copie_tube_int.result().setName("Tube_int_"+str(i))
 238         Rotation = model.addRotation(Part_1_doc, \
 239                 [Copie_tube_ext.result(), Copie_tube_int.result()], \
 240                 center = model.selection("VERTEX", "PartSet/Origin"), \
 241                 start = Point_init.result(), \
 242                 end = CurPoint.result(), keepSubResults = True)
 243         Copie_tube_ext = Rotation.results()[0]
 244         Copie_tube_int = Rotation.results()[1]
 245
 246         model.do()
 247
 248         CutTools.append(model.selection("SOLID", "Tube_int_"+str(i)))
 249         IntTools.append(Copie_tube_ext)
 250
 251 #======================================================================
 252 # Création d'un congé de raccordement à la base du tube courant (n°i) :
 253 #======================================================================
 254 # Calcul de l'intersection entre tous les tubes ext. et la sphère extérieure
 255 # de manière à identifier l'arête sur  laquelle on va faire un congé (fillet) :
 256 Intersection_1 = model.addIntersection(Part_1_doc, IntTools, keepSubResults = True)
 257 maxRad=-999
 258 #print("=============================================================")
 259 TOLERANCE = 1.e-5
 260 EdgesForFillet = []
 261 for j in range(Intersection_1.result().numberOfSubs()):
 262         edge = GeomAPI_Edge(Intersection_1.result().subResult(j).resultSubShapePair()[1])
 263         rad = radius(edge)
 264 #       print("Edge ",j,":",Intersection_1.result().subResult(j).name(),"-- rayon =",rad)
 265         if rad > maxRad + TOLERANCE:
 266                 jmax=j
 267                 maxRad = rad
 268                 EdgesForFillet = [edge]
 269         elif rad + TOLERANCE > maxRad:
 270                 EdgesForFillet.append(edge)
 271 # Restauration de la sphère ext. et de tous les tubes ext. :
 272 Recover_1 = model.addRecover(Part_1_doc, Intersection_1, IntTools)
 273 Sphere_ext = Recover_1.results()[0]
 274 FuseTools = Recover_1.results()[1:]
 275
 276 # Coupe de la sphère par tous les tubes int. et fusion du résultat avec tous les tubes ext. :
 277 Cut_2 = model.addCut(Part_1_doc, [Sphere_ext], CutTools, keepSubResults = False)
 278 Fuse_2 = model.addFuse(Part_1_doc, [Cut_2.result()], FuseTools, keepSubResults = False)
 279 model.do()
 280
 281 # Détermination des arêtes du bas des tubes sur lesquelles mettre des "fillets" :
 282 #print("Edge n°",jmax,"avec le plus grand rayon :",Intersection_1.result().subResult(jmax).name(),"-- rayon=",maxRad)
 283 FilletEdges = []
 284 exp = GeomAPI_ShapeExplorer(Fuse_2.defaultResult().shape(), GeomAPI_Shape.EDGE)
 285 while exp.more():
 286         cur = exp.current().edge()
 287         for edge in EdgesForFillet:
 288                 if edge.isEqual(cur) : # and edge.isSameGeometry(cur):
 289                         FilletEdges.append(model.selection(Fuse_2.defaultResult(), cur))
 290                         EdgesForFillet.remove(edge)
 291                         break
 292         exp.next()
 293 #print(type(res))
 294 Fillet_2 = model.addFillet(Part_1_doc, FilletEdges, "fillet_radius_bottom", keepSubResults = False)
 295 Fillet_2.setName("Fillets_bottom_tubes")
 296 Fillet_2.result().setName("Resultat_1_tube")
 297 Sphere_ext = Fillet_2
 298 model.addRemoveResults(Part_1_doc, [Intersection_1.result()])
 299 model.do()
 300 ##      #======================================================================
 301 ##      # Fin de la création du congé de raccordement
 302 ##      #======================================================================
 303 ##      Sphere_ext.result().setName("Resultat_"+str(i))
 304 ##      i+=1
 305 ##      model.do()
 306 ###=====================================================
 307 ### Fin de la boucle de répétition des tubes.
 308 ### A ce stade, on a une sphère fusionnée avec les num_pts tubes.
 309 ###=====================================================
 310
 311 # Pour faire le ménage dans les résultats et ne laisser que l'unique solide résultant,
 312 # on supprime les 2 tubes initiaux qui ont servi à la copie dans la boucle de répétition :
 313 model.addRemoveResults(Part_1_doc, [model.selection("SOLID", "Tube_ext_1")])
 314 model.addRemoveResults(Part_1_doc, [model.selection("SOLID", "Tube_int_1")])
 315
 316 ## Finalement, on crée la sphère intérieure avec laquelle on va creuser (cut) la sphère fusionnée
 317 ## avec les num_pts tubes :
 318 Sphere_int = model.addSphere(Part_1_doc, model.selection("VERTEX", "PartSet/Origin"), "Rint")
 319 Sphere_int.setName("Sphere_int")
 320 Cut_2 = model.addCut(Part_1_doc, [Sphere_ext.result()], [Sphere_int.result()], keepSubResults = False)
 321 Cut_2.result().setName("Resultat_final")
 322 Cut_2.result().setColor(0,0,150)
 323 #Cut_2.result().setTransparency(0.5)
 324 model.do()
 325
 326 ## Création de répertoires pour y ranger les features (et compacter l'OB) :
 327 Folder_1 = model.addFolder(Part_1_doc, Point_init, f2)
 328 Folder_1.setName("Tube_initial")
 329 Folder_2 = model.addFolder(Part_1_doc, f3, Fillet_2)
 330 Folder_2.setName("Boucle_copie_"+str(num_pts)+"_tubes")
 331 model.do()
 332 model.end()
 333
 334 #=====================================================
 335 # Tests divers (paramétrage, vérification nombre de solides, faces) :
 336 #=====================================================
 337 model.begin()
 338 Rext.setValue(11);model.do() # 10
 339 longueur.setValue(6);model.do() # 4
 340 Rtube.setValue(0.4);model.do() # 0.8
 341 coef1.setValue(0.3);model.do() # 0.5
 342 coef2.setValue(2);model.do() # 1.2
 343 ep_tube.setValue(0.2);model.do() # 0.3
 344 fillet_radius_top.setValue(0.06);model.do() # 0.12
 345 fillet_radius_bottom.setValue(1);model.do() # 0.8
 346
 347 # Retour aux valeurs initiales sauf coef2 et fillet_radius_bottom :
 348 Rext.setValue(10);model.do() # 10
 349 longueur.setValue(4);model.do() # 4
 350 Rtube.setValue(0.8);model.do() # 0.8
 351 coef1.setValue(0.5);model.do() # 0.5
 352 coef2.setValue(1.5);model.do() # 1.2
 353 ep_tube.setValue(0.3);model.do() # 0.3
 354 fillet_radius_top.setValue(0.12);model.do() # 0.12
 355 fillet_radius_bottom.setValue(0.7);model.do() # 0.8
 356
 357 model.end()
 358
 359 #model.generateTests(Part_1, "Part_1")
 360 model.testNbResults(Part_1, 1)
 361 model.testNbSubResults(Part_1, [0])
 362 model.testNbSubShapes(Part_1, GeomAPI_Shape.SOLID, [1])
 363 #model.testNbSubShapes(Part_1, GeomAPI_Shape.FACE, [358])
 364 model.testResultsVolumes(Part_1, [1000.4175], 5)