src/Tools/t_shape/t_shape_builder.py

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