src/Tools/t_shape/t_shape_builder.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 sys
  22 import salome
  23
  24 import GEOM
  25 from salome.geom import geomBuilder
  26 import math
  27 import SALOMEDS
  28 import time
  29
  30 geompy = None
  31
  32 def demidisk(r1, a1, roty=0, solid_thickness=0):
  33   if solid_thickness < 1e-7:
  34     with_solid = False
  35   else:
  36     with_solid = True
  37
  38   O = geompy.MakeVertex(0, 0, 0)
  39   OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
  40   OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
  41   OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
  42
  43   v=list(range(8))
  44   l=list(range(8))
  45   v0 = geompy.MakeVertex(0, 0, 0)
  46   v[0] = geompy.MakeVertex(0, r1/2.0, 0)
  47   v[1] = geompy.MakeVertex(0, r1, 0)
  48   l[1] = geompy.MakeLineTwoPnt(v[0], v[1])
  49   l[2] = geompy.MakeRotation(l[1], OX, a1*math.pi/180.0)
  50   v[4] = geompy.MakeRotation(v[0], OX, a1*math.pi/180.0)
  51   v[6] = geompy.MakeRotation(v[1], OX, a1*math.pi/180.0)
  52
  53   v[2] = geompy.MakeVertex(0, -r1/2.0, 0)
  54   v[3] = geompy.MakeVertex(0, -r1, 0)
  55   l[3] = geompy.MakeLineTwoPnt(v[2], v[3])
  56   l[4] = geompy.MakeRotation(l[3], OX, -a1*math.pi/180.0)
  57   v[5] = geompy.MakeRotation(v[2], OX, -a1*math.pi/180.0)
  58   v[7] = geompy.MakeRotation(v[3], OX, -a1*math.pi/180.0)
  59
  60   l[5] = geompy.MakeLineTwoPnt(v[4], v[5])
  61   l[6] = geompy.MakeLineTwoPnt(v[0], v[4])
  62   l[7] = geompy.MakeLineTwoPnt(v[2], v[5])
  63
  64   v7 = geompy.MakeVertex(0, 0, r1)
  65   arc1 = geompy.MakeArc(v[1], v7, v[3])
  66   l[0] = geompy.MakeLineTwoPnt(v[1], v[3])
  67   face1 = geompy.MakeFaceWires([arc1, l[0]], 1)
  68   part1 = geompy.MakePartition([face1], [l[2], l[4], l[5], l[6], l[7]], [], [], geompy.ShapeType["FACE"], 0, [], 0)
  69
  70   if with_solid:
  71     # Add some faces corresponding to the solid layer outside
  72     # the fluid part
  73
  74     # --- Vertices
  75     v0 = geompy.MakeVertex(0, r1 + solid_thickness, 0)
  76     v1 = geompy.MakeRotation(v0, OX, a1*math.pi/180.0)
  77     v2 = geompy.MakeRotation(v0, OX, math.pi - (a1*math.pi/180.0))
  78     v3 = geompy.MakeRotation(v0, OX, math.pi)
  79     v.extend([v0,v1,v3,v2]) # The order is important for use in pointsProjetes
  80     # --- Lines
  81     l0 = geompy.MakeLineTwoPnt(v[1], v0)
  82     l2 = geompy.MakeRotation(l0, OX, a1*math.pi/180.0)
  83     l3 = geompy.MakeRotation(l0, OX, math.pi - (a1*math.pi/180.0))
  84     # --- Faces
  85     face2 = geompy.MakeRevolution(l0, OX, a1*math.pi/180.0)
  86     face3 = geompy.MakeRevolution(l2, OX, math.pi - 2*a1*math.pi/180.0)
  87     face4 = geompy.MakeRevolution(l3, OX, a1*math.pi/180.0)
  88     # --- Compound of the "fluid part" of the divided disk and the additional faces
  89     compound1 = geompy.MakeCompound([part1, face2, face3, face4])
  90     # --- Glue edges
  91     part1 = geompy.MakeGlueEdges(compound1,1e-7)
  92
  93   if roty != 0:
  94     vrot = [ geompy.MakeRotation(vert, OY, roty*math.pi/180.0) for vert in v ]
  95     lrot = [ geompy.MakeRotation(lin, OY, roty*math.pi/180.0) for lin in l ]
  96     arc = geompy.MakeRotation(arc1, OY, roty*math.pi/180.0)
  97     part = geompy.MakeRotation(part1, OY, roty*math.pi/180.0)
  98     return vrot, lrot, arc, part
  99   else:
 100     return v, l, arc1, part1
 101
 102 def pointsProjetes(vref, face):
 103   vface = geompy.ExtractShapes(face, geompy.ShapeType["VERTEX"], True)
 104   vord = list(range(len(vref)))
 105   plan = geompy.MakePlaneThreePnt(vref[0], vref[1], vref[-1], 10000)
 106   vproj = [ geompy.MakeProjection(vert, plan) for vert in vface ]
 107   for i,v in enumerate(vproj):
 108     dist = [ (geompy.MinDistance(v, vr), j) for j,vr in enumerate(vref) ]
 109     dist.sort()
 110     if dist[0][0] < 1.e-3:
 111       vord[dist[0][1]] = vface[i]
 112   return vord
 113
 114 def arcsProjetes(vf, face):
 115   lface = geompy.ExtractShapes(face, geompy.ShapeType["EDGE"], True)
 116   lord = list(range(3))
 117   ends = [vf[1], vf[6], vf[7], vf[3]]
 118   for i in range(3):
 119     for lf in lface:
 120       pts = geompy.ExtractShapes(lf, geompy.ShapeType["VERTEX"], True)
 121       if (((geompy.MinDistance(pts[0], ends[i]) < 0.001) and (geompy.MinDistance(pts[1], ends[i+1]) < 0.001)) or
 122           ((geompy.MinDistance(pts[1], ends[i]) < 0.001) and (geompy.MinDistance(pts[0], ends[i+1]) < 0.001))):
 123         lord[i] = lf
 124         #print "arc_%d OK"%i
 125         break
 126     pass
 127   return lord
 128
 129 def build_shape(r1, r2, h1, h2, solid_thickness=0, progressBar=None ):
 130   """ Builds the final shape """
 131
 132   if progressBar is not None:
 133     time0 = time.time()
 134     print(time.time() -time0)
 135
 136   if solid_thickness < 1e-7:
 137     with_solid = False
 138   else:
 139     with_solid = True
 140
 141   global geompy
 142   geompy = geomBuilder.New()
 143
 144   O = geompy.MakeVertex(0, 0, 0)
 145   OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
 146   OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
 147   OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
 148
 149   a1 = 45.0
 150   seuilmax = 0.1
 151   ratio = float(r2)/float(r1)
 152   if ratio > (1.0 -seuilmax):
 153     a1 = 45.0*(1.0 -ratio)/seuilmax
 154
 155   # --- Creation of the jonction faces
 156   [faci, sect45, arc1, l1, lord90, lord45, edges, arcextru] = jonction(r1, r2,\
 157                                                                        h1, h2, a1)
 158   if progressBar is not None:
 159     progressBar.addSteps(2)
 160     print(time.time() -time0)
 161
 162   if with_solid:
 163     # The same code is executed again with different external radiuses in order
 164     # to get the needed faces and edges to build the solid layer of the pipe
 165     [faci_ext, sect45_ext, arc1_ext, l1_ext, \
 166      lord90_ext, lord45_ext, edges_ext, arcextru_ext] = jonction(r1 + solid_thickness, r2 + solid_thickness,\
 167                                                                  h1, h2, a1)
 168     faces_jonction_ext = []
 169     for i,l in enumerate(lord90):
 170       faces_jonction_ext.append(geompy.MakeQuad2Edges(lord90[i],lord90_ext[i]))
 171     for i in [1, 3, 6, 7]:
 172       faces_jonction_ext.append(geompy.MakeQuad2Edges(edges[i],edges_ext[i]))
 173     for i,l in enumerate(lord45):
 174       faces_jonction_ext.append(geompy.MakeQuad2Edges(lord45[i],lord45_ext[i]))
 175
 176   if progressBar is not None:
 177     progressBar.addSteps(4)
 178     print(time.time() -time0)
 179
 180   # --- extrusion droite des faces de jonction, pour reconstituer les demi cylindres
 181   if with_solid:
 182     sect45 = geompy.MakeCompound([sect45]+faces_jonction_ext[-3:])
 183     sect45 = geompy.MakeGlueEdges(sect45, 1e-7)
 184
 185   if progressBar is not None:
 186     progressBar.addSteps(1)
 187     print(time.time() -time0)
 188
 189   extru1 = geompy.MakePrismVecH(sect45, OX, h1+10)
 190
 191   faces_coupe = faci[5:]
 192   if with_solid:
 193     faces_coupe = faci[5:]+faces_jonction_ext[:3]
 194   base2 = geompy.MakePartition(faces_coupe, [], [], [], geompy.ShapeType["FACE"], 0, [], 0)
 195   extru2 = geompy.MakePrismVecH(base2, OZ, h2)
 196
 197   if progressBar is not None:
 198     progressBar.addSteps(1)
 199     print(time.time() -time0)
 200
 201   # --- partition et coupe
 202
 203   if with_solid:
 204      demiDisque = geompy.MakeFaceWires([arc1_ext, l1_ext[0]], 1)
 205   else:
 206      demiDisque = geompy.MakeFaceWires([arc1, l1[0]], 1)
 207   demiCylindre = geompy.MakePrismVecH(demiDisque, OX, h1)
 208
 209   if progressBar is not None:
 210     progressBar.addSteps(1)
 211     print(time.time() -time0)
 212
 213   box = geompy.MakeBox(0, -2*(r1+h1), -2*(r1+h1), 2*(r1+h1), 2*(r1+h1), 2*(r1+h1))
 214   rot = geompy.MakeRotation(box, OY, 45*math.pi/180.0)
 215
 216   # NOTE: The following Cut takes almost half of the total execution time
 217   garder = geompy.MakeCutList(demiCylindre, [extru2, rot], True)
 218
 219   if progressBar is not None:
 220     progressBar.addSteps(9)
 221     print(time.time() -time0)
 222
 223   faces_coupe = faci[:5]
 224   if with_solid:
 225     faces_coupe.extend(faces_jonction_ext[-7:])
 226   raccord = geompy.MakePartition([garder], faces_coupe + [arcextru], [], [], geompy.ShapeType["SOLID"], 0, [], 0)
 227   assemblage = geompy.MakeCompound([raccord, extru1, extru2])
 228   assemblage = geompy.MakeGlueFaces(assemblage, 1e-7)
 229
 230   if progressBar is not None:
 231     progressBar.addSteps(3)
 232     print(time.time() -time0)
 233
 234   box = geompy.MakeBox(-1, -(r1+r2+2*solid_thickness), -1, h1, r1+r2+2*solid_thickness, h2)
 235
 236   # NOTE: This operation takes about 1/4 of the total execution time
 237   final = geompy.MakeCommonList([box, assemblage], True)
 238
 239   if progressBar is not None:
 240     progressBar.addSteps(5)
 241     print(time.time() -time0)
 242
 243   # --- Partie inférieure
 244
 245   v3, l3, arc3, part3 = demidisk(r1, a1, 180.0, solid_thickness)
 246   extru3 = geompy.MakePrismVecH(part3, OX, h1)
 247
 248   # --- Symétrie
 249
 250   compound = geompy.MakeCompound([final, extru3])
 251   plane = geompy.MakePlane(O,OX,2000)
 252   compound_mirrored = geompy.MakeMirrorByPlane(compound, plane)
 253   compound_total = geompy.MakeCompound([compound, compound_mirrored])
 254   final = geompy.MakeGlueFaces(compound_total, 1e-07)
 255
 256   if progressBar is not None:
 257     progressBar.addSteps(1)
 258     print(time.time() -time0)
 259
 260   return final
 261
 262
 263 def jonction(r1, r2, h1, h2, a1):
 264   """ Builds the jonction faces and
 265   returns what is needed to build the whole pipe
 266   """
 267
 268   O = geompy.MakeVertex(0, 0, 0)
 269   OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
 270   OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
 271   OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
 272
 273   # --- sections droites des deux demi cylindres avec le partionnement
 274   v1, l1, arc1, part1 = demidisk(r1, a1, 0.)
 275   v2, l2, arc2, part2 = demidisk(r2, a1, 90.0)
 276
 277   # --- extrusion des sections --> demi cylindres de travail, pour en extraire les sections utilisées au niveau du Té
 278   #     et enveloppe cylindrique du cylindre principal
 279
 280   demicyl1 = geompy.MakePrismVecH(part1, OX, h1)
 281   demicyl2 = geompy.MakePrismVecH(part2, OZ, h2)
 282   arcextru = geompy.MakePrismVecH(arc1, OX, h1)
 283
 284   # --- plan de coupe à 45° sur le cylindre principal,
 285   #     section à 45° du cylndre principal,
 286   #     section du cylindre secondaire par l'enveloppe cylindrique du cylindre principal
 287
 288   plan1 = geompy.MakePlane(O, OX, 4*r1)
 289   planr = geompy.MakeRotation(plan1, OY, 45*math.pi/180.0)
 290
 291   sect45 = geompy.MakeCommonList([demicyl1, planr], True)
 292   sect90 = geompy.MakeCommonList([demicyl2, arcextru], True)
 293   #geompy.addToStudy(sect90, "sect90")
 294
 295   # --- liste ordonnée des points projetés sur les deux sections
 296
 297   vord45 = pointsProjetes(v1, sect45)
 298   vord90 = pointsProjetes(v2, sect90)
 299
 300   # --- identification des projections des trois arcs de cercle, sur les deux sections.
 301
 302   lord45 = arcsProjetes(vord45, sect45)
 303   lord90 = arcsProjetes(vord90, sect90)
 304
 305   # --- abaissement des quatre points centraux de la section du cylindre secondaire
 306
 307   dz = -r2/2.0
 308   for i in (0, 2, 4, 5):
 309     vord90[i] = geompy.TranslateDXDYDZ(vord90[i], 0, 0, dz, True)
 310     #geompy.addToStudyInFather(sect90, vord90[i], 'vm%d'%i)
 311
 312   # --- création des deux arêtes curvilignes sur l'enveloppe cylindrique du cylindre principal, à la jonction
 313
 314   curv = [None for i in range(4)] # liaisons entre les points 1, 3, 6 et 7 des 2 sections
 315
 316   curv[0] = geompy.MakeArcCenter(O, vord90[1] , vord45[1], False)
 317   curv[1] = geompy.MakeArcCenter(O, vord90[3] , vord45[3], False)
 318
 319   lipts = ((6, 6, 4), (7, 7, 5))
 320   for i, ipts in enumerate(lipts):
 321     #print i, ipts
 322     p0 = vord90[ipts[0]]
 323     p1 = vord45[ipts[1]]
 324     p2 = vord45[ipts[2]]
 325     plan = geompy.MakePlaneThreePnt(p0, p1, p2, 10000)
 326     #geompy.addToStudy(plan, "plan%d"%i)
 327     section = geompy.MakeSection(plan, arcextru, True)
 328     secpart = geompy.MakePartition([section], [sect45, sect90], [], [], geompy.ShapeType["EDGE"], 0, [], 0)
 329     #geompy.addToStudy(secpart, "secpart%d"%i)
 330     lsec = geompy.ExtractShapes(secpart, geompy.ShapeType["EDGE"], True)
 331
 332     for l in lsec:
 333       pts = geompy.ExtractShapes(l, geompy.ShapeType["VERTEX"], True)
 334       if (((geompy.MinDistance(pts[0], p0) < 0.001) and (geompy.MinDistance(pts[1], p1) < 0.001)) or
 335           ((geompy.MinDistance(pts[1], p0) < 0.001) and (geompy.MinDistance(pts[0], p1) < 0.001))):
 336         curv[i+2] =l
 337         #print "curv_%d OK"%i
 338         break
 339
 340   # --- creation des arêtes droites manquantes, des faces et volumes pour les quatre volumes de la jonction
 341
 342   edges = [None for i in range(8)]
 343   edges[0] = geompy.MakeLineTwoPnt(vord45[0], vord90[0])
 344   edges[1] = curv[0]
 345   edges[2] = geompy.MakeLineTwoPnt(vord45[2], vord90[2])
 346   edges[3] = curv[1]
 347   edges[4] = geompy.MakeLineTwoPnt(vord45[4], vord90[4])
 348   edges[5] = geompy.MakeLineTwoPnt(vord45[5], vord90[5])
 349   edges[6] = curv[2]
 350   edges[7] = curv[3]
 351
 352   ed45 = [None for i in range(8)]
 353   ed45[0] = geompy.MakeLineTwoPnt(vord45[0], vord45[2])
 354   ed45[1] = geompy.MakeLineTwoPnt(vord45[0], vord45[1])
 355   ed45[2] = geompy.MakeLineTwoPnt(vord45[4], vord45[6])
 356   ed45[3] = geompy.MakeLineTwoPnt(vord45[2], vord45[3])
 357   ed45[4] = geompy.MakeLineTwoPnt(vord45[5], vord45[7])
 358   ed45[5] = geompy.MakeLineTwoPnt(vord45[4], vord45[5])
 359   ed45[6] = geompy.MakeLineTwoPnt(vord45[0], vord45[4])
 360   ed45[7] = geompy.MakeLineTwoPnt(vord45[2], vord45[5])
 361
 362   ed90 = [None for i in range(8)]
 363   ed90[0] = geompy.MakeLineTwoPnt(vord90[0], vord90[2])
 364   ed90[1] = geompy.MakeLineTwoPnt(vord90[0], vord90[1])
 365   ed90[2] = geompy.MakeLineTwoPnt(vord90[4], vord90[6])
 366   ed90[3] = geompy.MakeLineTwoPnt(vord90[2], vord90[3])
 367   ed90[4] = geompy.MakeLineTwoPnt(vord90[5], vord90[7])
 368   ed90[5] = geompy.MakeLineTwoPnt(vord90[4], vord90[5])
 369   ed90[6] = geompy.MakeLineTwoPnt(vord90[0], vord90[4])
 370   ed90[7] = geompy.MakeLineTwoPnt(vord90[2], vord90[5])
 371
 372   faci = []
 373   faci.append(geompy.MakeFaceWires([ed45[6], edges[0], ed90[6], edges[4]], 0))
 374   faci.append(geompy.MakeFaceWires([ed45[7], edges[2], ed90[7], edges[5]], 0))
 375   faci.append(geompy.MakeFaceWires([ed45[2], edges[4], ed90[2], edges[6]], 0))
 376   faci.append(geompy.MakeFaceWires([ed45[5], edges[4], ed90[5], edges[5]], 0))
 377   faci.append(geompy.MakeFaceWires([ed45[4], edges[5], ed90[4], edges[7]], 0))
 378   faci.append(geompy.MakeFaceWires([ed90[0], ed90[6],  ed90[5], ed90[7]], 0))
 379   faci.append(geompy.MakeFaceWires([ed90[1], ed90[6],  ed90[2], lord90[0]], 0))
 380   faci.append(geompy.MakeFaceWires([ed90[2], ed90[5],  ed90[4], lord90[1]], 0))
 381   faci.append(geompy.MakeFaceWires([ed90[3], ed90[7],  ed90[4], lord90[2]], 0))
 382
 383   return faci, sect45, arc1, l1, lord90, lord45, edges, arcextru
 384
 385 def test_t_shape_builder():
 386   """For testing purpose"""
 387   salome.salome_init()
 388   geompy = geomBuilder.New()
 389   for r1 in [1., 100.]:
 390     for r2 in [0.9*r1, 0.5*r1, 0.1*r1, 0.05*r1]:
 391       for thickness in [r1/100., r1/10., r1/2.]:
 392         print(r1, r2, thickness)
 393         h1 = r1 * 2.0
 394         h2 = h1
 395         try:
 396           res = build_shape(r1, r2, h1, h2, thickness)
 397           geompy.addToStudy(res, "res_%f_%f_%f"%(r1,r2, thickness))
 398         except:
 399           print("problem with res_%f_%f_%f"%(r1,r2, thickness))
 400
 401 if __name__=="__main__":
 402   """For testing purpose"""
 403   test_t_shape_builder()