1 # Copyright (C) 2007-2014 CEA/DEN, EDF R&D, OPEN CASCADE
3 # Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 # CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 # This library is free software; you can redistribute it and/or
7 # modify it under the terms of the GNU Lesser General Public
8 # License as published by the Free Software Foundation; either
9 # version 2.1 of the License, or (at your option) any later version.
11 # This library is distributed in the hope that it will be useful,
12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 # Lesser General Public License for more details.
16 # You should have received a copy of the GNU Lesser General Public
17 # License along with this library; if not, write to the Free Software
18 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
24 # In this file are all the generation functions for manipulating the different created macro-objects
31 from salome.geom import geomBuilder
32 geompy = geomBuilder.New( Config.theStudy )
34 from salome.smesh import smeshBuilder
35 smesh = smeshBuilder.New( Config.theStudy )
37 ##########################################################################################################
39 def Box11 (MacObject):
40 if Config.debug : print "Generating regular box"
42 dummy1 = geompy.MakeScaleAlongAxes( ElemBox11 (), None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
43 RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
45 MacObject.GeoChildren.append(RectFace)
46 MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
48 if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
51 MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh
52 Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis
54 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
55 Reg1D = MacObject.Mesh[0].Segment()
56 Reg1D.NumberOfSegments(MacObject.MeshPar[0])
58 MacObject.Mesh[0].Compute() # Generates the mesh
60 MacObject.DirectionalMeshParams = [MacObject.MeshPar[0],MacObject.MeshPar[0],MacObject.MeshPar[0],MacObject.MeshPar[0]]
63 Config.ListObj.append(MacObject)
66 ##########################################################################################################
68 def Box42 (MacObject):
69 if Config.debug : print "Generating box 4-2 reducer"
71 Z_Axis = geompy.MakeVectorDXDYDZ(0., 0., 1.)
72 RotAngle = {'SN' : lambda : 0,
73 'NS' : lambda : math.pi,
74 'EW' : lambda : math.pi/2,
75 'WE' : lambda : -math.pi/2, }[MacObject.MeshPar[1]]()
76 dummy0 = geompy.MakeRotation( ElemBox42 () , Z_Axis, RotAngle )
77 dummy1 = geompy.MakeScaleAlongAxes( dummy0, None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
78 RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
80 MacObject.GeoChildren.append(RectFace)
81 MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
83 if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
86 MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh
87 Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis
89 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
90 Reg1D = MacObject.Mesh[0].Segment()
91 Reg1D.NumberOfSegments(MacObject.MeshPar[0])
93 MacObject.Mesh[0].Compute() # Generates the mesh
97 x = MacObject.MeshPar[0]
98 MacObject.DirectionalMeshParams = {'SN' : lambda : [3*x, 3*x, 4*x, 2*x],
99 'NS' : lambda : [3*x, 3*x, 2*x, 4*x],
100 'EW' : lambda : [2*x, 4*x, 3*x, 3*x],
101 'WE' : lambda : [4*x, 2*x, 3*x, 3*x], }[MacObject.MeshPar[1]]()
103 Config.ListObj.append(MacObject)
107 ##########################################################################################################
109 def BoxAng32 (MacObject):
110 if Config.debug : print "Generating sharp angle"
111 Z_Axis = geompy.MakeVectorDXDYDZ(0., 0., 1.)
112 RotAngle = {'NE' : lambda : 0,
113 'NW' : lambda : math.pi/2,
114 'SW' : lambda : math.pi,
115 'SE' : lambda : -math.pi/2, }[MacObject.MeshPar[1]]()
116 dummy0 = geompy.MakeRotation( ElemEdge32 () , Z_Axis, RotAngle )
117 dummy1 = geompy.MakeScaleAlongAxes( dummy0, None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
118 RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
120 MacObject.GeoChildren.append(RectFace)
121 MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
123 if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
126 MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh
127 Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis
129 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
130 Reg1D = MacObject.Mesh[0].Segment()
131 Reg1D.NumberOfSegments(MacObject.MeshPar[0])
133 MacObject.Mesh[0].Compute() # Generates the mesh
137 x = MacObject.MeshPar[0]
138 MacObject.DirectionalMeshParams = {'NE' : lambda : [3*x, 2*x, 3*x, 2*x],
139 'NW' : lambda : [2*x, 3*x, 3*x, 2*x],
140 'SW' : lambda : [2*x, 3*x, 2*x, 3*x],
141 'SE' : lambda : [3*x, 2*x, 2*x, 3*x], }[MacObject.MeshPar[1]]()
143 Config.ListObj.append(MacObject)
145 ##########################################################################################################
146 def CompBox (MacObject):
147 if Config.debug : print "Generating composite box"
149 dummy1 = geompy.MakeScaleAlongAxes( ElemBox11 (), None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
150 RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
152 MacObject.GeoChildren.append(RectFace)
153 MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
155 if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
158 MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh
159 Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis
161 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
163 ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1])
166 Vec = [(1,0,0),(0,1,0)]
169 if IsParallel(Edge,Vec[i]):
170 if not Reference[i]: # If this is the first found edge to be parallel to this direction, apply user preferences for meshing
172 ApplyConstant1DMesh(MacObject.Mesh[0],Edge,int(round(ReducedRatio[i]*MacObject.MeshPar[0])))
174 else: # If there already exists an edge parallel to this direction, then use a 1D projection
175 Apply1DProjMesh(MacObject.Mesh[0],Edge,Reference[i])
178 MacObject.Mesh[0].Compute() # Generates the mesh
180 MacObject.DirectionalMeshParams = [MacObject.MeshPar[0]*ReducedRatio[1],MacObject.MeshPar[0]*ReducedRatio[1],MacObject.MeshPar[0]*ReducedRatio[0],MacObject.MeshPar[0]*ReducedRatio[0]]
183 Config.ListObj.append(MacObject)
186 ##########################################################################################################
188 def CompBoxF (MacObject):
189 if Config.debug : print "Generating composite box"
191 dummy1 = geompy.MakeScaleAlongAxes( ElemBox11 (), None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
192 RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
194 MacObject.GeoChildren.append(RectFace)
195 MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
197 if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
200 MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh
201 Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis
203 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
205 #ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1])
208 Vec = [(1,0,0),(0,1,0)]
211 if IsParallel(Edge,Vec[i]):
212 if not Reference[i]: # If this is the first found edge to be parallel to this direction, apply user preferences for meshing
214 ApplyConstant1DMesh(MacObject.Mesh[0],Edge,int(round(MacObject.MeshPar[0][i])))
216 else: # If there already exists an edge parallel to this direction, then use a 1D projection
217 Apply1DProjMesh(MacObject.Mesh[0],Edge,Reference[i])
220 MacObject.Mesh[0].Compute() # Generates the mesh
222 MacObject.DirectionalMeshParams = [MacObject.MeshPar[0][1],MacObject.MeshPar[0][1],MacObject.MeshPar[0][0],MacObject.MeshPar[0][0]]
225 Config.ListObj.append(MacObject)
227 ##########################################################################################################
229 def NonOrtho (MacObject):
230 if Config.debug : print "Generating Non-orthogonal quadrangle"
232 RectFace = Quadrangler (MacObject.PtCoor)
234 MacObject.GeoChildren.append(RectFace)
235 MacObject.GeoChildrenNames.append("Quad_"+ str(len(Config.ListObj)+1))
238 if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
241 MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh
242 Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis
244 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
246 #ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1])
248 Vec = [MacObject.DirVectors(i) for i in range(4)]
250 Dir = [IsParallel(Edge,Vec[j]) for j in range(4)].index(True)
251 DirConv = [0,0,1,1][Dir]
252 ApplyConstant1DMesh(MacObject.Mesh[0],Edge,int(round(MacObject.MeshPar[0][DirConv])))
254 MacObject.Mesh[0].Compute() # Generates the mesh
256 MacObject.DirectionalMeshParams = [MacObject.MeshPar[0][1],MacObject.MeshPar[0][1],MacObject.MeshPar[0][0],MacObject.MeshPar[0][0]]
259 Config.ListObj.append(MacObject)
262 ##########################################################################################################
264 def QuartCyl (MacObject):
265 if Config.debug : print "Generating quarter cylinder"
266 Z_Axis = geompy.MakeVectorDXDYDZ(0., 0., 1.)
267 RotAngle = {'NE' : lambda : 0,
268 'NW' : lambda : math.pi/2,
269 'SW' : lambda : math.pi,
270 'SE' : lambda : -math.pi/2, }[MacObject.MeshPar[1]]()
271 dummy0 = geompy.MakeRotation( ElemQuartCyl(MacObject.MeshPar[2]) , Z_Axis, RotAngle )
272 dummy1 = geompy.MakeScaleAlongAxes( dummy0, None , MacObject.GeoPar[1][0]/10., MacObject.GeoPar[1][1]/10., 1)
273 RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
275 MacObject.GeoChildren.append(RectFace)
276 MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
278 if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
281 MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh
282 Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis
284 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
285 Reg1D = MacObject.Mesh[0].Segment()
287 #if MacObject.MeshPar[0] == 2 and MacObject.MeshPar[2] <= 2.:
288 # print("Due to a bug in Salome 6.3, we are forced to either increase or decrease the local refinement by 50%, we choose in this case to increase the model's refinement.")
289 # MacObject.MeshPar[0] = 3
291 Reg1D.NumberOfSegments(MacObject.MeshPar[0])
293 MacObject.Mesh[0].Compute() # Generates the mesh
297 x = MacObject.MeshPar[0]
298 N = QuarCylParam(MacObject.MeshPar[2])+1
300 MacObject.DirectionalMeshParams = {'NE' : lambda : [2*x, N*x, 2*x, N*x],
301 'NW' : lambda : [N*x, 2*x, 2*x, N*x],
302 'SW' : lambda : [N*x, 2*x, N*x, 2*x],
303 'SE' : lambda : [2*x, N*x, N*x, 2*x], }[MacObject.MeshPar[1]]()
305 Config.ListObj.append(MacObject)
308 ##########################################################################################################
309 # Below this are the elementary calculation/visualization functions
310 ##########################################################################################################
312 def Publish (ObjToPublish,NamesToPublish):
314 for GeoObj in ObjToPublish :
315 geompy.addToStudy(GeoObj,NamesToPublish[i])
318 def IsParallel (Edge, Vector):
320 Function checks whether a given edge object is parallel to a reference vector.
321 Output can be 0 (not parallel) or 1 (parallel and same sense) or 2 (parallel and opposite sense).
322 If the reference vector is null, the function returns 0
324 if Vector == (0,0,0) : return 0
326 P1 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,0))
327 P2 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,1))
328 V0 = [ P1[0] - P2[0], P1[1] - P2[1], P1[2] - P2[2] ]
329 if Distance2Pt((0,0,0),CrossProd(V0,Vector))<1e-7 and DotProd(V0,Vector) > 0 : return 1
330 elif Distance2Pt((0,0,0),CrossProd(V0,Vector))<1e-7 and DotProd(V0,Vector) < 0 : return 2
333 def IsOnCircle (Edge, Center, Radius):
335 Function checks whether a given edge object belong to the periphery of a circle defined by its
337 Output can be 0 (does not belong) or 1 (belongs).
338 If the reference Radius is null, the function returns 0
339 Note that this function is basic in the sense that it only checks if the two border points of a
340 given edge belong to the arc of reference.
342 if Radius == 0 : return 0
344 P1 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,0))
345 P2 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,1))
346 if abs(Distance2Pt(Center,P1)-Radius) < 1e-6 and abs(Distance2Pt(Center,P2)-Radius) < 1e-6:
351 def CrossProd(V1,V2):
353 Determines the cross product of two 3D vectors
355 return ([V1[1]*V2[2]-V1[2]*V2[1], V1[2]*V2[0]-V1[0]*V2[2], V1[0]*V2[1]-V1[1]*V2[0]])
357 def QuarCylParam(PitchRatio):
358 R = float(PitchRatio)/(PitchRatio+1)
360 X = (R+Eps/2.)*math.sin(math.pi/4)+Eps/2.
361 N = int(math.floor((math.pi*R/4.)/(Eps/2.)))
366 Determines the dot product of two 3D vectors
368 if len(V1)==2 : V1.append(0)
369 if len(V2)==2 : V2.append(0)
371 return (V1[0]*V2[0]+V1[1]*V2[1]+V1[2]*V2[2])
373 def Distance2Pt(P1,P2):
375 Returns the distance between two points
377 return (math.sqrt((P1[0]-P2[0])**2+(P1[1]-P2[1])**2+(P1[2]-P2[2])**2))
379 def ApplyConstant1DMesh (ParentMsh, Edge, Nseg):
380 Reg1D = ParentMsh.Segment(geom=Edge)
381 Len = Reg1D.NumberOfSegments(Nseg)
383 def Apply1DProjMesh (ParentMsh, Edge, Ref):
384 Proj1D = ParentMsh.Projection1D(geom=Edge)
385 SrcEdge = Proj1D.SourceEdge(Ref,None,None,None)
387 def EdgeLength (Edge):
389 This function returns the edge object length.
391 P1 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,0))
392 P2 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,1))
393 return Distance2Pt(P1,P2)
397 return Angle*math.pi/180
400 return Angle*180/math.pi
402 def F2D (FloatNumber):
403 return round(FloatNumber*100.)/100.
405 def BezierGen (PointA, PointB, AngleA, AngleB):
407 if AngleA == 0 and AngleB == 0 : return (geompy.MakeEdge(PointA, PointB))
409 A = geompy.PointCoordinates(PointA)
410 B = geompy.PointCoordinates(PointB)
411 dAB = Distance2Pt(A,B)
412 dAC = dAB * (math.tan(AngleA)*math.tan(AngleB)) / (math.sin(AngleA) * ( math.tan(AngleA)+math.tan(AngleB) ) )
413 AngleOX_AB = math.acos((B[0]-A[0])/dAB)
414 PointC = geompy.MakeVertex(A[0]+math.cos(AngleA+AngleOX_AB)*dAC,A[1]+math.sin(AngleA+AngleOX_AB)*dAC,0)
415 CurveACB = geompy.MakeBezier([PointA,PointC,PointB])
418 def GetSideAngleForBezier (PointA , PointB):
420 This function takes for input two points A and B where the bezier line is needed. It calculates the incident
421 angle needed at point A so that the final curve is either at 0 or 90 degrees from the x'Ox axis
423 A = geompy.PointCoordinates(PointA)
424 B = geompy.PointCoordinates(PointB)
426 dAB = Distance2Pt(A,B)
427 Alpha = math.acos(ABx/dAB)
428 #print "New angle request"
429 #print ABx, dAB, R2D(Alpha)
430 if Alpha < math.pi/4 :
431 #print "returning", R2D(-Alpha)
433 elif Alpha < 3*math.pi/4 :
434 #print "returning", R2D(-(Alpha-math.pi/2))
435 return -(Alpha-math.pi/2)
437 #print "returning", R2D(-(Alpha-math.pi))
438 return -(Alpha-math.pi)
440 def VecDivRatio (Vec1, Vec2):
442 This function tries to find the ratio of Vec1 on Vec2 while neglecting any zero term in Vec1. This is used afterwards
443 for determining the global mesh parameter from automatically detected directional mesh params. If no compatibility is
444 possible, the function returns -1
447 for i in range(len(Vec1)) :
448 Vec3.append(float(Vec1[i])/Vec2[i])
451 if not (abs(i)<1e-7) : Ratio.append(i)
453 if min(Ratio) == max(Ratio) and min(Ratio)==int(min(Ratio)) : return(min(Ratio))
459 def ReduceRatio (dx, dy):
461 This function transforms a decimal ratio into a scale between two integers, for example : [0.2,0.05] --> [4,1] ;
465 if isinteger(ratio) : return [1,ratio]
466 elif dx == 1 : # when this function is called recursively!
467 for i in range(1,20) : # searches over 20 decimals
468 if isinteger(ratio * (10**i) ) :
469 Output = GetScale((10**i),int(round(ratio * (10**i) ) ) )
472 for n in range(0,i) :
473 if isinteger(ratio * ( 10**(i)-10**(n) )) :
474 Output = GetScale( 10**(i)-10**(n) , int(round(ratio * ( 10**(i)-10**(n) ) ) ) )
476 if not (Output==[0,0]) : break
479 for i in range(1,10) : # searches over 10 decimals
480 if isinteger(ratio * (10**i) ) :
481 Output = GetScale((10**i),int(round(ratio * (10**i) ) ) )
484 for n in range(0,i) :
485 if isinteger(ratio * ( 10**(i)-10**(n) )) :
486 Output = GetScale( 10**(i)-10**(n) , int(round(ratio * ( 10**(i)-10**(n) ) ) ) )
488 if not (Output==[0,0]) : break
491 print "We are having some trouble while interpreting the following ratio: ",ratio, "\nWe will try a recursive method which may in some cases take some time..."
493 A = ReduceRatio (dx, dy-dx)
494 return ([A[0],A[1]+A[0]])
496 A = ReduceRatio (dy, dx-dy)
497 return ([A[1]+A[0],A[0]])
503 This function is called within ReduceRatio and aims to reduce down two integers X and Y by dividing them with their common divisors;
504 Example: 25 and 5 ---> 5 and 1 / 63 and 12 ---> 21 and 4
507 Divisor = 2 # Initializing the divisor
508 while MaxDiv >= Divisor :
513 MaxDiv = max(MaxDiv,X0)
516 MaxDiv = max(MaxDiv,Y0)
521 Divisor = Divisor + 1
526 This functions applies a simple check if the entered value is an integer
528 x = float('%.5f' % (x)) #Truncate x to 5 digits after the decimal point
529 if math.ceil(x) == math.floor(x) : return True
531 ##########################################################################################
532 # Below this are the functions that create the elementary forms for the macro objects
533 ##########################################################################################
537 This function returns a simple square face of 1 side length
539 RectFace = geompy.MakeFaceHW(1, 1, 1)
544 This function returns a square face of 1 side length, partitioned
545 according to the elementary 4 to 2 reductor method
547 OrigRectFace = geompy.MakeFaceHW(1, 1, 1)
549 SouthPt1 = geompy.MakeVertex (-.25, -.5, 0)
550 SouthPt2 = geompy.MakeVertex (0, -.5, 0)
551 SouthPt3 = geompy.MakeVertex (.25, -.5, 0)
552 WestPt1 = geompy.MakeVertex (-.5, -.5+1./3, 0)
553 WestPt2 = geompy.MakeVertex (-.5, -.5+2./3, 0)
554 EastPt1 = geompy.MakeVertex (.5, -.5+1./3, 0)
555 EastPt2 = geompy.MakeVertex (.5, -.5+2./3, 0)
556 NorthPt = geompy.MakeVertex (0, .5, 0)
557 MidPt1 = geompy.MakeVertex (0, .05, 0)
558 MidPt2 = geompy.MakeVertex (.2, -.18, 0)
559 MidPt3 = geompy.MakeVertex (0, -.28, 0)
560 MidPt4 = geompy.MakeVertex (-.2, -.18, 0)
563 Cutter.append(geompy.MakeEdge(SouthPt2, MidPt3))
564 Cutter.append(geompy.MakeEdge(MidPt1, NorthPt))
565 Cutter.append(BezierGen(SouthPt1, MidPt4, GetSideAngleForBezier(SouthPt1,MidPt4), D2R(15)))
566 Cutter.append(BezierGen(SouthPt3, MidPt2, GetSideAngleForBezier(SouthPt3,MidPt2), D2R(-15)))
567 Cutter.append(BezierGen(WestPt1, MidPt4, GetSideAngleForBezier(WestPt1,MidPt4), D2R(-10)))
568 Cutter.append(BezierGen(EastPt1, MidPt2, GetSideAngleForBezier(EastPt1,MidPt2), D2R(10)))
569 Cutter.append(BezierGen(WestPt2, MidPt1, GetSideAngleForBezier(WestPt2,MidPt1), D2R(-10)))
570 Cutter.append(BezierGen(EastPt2, MidPt1, GetSideAngleForBezier(EastPt2,MidPt1), D2R(10)))
571 Cutter.append(BezierGen(MidPt2, MidPt1, D2R(-15), D2R(-15)))
572 Cutter.append(BezierGen(MidPt3, MidPt2, D2R(10), D2R(15)))
573 Cutter.append(BezierGen(MidPt3, MidPt4, D2R(-10), D2R(-15)))
574 Cutter.append(BezierGen(MidPt4, MidPt1, D2R(15), D2R(15)))
576 RectFace = geompy.MakePartition([OrigRectFace],Cutter, [], [],4, 0, [], 0) #Creating the partition object
578 #for SingleCut in Cutter :
579 # geompy.addToStudy(SingleCut,'Cutter'+str(i))
581 #geompy.addToStudy(RectFace,'RectFace')
586 This function returns a square face of 1 side length, partitioned
587 according to the elementary edge with 3 to 2 reductor
589 OrigRectFace = geompy.MakeFaceHW(1., 1., 1)
591 SouthPt1 = geompy.MakeVertex (-1./6, -0.5, 0.)
592 SouthPt2 = geompy.MakeVertex ( 1./6, -0.5, 0.)
593 WestPt1 = geompy.MakeVertex (-0.5, -1./6, 0.)
594 WestPt2 = geompy.MakeVertex (-0.5, 1./6, 0.)
595 EastPt = geompy.MakeVertex ( 0.5, 0., 0.)
596 NorthPt = geompy.MakeVertex (0., 0.5, 0.)
598 MidPt1 = geompy.MakeVertex (-0.2, -0.2, 0.)
599 MidPt2 = geompy.MakeVertex ( -0.02, -0.02, 0.)
602 Cutter.append(BezierGen(SouthPt1, MidPt1, GetSideAngleForBezier(SouthPt1,MidPt1) , D2R(-5)))
603 Cutter.append(BezierGen( WestPt1, MidPt1, GetSideAngleForBezier(WestPt1 ,MidPt1) , D2R(-5)))
604 Cutter.append(BezierGen(SouthPt2, MidPt2, GetSideAngleForBezier(SouthPt2,MidPt2) , D2R(-10)))
605 Cutter.append(BezierGen( EastPt, MidPt2, GetSideAngleForBezier(EastPt ,MidPt2) , D2R(5)))
606 Cutter.append(BezierGen( WestPt2, MidPt2, GetSideAngleForBezier(WestPt2 ,MidPt2) , D2R(-10)))
607 Cutter.append(BezierGen( MidPt2, NorthPt, GetSideAngleForBezier(NorthPt ,MidPt2) , D2R(-5)))
609 Cutter.append(geompy.MakeEdge(MidPt1, MidPt2))
611 RectFace = geompy.MakePartition([OrigRectFace],Cutter, [], [],4, 0, [], 0) #Creating the partition object
613 #for SingleCut in Cutter :
614 # geompy.addToStudy(SingleCut,'Cutter'+str(i))
616 #geompy.addToStudy(RectFace,'RectFace')
619 def Quadrangler (Points):
621 This function returns a quadranglar face based on four points, non of which 3 are non-colinear.
622 The points are defined by their 2D [(x1,y1),(x2,y2)..] coordinates.
623 Note that the list of points is already arranged upon the creation in MacObject
626 for Point in Points: Pt.append(geompy.MakeVertex(Point[0], Point[1], 0))
627 # The first point is added at the end of the list in order to facilitate the line creation
629 #Draw the lines in order to form the 4 side polygon
631 for i in range(4) : Ln.append(geompy.MakeLineTwoPnt(Pt[i],Pt[i+1]))
632 RectFace = geompy.MakeQuad (Ln[0],Ln[1],Ln[2],Ln[3])
637 This function returns a quarter cylinder to box relay of 1 side length, partitioned
638 with a pitch ratio of K, In other words the side of the box is R*(1+(1/K))
640 R = 10.*float(K)/(K+1)
643 Config.theStudy.SetReal("R" , R)
644 Config.theStudy.SetReal("minusR" , -R)
645 Config.theStudy.SetReal("Eps", Eps)
647 CylWire = geompy.MakeSketcher("Sketcher:F 'R' 0:R 0:L 'Eps':TT 10. 10.0:R 90:L 10.0:R 90:L 'Eps':R 90:C 'minusR' 90.0:WW", [0, 0, 0, 0, 0, 1, 1, 0, -0])
648 CylFace = geompy.MakeFace(CylWire, 1)
650 SouthPt = geompy.MakeVertex (R+Eps/2., 0., 0)
651 SouthWestPt = geompy.MakeVertex ( 0.,0., 0) #The origin can be used for practical partionning objectifs
652 WestPt = geompy.MakeVertex (0., R+Eps/2., 0)
654 N = int(math.floor((math.pi*R/4.)/(Eps/2.)))
655 X = 10.*(1.-1./(N+1))
658 EastPt = geompy.MakeVertex (10.0, X, 0.)
659 NorthPt = geompy.MakeVertex ( X, 10.0, 0.)
661 DivFactor = 8./(F2D(math.log(K))-0.223)
662 #MidPt = geompy.MakeVertex ((R+Eps)*math.cos(math.pi/4), (R+Eps)*math.sin(math.pi/4), 0.)
663 MidPt = geompy.MakeVertex (X-Eps/DivFactor, X-Eps/DivFactor, 0.)
666 Cutter.append(BezierGen(SouthWestPt, MidPt, GetSideAngleForBezier(SouthWestPt,MidPt) , D2R(-5)))
667 Cutter.append(BezierGen( EastPt, MidPt, GetSideAngleForBezier(EastPt,MidPt) , D2R(5)))
668 Cutter.append(BezierGen( MidPt, NorthPt, (-1)**((K<1.25)*1)*D2R(-5), GetSideAngleForBezier(NorthPt,MidPt)))
669 SMBezier = BezierGen( SouthPt, MidPt, GetSideAngleForBezier(SouthPt ,MidPt) , D2R((K<1.25)*180-5))
670 WMBezier = BezierGen( WestPt, MidPt, GetSideAngleForBezier(WestPt, MidPt) , D2R(-5))
671 Cutter.append(WMBezier)
672 Cutter.append(SMBezier)
674 for i in range(1,N) :
675 # Determining intermediate points on the bezier lines and then performing additional cuts
677 TempAnglePlus = (math.pi/4)*(1+float(i)/N)
678 SectionResult = CutnGroup.Go(WMBezier, [(0,0,0,math.sin(TempAnglePlus),-math.cos(TempAnglePlus),0)], [1], ['Dummy'], 0)
679 TempPt1 = SectionResult[1][0]
680 TempPt11 = geompy.MakeVertex ((N-i)*X/N, 10., 0)
682 TempAngleMinus = (math.pi/4)*(1-float(i)/N)
683 SectionResult = CutnGroup.Go(SMBezier, [(0,0,0,math.sin(TempAngleMinus),-math.cos(TempAngleMinus),0)], [1], ['Dummy'], 0)
684 TempPt2 = SectionResult[1][0]
685 TempPt21 = geompy.MakeVertex (10., (N-i)*X/N, 0)
687 Cutter.append(geompy.MakeEdge(SouthWestPt, TempPt1))
688 Cutter.append(geompy.MakeEdge(SouthWestPt, TempPt2))
689 Cutter.append(geompy.MakeEdge(TempPt1, TempPt11))
690 Cutter.append(geompy.MakeEdge(TempPt2, TempPt21))
692 CylFace = geompy.MakePartition([CylFace],Cutter, [], [],4, 0, [], 0) #Creating the partition object
693 CylFace = geompy.MakeTranslation(CylFace, -5., -5., 0.0)
697 def CompatibilityTest(MacObject):
698 Type = MacObject.Type
700 BaseDirPar = [1,1,1,1]
701 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
702 elif Type == 'Box42' :
703 BaseDirPar = {'SN' : lambda : [3, 3, 4, 2],
704 'NS' : lambda : [3, 3, 2, 4],
705 'EW' : lambda : [2, 4, 3, 3],
706 'WE' : lambda : [4, 2, 3, 3], }[MacObject.MeshPar[1]]()
707 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
708 elif Type == 'BoxAng32' :
709 BaseDirPar = {'NE' : lambda : [3, 2, 3, 2],
710 'NW' : lambda : [2, 3, 3, 2],
711 'SW' : lambda : [2, 3, 2, 3],
712 'SE' : lambda : [3, 2, 2, 3], }[MacObject.MeshPar[1]]()
713 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
714 elif Type == 'CompBox' :
715 #print "dx is: ", MacObject.GeoPar[1][1], ". dy is: ",MacObject.GeoPar[1][0]
716 ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0], MacObject.GeoPar[1][1])
718 BaseDirPar = [ReducedRatio[1], ReducedRatio[1], ReducedRatio[0], ReducedRatio[0]]
719 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
721 elif Type == 'QuartCyl' :
722 N = QuarCylParam(MacObject.MeshPar[2])+1
723 BaseDirPar = {'NE' : lambda : [2, N, 2, N],
724 'NW' : lambda : [N, 2, 2, N],
725 'SW' : lambda : [N, 2, N, 2],
726 'SE' : lambda : [2, N, N, 2], }[MacObject.MeshPar[1]]()
727 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
728 elif Type == 'CompBoxF' :
729 RealRatio = MacObject.GeoPar[1][1]/MacObject.GeoPar[1][0]
732 if MacObject.DirectionalMeshParams[2]+MacObject.DirectionalMeshParams[3] :
733 A = int(max(MacObject.DirectionalMeshParams[2:4]))
734 Xd = int(VecDivRatio([A,0,0,0], [1,1,1,1]))
735 if MacObject.DirectionalMeshParams[0]+MacObject.DirectionalMeshParams[1] :
736 A = int(max(MacObject.DirectionalMeshParams[0:2]))
737 Yd = int(VecDivRatio([0,0,A,0], [1,1,1,1]))
739 if Xd == 0 and Yd : Xd = int(round(Yd/RealRatio))
740 elif Yd == 0 : Yd = int(round(RealRatio*Xd))
743 elif Type == 'NonOrtho' :
744 MeanDX = 0.5*(IntLen(MacObject.DirBoundaries(0))+IntLen(MacObject.DirBoundaries(1)))
745 MeanDY = 0.5*(IntLen(MacObject.DirBoundaries(2))+IntLen(MacObject.DirBoundaries(3)))
746 RealRatio = MeanDY/MeanDX
749 if MacObject.DirectionalMeshParams[2]+MacObject.DirectionalMeshParams[3] :
750 A = int(max(MacObject.DirectionalMeshParams[2:4]))
751 Xd = int(VecDivRatio([A,0,0,0], [1,1,1,1]))
752 if MacObject.DirectionalMeshParams[0]+MacObject.DirectionalMeshParams[1] :
753 A = int(max(MacObject.DirectionalMeshParams[0:2]))
754 Yd = int(VecDivRatio([0,0,A,0], [1,1,1,1]))
756 if Xd == 0 and Yd : Xd = int(round(Yd/RealRatio))
757 elif Yd == 0 : Yd = int(round(RealRatio*Xd))
761 def IntLen (Interval) :
763 This function returns the length of a given interval even if the latter is not sorted correctly.
765 return abs(Interval[1]-Interval[0])
767 def NextTo (RefBox, Direction, Extension):
769 This functions returns geometrical parameters for easy positioning of neighbouring objects.
770 The input (RefBox) and output are in the form : [(X0,Y0),(DX,DY)]
777 DirectionalCoef = {'Above' : lambda : [ 0, 1],
778 'Below' : lambda : [ 0,-1],
779 'Right' : lambda : [ 1, 0],
780 'Left ' : lambda : [-1, 0], }[Direction]()
782 X0_1 = X0_0+ DirectionalCoef[0] * (DX_0/2.+Extension/2.)
783 DX_1 = abs(DirectionalCoef[0]) * (Extension) + abs(DirectionalCoef[1])*DX_0
784 Y0_1 = Y0_0+ DirectionalCoef[1] * (DY_0/2.+Extension/2.)
785 DY_1 = abs(DirectionalCoef[1]) * (Extension) + abs(DirectionalCoef[0])*DY_0
787 return [(X0_1,Y0_1),(DX_1,DY_1)]
789 def GeomMinMax (PtA, PtB):
791 This function returns geometrical parameters in the format [(X0,Y0),(DX,DY)]. The input being
792 the coordinates of two points (Xa,Ya), (Xb,Yb).
794 # First test that the vector relying the two points is oblique
795 AB = [PtB[0]- PtA[0],PtB[1]- PtA[1]]
797 print ("Error: the two points are not correctly defined. In the orthonormal system XOY, it is impossible to define a rectangle with these two points")
800 X0 = 0.5*(PtA[0]+PtB[0])
801 Y0 = 0.5*(PtA[1]+PtB[1])
804 return [(X0,Y0),(DX,DY)]
806 def AddIfDifferent (List, Element):
807 if not(Element in List):
808 List = List+(Element,)
811 def IndexMultiOcc (Array,Element) :
813 This functions returns the occurrences indices of Element in Array.
814 As opposed to Array.index(Element) method, this allows determining
815 multiple entries rather than just the first one!
818 try : Array.index(Element)
819 except ValueError : print "No more occurrences"
820 else : Output.append(Array.index(Element))
822 if not(Output == []) and len(Array) > 1 :
823 for index, ArrElem in enumerate(Array[Output[0]+1:]) :
824 if ArrElem == Element : Output.append(index+Output[0]+1)
828 def SortList (ValList, CritList):
830 SortedCritList = copy.copy(CritList)
831 SortedCritList.sort()
832 for i in range(0,len(ValList)):
834 if not(SortedCritList[i]==SortedCritList[i-1]):
835 index = IndexMultiOcc(CritList,SortedCritList[i])
836 Output= Output + [ValList[j] for j in index]
838 index = IndexMultiOcc(CritList,SortedCritList[i])
839 Output= Output + [ValList[j] for j in index]
843 def SortPoints(Points):
845 This function sorts a list of the coordinates of N points as to start at
846 an origin that represents Xmin and Xmax and then proceed in a counter
850 Xmin = min([Points[i][0] for i in range(NbPts)])
851 Ymin = min([Points[i][1] for i in range(NbPts)])
852 Xmax = max([Points[i][0] for i in range(NbPts)])
853 Ymax = max([Points[i][1] for i in range(NbPts)])
854 Crit = [(abs(Point[0]-Xmin)+0.1*(Xmax-Xmin))*(abs(Point[1]-Ymin)+0.1*(Ymax-Ymin)) for Point in Points]
855 #print "Input Points : ", Points
856 #print "Sorting Criterion : ", Crit
857 Order = SortList (range(NbPts), Crit)
858 #print "Sorted Results : ", Order
860 Output.append(Points[Order[0]])
862 Point0 = Points[Order[0]]
863 #print "Reference point :", Point0
865 V = [[Point1[0]-Point0[0],Point1[1]-Point0[1]] for Point1 in Points]
866 Cosines = [-(vec[0]-1E-10)/(math.sqrt(DotProd(vec,vec)+1e-25)) for vec in V]
867 #print "Cosines criterion :", Cosines
868 Order = SortList(range(NbPts),Cosines)
869 #print "Ordered points:", Order
870 for PtIndex in Order[:-1]: Output.append(Points[PtIndex])