1 # Copyright (C) 2014-2022 EDF R&D
3 # This library is free software; you can redistribute it and/or
4 # modify it under the terms of the GNU Lesser General Public
5 # License as published by the Free Software Foundation; either
6 # version 2.1 of the License, or (at your option) any later version.
8 # This library is distributed in the hope that it will be useful,
9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 # Lesser General Public License for more details.
13 # You should have received a copy of the GNU Lesser General Public
14 # License along with this library; if not, write to the Free Software
15 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 # In this file are all the generation functions for manipulating the different created macro-objects
29 from salome.geom import geomBuilder
30 geompy = geomBuilder.New()
32 from salome.smesh import smeshBuilder
33 smesh = smeshBuilder.New()
35 ##########################################################################################################
37 def Box11 (MacObject):
38 if Config.debug : print("Generating regular box")
40 dummy1 = geompy.MakeScaleAlongAxes( ElemBox11 (), None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
41 RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
43 MacObject.GeoChildren.append(RectFace)
44 MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
46 if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
49 MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh
50 Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis
52 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
53 Reg1D = MacObject.Mesh[0].Segment()
54 Reg1D.NumberOfSegments(MacObject.MeshPar[0])
56 MacObject.Mesh[0].Compute() # Generates the mesh
58 MacObject.DirectionalMeshParams = [MacObject.MeshPar[0],MacObject.MeshPar[0],MacObject.MeshPar[0],MacObject.MeshPar[0]]
61 Config.ListObj.append(MacObject)
64 ##########################################################################################################
66 def Box42 (MacObject):
67 if Config.debug : print("Generating box 4-2 reducer")
69 Z_Axis = geompy.MakeVectorDXDYDZ(0., 0., 1.)
70 RotAngle = {'SN' : lambda : 0,
71 'NS' : lambda : math.pi,
72 'EW' : lambda : math.pi/2,
73 'WE' : lambda : -math.pi/2, }[MacObject.MeshPar[1]]()
74 dummy0 = geompy.MakeRotation( ElemBox42 () , Z_Axis, RotAngle )
75 dummy1 = geompy.MakeScaleAlongAxes( dummy0, None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
76 RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
78 MacObject.GeoChildren.append(RectFace)
79 MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
81 if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
84 MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh
85 Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis
87 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
88 Reg1D = MacObject.Mesh[0].Segment()
89 Reg1D.NumberOfSegments(MacObject.MeshPar[0])
91 MacObject.Mesh[0].Compute() # Generates the mesh
95 x = MacObject.MeshPar[0]
96 MacObject.DirectionalMeshParams = {'SN' : lambda : [3*x, 3*x, 4*x, 2*x],
97 'NS' : lambda : [3*x, 3*x, 2*x, 4*x],
98 'EW' : lambda : [2*x, 4*x, 3*x, 3*x],
99 'WE' : lambda : [4*x, 2*x, 3*x, 3*x], }[MacObject.MeshPar[1]]()
101 Config.ListObj.append(MacObject)
105 ##########################################################################################################
107 def BoxAng32 (MacObject):
108 if Config.debug : print("Generating sharp angle")
109 Z_Axis = geompy.MakeVectorDXDYDZ(0., 0., 1.)
110 RotAngle = {'NE' : lambda : 0,
111 'NW' : lambda : math.pi/2,
112 'SW' : lambda : math.pi,
113 'SE' : lambda : -math.pi/2, }[MacObject.MeshPar[1]]()
114 dummy0 = geompy.MakeRotation( ElemEdge32 () , Z_Axis, RotAngle )
115 dummy1 = geompy.MakeScaleAlongAxes( dummy0, None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
116 RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
118 MacObject.GeoChildren.append(RectFace)
119 MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
121 if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
124 MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh
125 Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis
127 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
128 Reg1D = MacObject.Mesh[0].Segment()
129 Reg1D.NumberOfSegments(MacObject.MeshPar[0])
131 MacObject.Mesh[0].Compute() # Generates the mesh
135 x = MacObject.MeshPar[0]
136 MacObject.DirectionalMeshParams = {'NE' : lambda : [3*x, 2*x, 3*x, 2*x],
137 'NW' : lambda : [2*x, 3*x, 3*x, 2*x],
138 'SW' : lambda : [2*x, 3*x, 2*x, 3*x],
139 'SE' : lambda : [3*x, 2*x, 2*x, 3*x], }[MacObject.MeshPar[1]]()
141 Config.ListObj.append(MacObject)
143 ##########################################################################################################
144 def CompBox (MacObject):
145 if Config.debug : print("Generating composite box")
147 dummy1 = geompy.MakeScaleAlongAxes( ElemBox11 (), None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
148 RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
150 MacObject.GeoChildren.append(RectFace)
151 MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
153 if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
156 MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh
157 Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis
159 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
161 ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1])
164 Vec = [(1,0,0),(0,1,0)]
167 if IsParallel(Edge,Vec[i]):
168 if not Reference[i]: # If this is the first found edge to be parallel to this direction, apply user preferences for meshing
170 ApplyConstant1DMesh(MacObject.Mesh[0],Edge,int(round(ReducedRatio[i]*MacObject.MeshPar[0])))
172 else: # If there already exists an edge parallel to this direction, then use a 1D projection
173 Apply1DProjMesh(MacObject.Mesh[0],Edge,Reference[i])
176 MacObject.Mesh[0].Compute() # Generates the mesh
178 MacObject.DirectionalMeshParams = [MacObject.MeshPar[0]*ReducedRatio[1],MacObject.MeshPar[0]*ReducedRatio[1],MacObject.MeshPar[0]*ReducedRatio[0],MacObject.MeshPar[0]*ReducedRatio[0]]
181 Config.ListObj.append(MacObject)
184 ##########################################################################################################
186 def CompBoxF (MacObject):
187 if Config.debug : print("Generating composite box")
189 dummy1 = geompy.MakeScaleAlongAxes( ElemBox11 (), None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
190 RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
192 MacObject.GeoChildren.append(RectFace)
193 MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
195 if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
198 MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh
199 Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis
201 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
203 #ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1])
206 Vec = [(1,0,0),(0,1,0)]
209 if IsParallel(Edge,Vec[i]):
210 if not Reference[i]: # If this is the first found edge to be parallel to this direction, apply user preferences for meshing
212 ApplyConstant1DMesh(MacObject.Mesh[0],Edge,int(round(MacObject.MeshPar[0][i])))
214 else: # If there already exists an edge parallel to this direction, then use a 1D projection
215 Apply1DProjMesh(MacObject.Mesh[0],Edge,Reference[i])
218 MacObject.Mesh[0].Compute() # Generates the mesh
220 MacObject.DirectionalMeshParams = [MacObject.MeshPar[0][1],MacObject.MeshPar[0][1],MacObject.MeshPar[0][0],MacObject.MeshPar[0][0]]
223 Config.ListObj.append(MacObject)
225 ##########################################################################################################
227 def NonOrtho (MacObject):
228 if Config.debug : print("Generating Non-orthogonal quadrangle")
230 RectFace = Quadrangler (MacObject.PtCoor)
232 MacObject.GeoChildren.append(RectFace)
233 MacObject.GeoChildrenNames.append("Quad_"+ str(len(Config.ListObj)+1))
236 if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
239 MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh
240 Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis
242 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
244 #ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1])
246 Vec = [MacObject.DirVectors(i) for i in range(4)]
248 Dir = [IsParallel(Edge,Vec[j]) for j in range(4)].index(True)
249 DirConv = [0,0,1,1][Dir]
250 ApplyConstant1DMesh(MacObject.Mesh[0],Edge,int(round(MacObject.MeshPar[0][DirConv])))
252 MacObject.Mesh[0].Compute() # Generates the mesh
254 MacObject.DirectionalMeshParams = [MacObject.MeshPar[0][1],MacObject.MeshPar[0][1],MacObject.MeshPar[0][0],MacObject.MeshPar[0][0]]
257 Config.ListObj.append(MacObject)
260 ##########################################################################################################
262 def QuartCyl (MacObject):
263 if Config.debug : print("Generating quarter cylinder")
264 Z_Axis = geompy.MakeVectorDXDYDZ(0., 0., 1.)
265 RotAngle = {'NE' : lambda : 0,
266 'NW' : lambda : math.pi/2,
267 'SW' : lambda : math.pi,
268 'SE' : lambda : -math.pi/2, }[MacObject.MeshPar[1]]()
269 dummy0 = geompy.MakeRotation( ElemQuartCyl(MacObject.MeshPar[2]) , Z_Axis, RotAngle )
270 dummy1 = geompy.MakeScaleAlongAxes( dummy0, None , MacObject.GeoPar[1][0]/10., MacObject.GeoPar[1][1]/10., 1)
271 RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
273 MacObject.GeoChildren.append(RectFace)
274 MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
276 if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
279 MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh
280 Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis
282 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
283 Reg1D = MacObject.Mesh[0].Segment()
285 #if MacObject.MeshPar[0] == 2 and MacObject.MeshPar[2] <= 2.:
286 # 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.")
287 # MacObject.MeshPar[0] = 3
289 Reg1D.NumberOfSegments(MacObject.MeshPar[0])
291 MacObject.Mesh[0].Compute() # Generates the mesh
295 x = MacObject.MeshPar[0]
296 N = QuarCylParam(MacObject.MeshPar[2])+1
298 MacObject.DirectionalMeshParams = {'NE' : lambda : [2*x, N*x, 2*x, N*x],
299 'NW' : lambda : [N*x, 2*x, 2*x, N*x],
300 'SW' : lambda : [N*x, 2*x, N*x, 2*x],
301 'SE' : lambda : [2*x, N*x, N*x, 2*x], }[MacObject.MeshPar[1]]()
303 Config.ListObj.append(MacObject)
306 ##########################################################################################################
307 # Below this are the elementary calculation/visualization functions
308 ##########################################################################################################
310 def Publish (ObjToPublish,NamesToPublish):
312 for GeoObj in ObjToPublish :
313 geompy.addToStudy(GeoObj,NamesToPublish[i])
316 def IsParallel (Edge, Vector):
318 Function checks whether a given edge object is parallel to a reference vector.
319 Output can be 0 (not parallel) or 1 (parallel and same sense) or 2 (parallel and opposite sense).
320 If the reference vector is null, the function returns 0
322 if Vector == (0,0,0) : return 0
324 P1 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,0))
325 P2 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,1))
326 V0 = [ P1[0] - P2[0], P1[1] - P2[1], P1[2] - P2[2] ]
327 if Distance2Pt((0,0,0),CrossProd(V0,Vector))<1e-7 and DotProd(V0,Vector) > 0 : return 1
328 elif Distance2Pt((0,0,0),CrossProd(V0,Vector))<1e-7 and DotProd(V0,Vector) < 0 : return 2
331 def IsOnCircle (Edge, Center, Radius):
333 Function checks whether a given edge object belong to the periphery of a circle defined by its
335 Output can be 0 (does not belong) or 1 (belongs).
336 If the reference Radius is null, the function returns 0
337 Note that this function is basic in the sense that it only checks if the two border points of a
338 given edge belong to the arc of reference.
340 if Radius == 0 : return 0
342 P1 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,0))
343 P2 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,1))
344 if abs(Distance2Pt(Center,P1)-Radius) < 1e-6 and abs(Distance2Pt(Center,P2)-Radius) < 1e-6:
349 def CrossProd(V1,V2):
351 Determines the cross product of two 3D vectors
353 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]])
355 def QuarCylParam(PitchRatio):
356 R = float(PitchRatio)/(PitchRatio+1)
358 X = (R+Eps/2.)*math.sin(math.pi/4)+Eps/2.
359 N = int(math.floor((math.pi*R/4.)/(Eps/2.)))
364 Determines the dot product of two 3D vectors
366 if len(V1)==2 : V1.append(0)
367 if len(V2)==2 : V2.append(0)
369 return (V1[0]*V2[0]+V1[1]*V2[1]+V1[2]*V2[2])
371 def Distance2Pt(P1,P2):
373 Returns the distance between two points
375 return (math.sqrt((P1[0]-P2[0])**2+(P1[1]-P2[1])**2+(P1[2]-P2[2])**2))
377 def ApplyConstant1DMesh (ParentMsh, Edge, Nseg):
378 Reg1D = ParentMsh.Segment(geom=Edge)
379 Len = Reg1D.NumberOfSegments(Nseg)
381 def Apply1DProjMesh (ParentMsh, Edge, Ref):
382 Proj1D = ParentMsh.Projection1D(geom=Edge)
383 SrcEdge = Proj1D.SourceEdge(Ref,None,None,None)
385 def EdgeLength (Edge):
387 This function returns the edge object length.
389 P1 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,0))
390 P2 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,1))
391 return Distance2Pt(P1,P2)
395 return Angle*math.pi/180
398 return Angle*180/math.pi
400 def F2D (FloatNumber):
401 return round(FloatNumber*100.)/100.
403 def BezierGen (PointA, PointB, AngleA, AngleB):
405 if AngleA == 0 and AngleB == 0 : return (geompy.MakeEdge(PointA, PointB))
407 A = geompy.PointCoordinates(PointA)
408 B = geompy.PointCoordinates(PointB)
409 dAB = Distance2Pt(A,B)
410 dAC = dAB * (math.tan(AngleA)*math.tan(AngleB)) / (math.sin(AngleA) * ( math.tan(AngleA)+math.tan(AngleB) ) )
411 AngleOX_AB = math.acos((B[0]-A[0])/dAB)
412 PointC = geompy.MakeVertex(A[0]+math.cos(AngleA+AngleOX_AB)*dAC,A[1]+math.sin(AngleA+AngleOX_AB)*dAC,0)
413 CurveACB = geompy.MakeBezier([PointA,PointC,PointB])
416 def GetSideAngleForBezier (PointA , PointB):
418 This function takes for input two points A and B where the bezier line is needed. It calculates the incident
419 angle needed at point A so that the final curve is either at 0 or 90 degrees from the x'Ox axis
421 A = geompy.PointCoordinates(PointA)
422 B = geompy.PointCoordinates(PointB)
424 dAB = Distance2Pt(A,B)
425 Alpha = math.acos(ABx/dAB)
426 #print "New angle request"
427 #print ABx, dAB, R2D(Alpha)
428 if Alpha < math.pi/4 :
429 #print "returning", R2D(-Alpha)
431 elif Alpha < 3*math.pi/4 :
432 #print "returning", R2D(-(Alpha-math.pi/2))
433 return -(Alpha-math.pi/2)
435 #print "returning", R2D(-(Alpha-math.pi))
436 return -(Alpha-math.pi)
438 def VecDivRatio (Vec1, Vec2):
440 This function tries to find the ratio of Vec1 on Vec2 while neglecting any zero term in Vec1. This is used afterwards
441 for determining the global mesh parameter from automatically detected directional mesh params. If no compatibility is
442 possible, the function returns -1
445 for i in range(len(Vec1)) :
446 Vec3.append(float(Vec1[i])/Vec2[i])
449 if not (abs(i)<1e-7) : Ratio.append(i)
451 if min(Ratio) == max(Ratio) and min(Ratio)==int(min(Ratio)) : return(min(Ratio))
457 def ReduceRatio (dx, dy):
459 This function transforms a decimal ratio into a scale between two integers, for example : [0.2,0.05] --> [4,1] ;
463 if isinteger(ratio) : return [1,ratio]
464 elif dx == 1 : # when this function is called recursively!
465 for i in range(1,20) : # searches over 20 decimals
466 if isinteger(ratio * (10**i) ) :
467 Output = GetScale((10**i),int(round(ratio * (10**i) ) ) )
470 for n in range(0,i) :
471 if isinteger(ratio * ( 10**(i)-10**(n) )) :
472 Output = GetScale( 10**(i)-10**(n) , int(round(ratio * ( 10**(i)-10**(n) ) ) ) )
474 if not (Output==[0,0]) : break
477 for i in range(1,10) : # searches over 10 decimals
478 if isinteger(ratio * (10**i) ) :
479 Output = GetScale((10**i),int(round(ratio * (10**i) ) ) )
482 for n in range(0,i) :
483 if isinteger(ratio * ( 10**(i)-10**(n) )) :
484 Output = GetScale( 10**(i)-10**(n) , int(round(ratio * ( 10**(i)-10**(n) ) ) ) )
486 if not (Output==[0,0]) : break
489 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...")
491 A = ReduceRatio (dx, dy-dx)
492 return ([A[0],A[1]+A[0]])
494 A = ReduceRatio (dy, dx-dy)
495 return ([A[1]+A[0],A[0]])
501 This function is called within ReduceRatio and aims to reduce down two integers X and Y by dividing them with their common divisors;
502 Example: 25 and 5 ---> 5 and 1 / 63 and 12 ---> 21 and 4
505 Divisor = 2 # Initializing the divisor
506 while MaxDiv >= Divisor :
511 MaxDiv = max(MaxDiv,X0)
514 MaxDiv = max(MaxDiv,Y0)
519 Divisor = Divisor + 1
524 This functions applies a simple check if the entered value is an integer
526 x = float('%.5f' % (x)) #Truncate x to 5 digits after the decimal point
527 if math.ceil(x) == math.floor(x) : return True
529 ##########################################################################################
530 # Below this are the functions that create the elementary forms for the macro objects
531 ##########################################################################################
535 This function returns a simple square face of 1 side length
537 RectFace = geompy.MakeFaceHW(1, 1, 1)
542 This function returns a square face of 1 side length, partitioned
543 according to the elementary 4 to 2 reductor method
545 OrigRectFace = geompy.MakeFaceHW(1, 1, 1)
547 SouthPt1 = geompy.MakeVertex (-.25, -.5, 0)
548 SouthPt2 = geompy.MakeVertex (0, -.5, 0)
549 SouthPt3 = geompy.MakeVertex (.25, -.5, 0)
550 WestPt1 = geompy.MakeVertex (-.5, -.5+1./3, 0)
551 WestPt2 = geompy.MakeVertex (-.5, -.5+2./3, 0)
552 EastPt1 = geompy.MakeVertex (.5, -.5+1./3, 0)
553 EastPt2 = geompy.MakeVertex (.5, -.5+2./3, 0)
554 NorthPt = geompy.MakeVertex (0, .5, 0)
555 MidPt1 = geompy.MakeVertex (0, .05, 0)
556 MidPt2 = geompy.MakeVertex (.2, -.18, 0)
557 MidPt3 = geompy.MakeVertex (0, -.28, 0)
558 MidPt4 = geompy.MakeVertex (-.2, -.18, 0)
561 Cutter.append(geompy.MakeEdge(SouthPt2, MidPt3))
562 Cutter.append(geompy.MakeEdge(MidPt1, NorthPt))
563 Cutter.append(BezierGen(SouthPt1, MidPt4, GetSideAngleForBezier(SouthPt1,MidPt4), D2R(15)))
564 Cutter.append(BezierGen(SouthPt3, MidPt2, GetSideAngleForBezier(SouthPt3,MidPt2), D2R(-15)))
565 Cutter.append(BezierGen(WestPt1, MidPt4, GetSideAngleForBezier(WestPt1,MidPt4), D2R(-10)))
566 Cutter.append(BezierGen(EastPt1, MidPt2, GetSideAngleForBezier(EastPt1,MidPt2), D2R(10)))
567 Cutter.append(BezierGen(WestPt2, MidPt1, GetSideAngleForBezier(WestPt2,MidPt1), D2R(-10)))
568 Cutter.append(BezierGen(EastPt2, MidPt1, GetSideAngleForBezier(EastPt2,MidPt1), D2R(10)))
569 Cutter.append(BezierGen(MidPt2, MidPt1, D2R(-15), D2R(-15)))
570 Cutter.append(BezierGen(MidPt3, MidPt2, D2R(10), D2R(15)))
571 Cutter.append(BezierGen(MidPt3, MidPt4, D2R(-10), D2R(-15)))
572 Cutter.append(BezierGen(MidPt4, MidPt1, D2R(15), D2R(15)))
574 RectFace = geompy.MakePartition([OrigRectFace],Cutter, [], [],4, 0, [], 0) #Creating the partition object
576 #for SingleCut in Cutter :
577 # geompy.addToStudy(SingleCut,'Cutter'+str(i))
579 #geompy.addToStudy(RectFace,'RectFace')
584 This function returns a square face of 1 side length, partitioned
585 according to the elementary edge with 3 to 2 reductor
587 OrigRectFace = geompy.MakeFaceHW(1., 1., 1)
589 SouthPt1 = geompy.MakeVertex (-1./6, -0.5, 0.)
590 SouthPt2 = geompy.MakeVertex ( 1./6, -0.5, 0.)
591 WestPt1 = geompy.MakeVertex (-0.5, -1./6, 0.)
592 WestPt2 = geompy.MakeVertex (-0.5, 1./6, 0.)
593 EastPt = geompy.MakeVertex ( 0.5, 0., 0.)
594 NorthPt = geompy.MakeVertex (0., 0.5, 0.)
596 MidPt1 = geompy.MakeVertex (-0.2, -0.2, 0.)
597 MidPt2 = geompy.MakeVertex ( -0.02, -0.02, 0.)
600 Cutter.append(BezierGen(SouthPt1, MidPt1, GetSideAngleForBezier(SouthPt1,MidPt1) , D2R(-5)))
601 Cutter.append(BezierGen( WestPt1, MidPt1, GetSideAngleForBezier(WestPt1 ,MidPt1) , D2R(-5)))
602 Cutter.append(BezierGen(SouthPt2, MidPt2, GetSideAngleForBezier(SouthPt2,MidPt2) , D2R(-10)))
603 Cutter.append(BezierGen( EastPt, MidPt2, GetSideAngleForBezier(EastPt ,MidPt2) , D2R(5)))
604 Cutter.append(BezierGen( WestPt2, MidPt2, GetSideAngleForBezier(WestPt2 ,MidPt2) , D2R(-10)))
605 Cutter.append(BezierGen( MidPt2, NorthPt, GetSideAngleForBezier(NorthPt ,MidPt2) , D2R(-5)))
607 Cutter.append(geompy.MakeEdge(MidPt1, MidPt2))
609 RectFace = geompy.MakePartition([OrigRectFace],Cutter, [], [],4, 0, [], 0) #Creating the partition object
611 #for SingleCut in Cutter :
612 # geompy.addToStudy(SingleCut,'Cutter'+str(i))
614 #geompy.addToStudy(RectFace,'RectFace')
617 def Quadrangler (Points):
619 This function returns a quadranglar face based on four points, non of which 3 are non-colinear.
620 The points are defined by their 2D [(x1,y1),(x2,y2)..] coordinates.
621 Note that the list of points is already arranged upon the creation in MacObject
624 for Point in Points: Pt.append(geompy.MakeVertex(Point[0], Point[1], 0))
625 # The first point is added at the end of the list in order to facilitate the line creation
627 #Draw the lines in order to form the 4 side polygon
629 for i in range(4) : Ln.append(geompy.MakeLineTwoPnt(Pt[i],Pt[i+1]))
630 RectFace = geompy.MakeQuad (Ln[0],Ln[1],Ln[2],Ln[3])
635 This function returns a quarter cylinder to box relay of 1 side length, partitioned
636 with a pitch ratio of K, In other words the side of the box is R*(1+(1/K))
638 R = 10.*float(K)/(K+1)
641 Config.theStudy.SetReal("R" , R)
642 Config.theStudy.SetReal("minusR" , -R)
643 Config.theStudy.SetReal("Eps", Eps)
645 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])
646 CylFace = geompy.MakeFace(CylWire, 1)
648 SouthPt = geompy.MakeVertex (R+Eps/2., 0., 0)
649 SouthWestPt = geompy.MakeVertex ( 0.,0., 0) #The origin can be used for practical partionning objectifs
650 WestPt = geompy.MakeVertex (0., R+Eps/2., 0)
652 N = int(math.floor((math.pi*R/4.)/(Eps/2.)))
653 X = 10.*(1.-1./(N+1))
656 EastPt = geompy.MakeVertex (10.0, X, 0.)
657 NorthPt = geompy.MakeVertex ( X, 10.0, 0.)
659 DivFactor = 8./(F2D(math.log(K))-0.223)
660 #MidPt = geompy.MakeVertex ((R+Eps)*math.cos(math.pi/4), (R+Eps)*math.sin(math.pi/4), 0.)
661 MidPt = geompy.MakeVertex (X-Eps/DivFactor, X-Eps/DivFactor, 0.)
664 Cutter.append(BezierGen(SouthWestPt, MidPt, GetSideAngleForBezier(SouthWestPt,MidPt) , D2R(-5)))
665 Cutter.append(BezierGen( EastPt, MidPt, GetSideAngleForBezier(EastPt,MidPt) , D2R(5)))
666 Cutter.append(BezierGen( MidPt, NorthPt, (-1)**((K<1.25)*1)*D2R(-5), GetSideAngleForBezier(NorthPt,MidPt)))
667 SMBezier = BezierGen( SouthPt, MidPt, GetSideAngleForBezier(SouthPt ,MidPt) , D2R((K<1.25)*180-5))
668 WMBezier = BezierGen( WestPt, MidPt, GetSideAngleForBezier(WestPt, MidPt) , D2R(-5))
669 Cutter.append(WMBezier)
670 Cutter.append(SMBezier)
672 for i in range(1,N) :
673 # Determining intermediate points on the bezier lines and then performing additional cuts
675 TempAnglePlus = (math.pi/4)*(1+float(i)/N)
676 SectionResult = CutnGroup.Go(WMBezier, [(0,0,0,math.sin(TempAnglePlus),-math.cos(TempAnglePlus),0)], [1], ['Dummy'], 0)
677 TempPt1 = SectionResult[1][0]
678 TempPt11 = geompy.MakeVertex ((N-i)*X/N, 10., 0)
680 TempAngleMinus = (math.pi/4)*(1-float(i)/N)
681 SectionResult = CutnGroup.Go(SMBezier, [(0,0,0,math.sin(TempAngleMinus),-math.cos(TempAngleMinus),0)], [1], ['Dummy'], 0)
682 TempPt2 = SectionResult[1][0]
683 TempPt21 = geompy.MakeVertex (10., (N-i)*X/N, 0)
685 Cutter.append(geompy.MakeEdge(SouthWestPt, TempPt1))
686 Cutter.append(geompy.MakeEdge(SouthWestPt, TempPt2))
687 Cutter.append(geompy.MakeEdge(TempPt1, TempPt11))
688 Cutter.append(geompy.MakeEdge(TempPt2, TempPt21))
690 CylFace = geompy.MakePartition([CylFace],Cutter, [], [],4, 0, [], 0) #Creating the partition object
691 CylFace = geompy.MakeTranslation(CylFace, -5., -5., 0.0)
695 def CompatibilityTest(MacObject):
696 Type = MacObject.Type
698 BaseDirPar = [1,1,1,1]
699 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
700 elif Type == 'Box42' :
701 BaseDirPar = {'SN' : lambda : [3, 3, 4, 2],
702 'NS' : lambda : [3, 3, 2, 4],
703 'EW' : lambda : [2, 4, 3, 3],
704 'WE' : lambda : [4, 2, 3, 3], }[MacObject.MeshPar[1]]()
705 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
706 elif Type == 'BoxAng32' :
707 BaseDirPar = {'NE' : lambda : [3, 2, 3, 2],
708 'NW' : lambda : [2, 3, 3, 2],
709 'SW' : lambda : [2, 3, 2, 3],
710 'SE' : lambda : [3, 2, 2, 3], }[MacObject.MeshPar[1]]()
711 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
712 elif Type == 'CompBox' :
713 #print "dx is: ", MacObject.GeoPar[1][1], ". dy is: ",MacObject.GeoPar[1][0]
714 ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0], MacObject.GeoPar[1][1])
716 BaseDirPar = [ReducedRatio[1], ReducedRatio[1], ReducedRatio[0], ReducedRatio[0]]
717 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
719 elif Type == 'QuartCyl' :
720 N = QuarCylParam(MacObject.MeshPar[2])+1
721 BaseDirPar = {'NE' : lambda : [2, N, 2, N],
722 'NW' : lambda : [N, 2, 2, N],
723 'SW' : lambda : [N, 2, N, 2],
724 'SE' : lambda : [2, N, N, 2], }[MacObject.MeshPar[1]]()
725 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
726 elif Type == 'CompBoxF' :
727 RealRatio = MacObject.GeoPar[1][1]/MacObject.GeoPar[1][0]
730 if MacObject.DirectionalMeshParams[2]+MacObject.DirectionalMeshParams[3] :
731 A = int(max(MacObject.DirectionalMeshParams[2:4]))
732 Xd = int(VecDivRatio([A,0,0,0], [1,1,1,1]))
733 if MacObject.DirectionalMeshParams[0]+MacObject.DirectionalMeshParams[1] :
734 A = int(max(MacObject.DirectionalMeshParams[0:2]))
735 Yd = int(VecDivRatio([0,0,A,0], [1,1,1,1]))
737 if Xd == 0 and Yd : Xd = int(round(Yd/RealRatio))
738 elif Yd == 0 : Yd = int(round(RealRatio*Xd))
741 elif Type == 'NonOrtho' :
742 MeanDX = 0.5*(IntLen(MacObject.DirBoundaries(0))+IntLen(MacObject.DirBoundaries(1)))
743 MeanDY = 0.5*(IntLen(MacObject.DirBoundaries(2))+IntLen(MacObject.DirBoundaries(3)))
744 RealRatio = MeanDY/MeanDX
747 if MacObject.DirectionalMeshParams[2]+MacObject.DirectionalMeshParams[3] :
748 A = int(max(MacObject.DirectionalMeshParams[2:4]))
749 Xd = int(VecDivRatio([A,0,0,0], [1,1,1,1]))
750 if MacObject.DirectionalMeshParams[0]+MacObject.DirectionalMeshParams[1] :
751 A = int(max(MacObject.DirectionalMeshParams[0:2]))
752 Yd = int(VecDivRatio([0,0,A,0], [1,1,1,1]))
754 if Xd == 0 and Yd : Xd = int(round(Yd/RealRatio))
755 elif Yd == 0 : Yd = int(round(RealRatio*Xd))
759 def IntLen (Interval) :
761 This function returns the length of a given interval even if the latter is not sorted correctly.
763 return abs(Interval[1]-Interval[0])
765 def NextTo (RefBox, Direction, Extension):
767 This functions returns geometrical parameters for easy positioning of neighbouring objects.
768 The input (RefBox) and output are in the form : [(X0,Y0),(DX,DY)]
775 DirectionalCoef = {'Above' : lambda : [ 0, 1],
776 'Below' : lambda : [ 0,-1],
777 'Right' : lambda : [ 1, 0],
778 'Left ' : lambda : [-1, 0], }[Direction]()
780 X0_1 = X0_0+ DirectionalCoef[0] * (DX_0/2.+Extension/2.)
781 DX_1 = abs(DirectionalCoef[0]) * (Extension) + abs(DirectionalCoef[1])*DX_0
782 Y0_1 = Y0_0+ DirectionalCoef[1] * (DY_0/2.+Extension/2.)
783 DY_1 = abs(DirectionalCoef[1]) * (Extension) + abs(DirectionalCoef[0])*DY_0
785 return [(X0_1,Y0_1),(DX_1,DY_1)]
787 def GeomMinMax (PtA, PtB):
789 This function returns geometrical parameters in the format [(X0,Y0),(DX,DY)]. The input being
790 the coordinates of two points (Xa,Ya), (Xb,Yb).
792 # First test that the vector relying the two points is oblique
793 AB = [PtB[0]- PtA[0],PtB[1]- PtA[1]]
795 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")
798 X0 = 0.5*(PtA[0]+PtB[0])
799 Y0 = 0.5*(PtA[1]+PtB[1])
802 return [(X0,Y0),(DX,DY)]
804 def AddIfDifferent (List, Element):
805 if not(Element in List):
806 List = List+(Element,)
809 def IndexMultiOcc (Array,Element) :
811 This functions returns the occurrences indices of Element in Array.
812 As opposed to Array.index(Element) method, this allows determining
813 multiple entries rather than just the first one!
816 try : Array.index(Element)
817 except ValueError : print("No more occurrences")
818 else : Output.append(Array.index(Element))
820 if not(Output == []) and len(Array) > 1 :
821 for index, ArrElem in enumerate(Array[Output[0]+1:]) :
822 if ArrElem == Element : Output.append(index+Output[0]+1)
826 def SortList (ValList, CritList):
828 SortedCritList = sorted(copy.copy(CritList))
829 for i in range(0,len(ValList)):
831 if not(SortedCritList[i]==SortedCritList[i-1]):
832 index = IndexMultiOcc(CritList,SortedCritList[i])
833 Output= Output + [ValList[j] for j in index]
835 index = IndexMultiOcc(CritList,SortedCritList[i])
836 Output= Output + [ValList[j] for j in index]
840 def SortPoints(Points):
842 This function sorts a list of the coordinates of N points as to start at
843 an origin that represents Xmin and Xmax and then proceed in a counter
847 Xmin = min([Points[i][0] for i in range(NbPts)])
848 Ymin = min([Points[i][1] for i in range(NbPts)])
849 Xmax = max([Points[i][0] for i in range(NbPts)])
850 Ymax = max([Points[i][1] for i in range(NbPts)])
851 Crit = [(abs(Point[0]-Xmin)+0.1*(Xmax-Xmin))*(abs(Point[1]-Ymin)+0.1*(Ymax-Ymin)) for Point in Points]
852 #print "Input Points : ", Points
853 #print "Sorting Criterion : ", Crit
854 Order = SortList (list(range(NbPts)), Crit)
855 #print "Sorted Results : ", Order
857 Output.append(Points[Order[0]])
859 Point0 = Points[Order[0]]
860 #print "Reference point :", Point0
862 V = [[Point1[0]-Point0[0],Point1[1]-Point0[1]] for Point1 in Points]
863 Cosines = [-(vec[0]-1E-10)/(math.sqrt(DotProd(vec,vec)+1e-25)) for vec in V]
864 #print "Cosines criterion :", Cosines
865 Order = SortList(list(range(NbPts)),Cosines)
866 #print "Ordered points:", Order
867 for PtIndex in Order[:-1]: Output.append(Points[PtIndex])