1 # -*- coding: utf-8 -*-
2 # Copyright (C) 2014-2016 EDF R&D
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.
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.
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
18 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
25 from salome.geom import geomBuilder
32 def demidisk(study, r1, a1, roty=0, solid_thickness=0):
33 if solid_thickness < 1e-7:
38 #geompy = geomBuilder.New(study)
40 O = geompy.MakeVertex(0, 0, 0)
41 OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
42 OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
43 OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
47 v0 = geompy.MakeVertex(0, 0, 0)
48 v[0] = geompy.MakeVertex(0, r1/2.0, 0)
49 v[1] = geompy.MakeVertex(0, r1, 0)
50 l[1] = geompy.MakeLineTwoPnt(v[0], v[1])
51 l[2] = geompy.MakeRotation(l[1], OX, a1*math.pi/180.0)
52 v[4] = geompy.MakeRotation(v[0], OX, a1*math.pi/180.0)
53 v[6] = geompy.MakeRotation(v[1], OX, a1*math.pi/180.0)
55 v[2] = geompy.MakeVertex(0, -r1/2.0, 0)
56 v[3] = geompy.MakeVertex(0, -r1, 0)
57 l[3] = geompy.MakeLineTwoPnt(v[2], v[3])
58 l[4] = geompy.MakeRotation(l[3], OX, -a1*math.pi/180.0)
59 v[5] = geompy.MakeRotation(v[2], OX, -a1*math.pi/180.0)
60 v[7] = geompy.MakeRotation(v[3], OX, -a1*math.pi/180.0)
62 l[5] = geompy.MakeLineTwoPnt(v[4], v[5])
63 l[6] = geompy.MakeLineTwoPnt(v[0], v[4])
64 l[7] = geompy.MakeLineTwoPnt(v[2], v[5])
66 v7 = geompy.MakeVertex(0, 0, r1)
67 arc1 = geompy.MakeArc(v[1], v7, v[3])
68 l[0] = geompy.MakeLineTwoPnt(v[1], v[3])
69 face1 = geompy.MakeFaceWires([arc1, l[0]], 1)
70 part1 = geompy.MakePartition([face1], [l[2], l[4], l[5], l[6], l[7]], [], [], geompy.ShapeType["FACE"], 0, [], 0, True)
73 # Add some faces corresponding to the solid layer outside
77 v0 = geompy.MakeVertex(0, r1 + solid_thickness, 0)
78 v1 = geompy.MakeRotation(v0, OX, a1*math.pi/180.0)
79 v2 = geompy.MakeRotation(v0, OX, math.pi - (a1*math.pi/180.0))
80 v3 = geompy.MakeRotation(v0, OX, math.pi)
81 v.extend([v0,v1,v3,v2]) # The order is important for use in pointsProjetes
83 l0 = geompy.MakeLineTwoPnt(v[1], v0)
84 l2 = geompy.MakeRotation(l0, OX, a1*math.pi/180.0)
85 l3 = geompy.MakeRotation(l0, OX, math.pi - (a1*math.pi/180.0))
87 face2 = geompy.MakeRevolution(l0, OX, a1*math.pi/180.0)
88 face3 = geompy.MakeRevolution(l2, OX, math.pi - 2*a1*math.pi/180.0)
89 face4 = geompy.MakeRevolution(l3, OX, a1*math.pi/180.0)
90 # --- Compound of the "fluid part" of the divided disk and the additional faces
91 compound1 = geompy.MakeCompound([part1, face2, face3, face4])
93 part1 = geompy.MakeGlueEdges(compound1,1e-7)
96 vrot = [ geompy.MakeRotation(vert, OY, roty*math.pi/180.0) for vert in v ]
97 lrot = [ geompy.MakeRotation(lin, OY, roty*math.pi/180.0) for lin in l ]
98 arc = geompy.MakeRotation(arc1, OY, roty*math.pi/180.0)
99 part = geompy.MakeRotation(part1, OY, roty*math.pi/180.0)
100 return vrot, lrot, arc, part
102 return v, l, arc1, part1
104 def pointsProjetes(study, vref, face):
105 #geompy = geomBuilder.New(study)
106 vface = geompy.ExtractShapes(face, geompy.ShapeType["VERTEX"], True)
107 vord = list(range(len(vref)))
108 plan = geompy.MakePlaneThreePnt(vref[0], vref[1], vref[-1], 10000)
109 vproj = [ geompy.MakeProjection(vert, plan) for vert in vface ]
110 for i,v in enumerate(vproj):
111 dist = [ (geompy.MinDistance(v, vr), j) for j,vr in enumerate(vref) ]
113 if dist[0][0] < 1.e-3:
114 vord[dist[0][1]] = vface[i]
117 def arcsProjetes(study, vf, face):
118 #geompy = geomBuilder.New(study)
119 lface = geompy.ExtractShapes(face, geompy.ShapeType["EDGE"], True)
120 lord = list(range(3))
121 ends = [vf[1], vf[6], vf[7], vf[3]]
124 pts = geompy.ExtractShapes(lf, geompy.ShapeType["VERTEX"], True)
125 if (((geompy.MinDistance(pts[0], ends[i]) < 0.001) and (geompy.MinDistance(pts[1], ends[i+1]) < 0.001)) or
126 ((geompy.MinDistance(pts[1], ends[i]) < 0.001) and (geompy.MinDistance(pts[0], ends[i+1]) < 0.001))):
133 def build_shape(study, r1, r2, h1, h2, solid_thickness=0, progressBar=None ):
134 """ Builds the final shape """
136 if progressBar is not None:
138 print(time.time() -time0)
140 if solid_thickness < 1e-7:
146 geompy = geomBuilder.New(study)
148 O = geompy.MakeVertex(0, 0, 0)
149 OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
150 OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
151 OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
155 ratio = float(r2)/float(r1)
156 if ratio > (1.0 -seuilmax):
157 a1 = 45.0*(1.0 -ratio)/seuilmax
159 # --- Creation of the jonction faces
160 [faci, sect45, arc1, l1, lord90, lord45, edges, arcextru] = jonction(study, r1, r2,\
162 if progressBar is not None:
163 progressBar.addSteps(2)
164 print(time.time() -time0)
167 # The same code is executed again with different external radiuses in order
168 # to get the needed faces and edges to build the solid layer of the pipe
169 [faci_ext, sect45_ext, arc1_ext, l1_ext, \
170 lord90_ext, lord45_ext, edges_ext, arcextru_ext] = jonction(study, r1 + solid_thickness, r2 + solid_thickness,\
172 faces_jonction_ext = []
173 for i,l in enumerate(lord90):
174 faces_jonction_ext.append(geompy.MakeQuad2Edges(lord90[i],lord90_ext[i]))
175 for i in [1, 3, 6, 7]:
176 faces_jonction_ext.append(geompy.MakeQuad2Edges(edges[i],edges_ext[i]))
177 for i,l in enumerate(lord45):
178 faces_jonction_ext.append(geompy.MakeQuad2Edges(lord45[i],lord45_ext[i]))
180 if progressBar is not None:
181 progressBar.addSteps(4)
182 print(time.time() -time0)
184 # --- extrusion droite des faces de jonction, pour reconstituer les demi cylindres
186 sect45 = geompy.MakeCompound([sect45]+faces_jonction_ext[-3:])
187 sect45 = geompy.MakeGlueEdges(sect45, 1e-7)
189 if progressBar is not None:
190 progressBar.addSteps(1)
191 print(time.time() -time0)
193 extru1 = geompy.MakePrismVecH(sect45, OX, h1+10)
195 faces_coupe = faci[5:]
197 faces_coupe = faci[5:]+faces_jonction_ext[:3]
198 base2 = geompy.MakePartition(faces_coupe, [], [], [], geompy.ShapeType["FACE"], 0, [], 0, True)
199 extru2 = geompy.MakePrismVecH(base2, OZ, h2)
201 if progressBar is not None:
202 progressBar.addSteps(1)
203 print(time.time() -time0)
205 # --- partition et coupe
208 demiDisque = geompy.MakeFaceWires([arc1_ext, l1_ext[0]], 1)
210 demiDisque = geompy.MakeFaceWires([arc1, l1[0]], 1)
211 demiCylindre = geompy.MakePrismVecH(demiDisque, OX, h1)
213 if progressBar is not None:
214 progressBar.addSteps(1)
215 print(time.time() -time0)
217 box = geompy.MakeBox(0, -2*(r1+h1), -2*(r1+h1), 2*(r1+h1), 2*(r1+h1), 2*(r1+h1))
218 rot = geompy.MakeRotation(box, OY, 45*math.pi/180.0)
220 # NOTE: The following Cut takes almost half of the total execution time
221 garder = geompy.MakeCutList(demiCylindre, [extru2, rot], True)
223 if progressBar is not None:
224 progressBar.addSteps(9)
225 print(time.time() -time0)
227 faces_coupe = faci[:5]
229 faces_coupe.extend(faces_jonction_ext[-7:])
230 raccord = geompy.MakePartition([garder], faces_coupe + [arcextru], [], [], geompy.ShapeType["SOLID"], 0, [], 0, True)
231 assemblage = geompy.MakeCompound([raccord, extru1, extru2])
232 assemblage = geompy.MakeGlueFaces(assemblage, 1e-7)
234 if progressBar is not None:
235 progressBar.addSteps(3)
236 print(time.time() -time0)
238 box = geompy.MakeBox(-1, -(r1+r2+2*solid_thickness), -1, h1, r1+r2+2*solid_thickness, h2)
240 # NOTE: This operation takes about 1/4 of the total execution time
241 final = geompy.MakeCommonList([box, assemblage], True)
243 if progressBar is not None:
244 progressBar.addSteps(5)
245 print(time.time() -time0)
247 # --- Partie inférieure
249 v3, l3, arc3, part3 = demidisk(study, r1, a1, 180.0, solid_thickness)
250 extru3 = geompy.MakePrismVecH(part3, OX, h1)
254 compound = geompy.MakeCompound([final, extru3])
255 plane = geompy.MakePlane(O,OX,2000)
256 compound_mirrored = geompy.MakeMirrorByPlane(compound, plane)
257 compound_total = geompy.MakeCompound([compound, compound_mirrored])
258 final = geompy.MakeGlueFaces(compound_total, 1e-07)
260 if progressBar is not None:
261 progressBar.addSteps(1)
262 print(time.time() -time0)
267 def jonction(study, r1, r2, h1, h2, a1):
268 """ Builds the jonction faces and
269 returns what is needed to build the whole pipe
271 #geompy = geomBuilder.New(study)
273 O = geompy.MakeVertex(0, 0, 0)
274 OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
275 OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
276 OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
278 # --- sections droites des deux demi cylindres avec le partionnement
279 v1, l1, arc1, part1 = demidisk(study, r1, a1, 0.)
280 v2, l2, arc2, part2 = demidisk(study, r2, a1, 90.0)
282 # --- extrusion des sections --> demi cylindres de travail, pour en extraire les sections utilisées au niveau du Té
283 # et enveloppe cylindrique du cylindre principal
285 demicyl1 = geompy.MakePrismVecH(part1, OX, h1)
286 demicyl2 = geompy.MakePrismVecH(part2, OZ, h2)
287 arcextru = geompy.MakePrismVecH(arc1, OX, h1)
289 # --- plan de coupe à 45° sur le cylindre principal,
290 # section à 45° du cylndre principal,
291 # section du cylindre secondaire par l'enveloppe cylindrique du cylindre principal
293 plan1 = geompy.MakePlane(O, OX, 4*r1)
294 planr = geompy.MakeRotation(plan1, OY, 45*math.pi/180.0)
296 sect45 = geompy.MakeCommonList([demicyl1, planr], True)
297 sect90 = geompy.MakeCommonList([demicyl2, arcextru], True)
298 #geompy.addToStudy(sect90, "sect90")
300 # --- liste ordonnée des points projetés sur les deux sections
302 vord45 = pointsProjetes(study, v1, sect45)
303 vord90 = pointsProjetes(study, v2, sect90)
305 # --- identification des projections des trois arcs de cercle, sur les deux sections.
307 lord45 = arcsProjetes(study, vord45, sect45)
308 lord90 = arcsProjetes(study, vord90, sect90)
310 # --- abaissement des quatre points centraux de la section du cylindre secondaire
313 for i in (0, 2, 4, 5):
314 vord90[i] = geompy.TranslateDXDYDZ(vord90[i], 0, 0, dz, True)
315 #geompy.addToStudyInFather(sect90, vord90[i], 'vm%d'%i)
317 # --- création des deux arêtes curvilignes sur l'enveloppe cylindrique du cylindre principal, à la jonction
319 curv = [None for i in range(4)] # liaisons entre les points 1, 3, 6 et 7 des 2 sections
321 curv[0] = geompy.MakeArcCenter(O, vord90[1] , vord45[1], False)
322 curv[1] = geompy.MakeArcCenter(O, vord90[3] , vord45[3], False)
324 lipts = ((6, 6, 4), (7, 7, 5))
325 for i, ipts in enumerate(lipts):
330 plan = geompy.MakePlaneThreePnt(p0, p1, p2, 10000)
331 #geompy.addToStudy(plan, "plan%d"%i)
332 section = geompy.MakeSection(plan, arcextru, True)
333 secpart = geompy.MakePartition([section], [sect45, sect90], [], [], geompy.ShapeType["EDGE"], 0, [], 0, True)
334 #geompy.addToStudy(secpart, "secpart%d"%i)
335 lsec = geompy.ExtractShapes(secpart, geompy.ShapeType["EDGE"], True)
338 pts = geompy.ExtractShapes(l, geompy.ShapeType["VERTEX"], True)
339 if (((geompy.MinDistance(pts[0], p0) < 0.001) and (geompy.MinDistance(pts[1], p1) < 0.001)) or
340 ((geompy.MinDistance(pts[1], p0) < 0.001) and (geompy.MinDistance(pts[0], p1) < 0.001))):
342 #print "curv_%d OK"%i
345 # --- creation des arêtes droites manquantes, des faces et volumes pour les quatre volumes de la jonction
347 edges = [None for i in range(8)]
348 edges[0] = geompy.MakeLineTwoPnt(vord45[0], vord90[0])
350 edges[2] = geompy.MakeLineTwoPnt(vord45[2], vord90[2])
352 edges[4] = geompy.MakeLineTwoPnt(vord45[4], vord90[4])
353 edges[5] = geompy.MakeLineTwoPnt(vord45[5], vord90[5])
357 ed45 = [None for i in range(8)]
358 ed45[0] = geompy.MakeLineTwoPnt(vord45[0], vord45[2])
359 ed45[1] = geompy.MakeLineTwoPnt(vord45[0], vord45[1])
360 ed45[2] = geompy.MakeLineTwoPnt(vord45[4], vord45[6])
361 ed45[3] = geompy.MakeLineTwoPnt(vord45[2], vord45[3])
362 ed45[4] = geompy.MakeLineTwoPnt(vord45[5], vord45[7])
363 ed45[5] = geompy.MakeLineTwoPnt(vord45[4], vord45[5])
364 ed45[6] = geompy.MakeLineTwoPnt(vord45[0], vord45[4])
365 ed45[7] = geompy.MakeLineTwoPnt(vord45[2], vord45[5])
367 ed90 = [None for i in range(8)]
368 ed90[0] = geompy.MakeLineTwoPnt(vord90[0], vord90[2])
369 ed90[1] = geompy.MakeLineTwoPnt(vord90[0], vord90[1])
370 ed90[2] = geompy.MakeLineTwoPnt(vord90[4], vord90[6])
371 ed90[3] = geompy.MakeLineTwoPnt(vord90[2], vord90[3])
372 ed90[4] = geompy.MakeLineTwoPnt(vord90[5], vord90[7])
373 ed90[5] = geompy.MakeLineTwoPnt(vord90[4], vord90[5])
374 ed90[6] = geompy.MakeLineTwoPnt(vord90[0], vord90[4])
375 ed90[7] = geompy.MakeLineTwoPnt(vord90[2], vord90[5])
378 faci.append(geompy.MakeFaceWires([ed45[6], edges[0], ed90[6], edges[4]], 0))
379 faci.append(geompy.MakeFaceWires([ed45[7], edges[2], ed90[7], edges[5]], 0))
380 faci.append(geompy.MakeFaceWires([ed45[2], edges[4], ed90[2], edges[6]], 0))
381 faci.append(geompy.MakeFaceWires([ed45[5], edges[4], ed90[5], edges[5]], 0))
382 faci.append(geompy.MakeFaceWires([ed45[4], edges[5], ed90[4], edges[7]], 0))
383 faci.append(geompy.MakeFaceWires([ed90[0], ed90[6], ed90[5], ed90[7]], 0))
384 faci.append(geompy.MakeFaceWires([ed90[1], ed90[6], ed90[2], lord90[0]], 0))
385 faci.append(geompy.MakeFaceWires([ed90[2], ed90[5], ed90[4], lord90[1]], 0))
386 faci.append(geompy.MakeFaceWires([ed90[3], ed90[7], ed90[4], lord90[2]], 0))
388 return faci, sect45, arc1, l1, lord90, lord45, edges, arcextru
390 def test_t_shape_builder():
391 """For testing purpose"""
393 theStudy = salome.myStudy
394 geompy = geomBuilder.New(theStudy)
395 for r1 in [1., 100.]:
396 for r2 in [0.9*r1, 0.5*r1, 0.1*r1, 0.05*r1]:
397 for thickness in [r1/100., r1/10., r1/2.]:
398 print(r1, r2, thickness)
402 res = build_shape(theStudy, r1, r2, h1, h2, thickness)
403 geompy.addToStudy(res, "res_%f_%f_%f"%(r1,r2, thickness))
405 print("problem with res_%f_%f_%f"%(r1,r2, thickness))
407 if __name__=="__main__":
408 """For testing purpose"""
409 test_t_shape_builder()