src/Tools/MacMesh/MacMesh/GenFunctions.py

   1 # In this file are all the generation functions for manipulating the different created macro-objects
   2 import geompy, smesh
   3 import math, copy
   4 import Config
   5 import CutnGroup
   6 import CompositeBox
   7
   8 ##########################################################################################################
   9
  10 def Box11 (MacObject):
  11         if Config.debug : print "Generating regular box"
  12
  13         dummy1 = geompy.MakeScaleAlongAxes( ElemBox11 (), None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
  14         RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
  15
  16         MacObject.GeoChildren.append(RectFace)
  17         MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
  18
  19         if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
  20
  21         if Config.publish :
  22                 MacObject.Mesh.append(smesh.Mesh(RectFace))                     # Creation of a new mesh
  23                 Quad2D = MacObject.Mesh[0].Quadrangle()                 # Applying a quadrangle hypothesis
  24
  25                 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6)  # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
  26                 Reg1D = MacObject.Mesh[0].Segment()
  27                 Reg1D.NumberOfSegments(MacObject.MeshPar[0])
  28
  29                 MacObject.Mesh[0].Compute()                     # Generates the mesh
  30
  31         MacObject.DirectionalMeshParams = [MacObject.MeshPar[0],MacObject.MeshPar[0],MacObject.MeshPar[0],MacObject.MeshPar[0]]
  32
  33         MacObject.status = 1
  34         Config.ListObj.append(MacObject)
  35         return MacObject
  36
  37 ##########################################################################################################
  38
  39 def Box42 (MacObject):
  40         if Config.debug : print "Generating box 4-2 reducer"
  41
  42         Z_Axis = geompy.MakeVectorDXDYDZ(0., 0., 1.)
  43         RotAngle = {'SN' : lambda : 0,
  44                     'NS' : lambda : math.pi,
  45                     'EW' : lambda : math.pi/2,
  46                     'WE' : lambda : -math.pi/2, }[MacObject.MeshPar[1]]()
  47         dummy0 = geompy.MakeRotation( ElemBox42 () , Z_Axis, RotAngle )
  48         dummy1 = geompy.MakeScaleAlongAxes( dummy0, None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
  49         RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
  50
  51         MacObject.GeoChildren.append(RectFace)
  52         MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
  53
  54         if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
  55
  56         if Config.publish :
  57                 MacObject.Mesh.append(smesh.Mesh(RectFace))                     # Creation of a new mesh
  58                 Quad2D = MacObject.Mesh[0].Quadrangle()                 # Applying a quadrangle hypothesis
  59
  60                 EdgeIDs = geompy.SubShapeAllSorted(RectFace,6)  # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
  61                 Reg1D = MacObject.Mesh[0].Segment()
  62                 Reg1D.NumberOfSegments(MacObject.MeshPar[0])
  63
  64                 MacObject.Mesh[0].Compute()                     # Generates the mesh
  65
  66         MacObject.status = 1
  67
  68         x = MacObject.MeshPar[0]
  69         MacObject.DirectionalMeshParams = {'SN' : lambda : [3*x, 3*x, 4*x, 2*x],
  70                                            'NS' : lambda : [3*x, 3*x, 2*x, 4*x],
  71                                            'EW' : lambda : [2*x, 4*x, 3*x, 3*x],
  72                                            'WE' : lambda : [4*x, 2*x, 3*x, 3*x], }[MacObject.MeshPar[1]]()
  73
  74         Config.ListObj.append(MacObject)
  75         return MacObject
  76
  77
  78 ##########################################################################################################
  79
  80 def BoxAng32 (MacObject):
  81         if Config.debug : print "Generating sharp angle"
  82         Z_Axis = geompy.MakeVectorDXDYDZ(0., 0., 1.)
  83         RotAngle = {'NE' : lambda : 0,
  84                     'NW' : lambda : math.pi/2,
  85                     'SW' : lambda : math.pi,
  86                     'SE' : lambda : -math.pi/2, }[MacObject.MeshPar[1]]()
  87         dummy0 = geompy.MakeRotation( ElemEdge32 () , Z_Axis, RotAngle )
  88         dummy1 = geompy.MakeScaleAlongAxes( dummy0, None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
  89         RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
  90
  91         MacObject.GeoChildren.append(RectFace)
  92         MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
  93
  94         if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
  95
  96         if Config.publish :
  97                MacObject.Mesh.append(smesh.Mesh(RectFace))                     # Creation of a new mesh
  98                Quad2D = MacObject.Mesh[0].Quadrangle()                 # Applying a quadrangle hypothesis
  99
 100                EdgeIDs = geompy.SubShapeAllSorted(RectFace,6)      # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
 101                Reg1D = MacObject.Mesh[0].Segment()
 102                Reg1D.NumberOfSegments(MacObject.MeshPar[0])
 103
 104                MacObject.Mesh[0].Compute()                     # Generates the mesh
 105
 106         MacObject.status = 1
 107
 108         x = MacObject.MeshPar[0]
 109         MacObject.DirectionalMeshParams = {'NE' : lambda : [3*x, 2*x, 3*x, 2*x],
 110                                            'NW' : lambda : [2*x, 3*x, 3*x, 2*x],
 111                                            'SW' : lambda : [2*x, 3*x, 2*x, 3*x],
 112                                            'SE' : lambda : [3*x, 2*x, 2*x, 3*x], }[MacObject.MeshPar[1]]()
 113
 114         Config.ListObj.append(MacObject)
 115         return MacObject
 116 ##########################################################################################################
 117 def CompBox (MacObject):
 118         if Config.debug : print "Generating composite box"
 119
 120         dummy1 = geompy.MakeScaleAlongAxes( ElemBox11 (), None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
 121         RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
 122
 123         MacObject.GeoChildren.append(RectFace)
 124         MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
 125
 126         if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
 127
 128         if Config.publish :
 129               MacObject.Mesh.append(smesh.Mesh(RectFace))                     # Creation of a new mesh
 130               Quad2D = MacObject.Mesh[0].Quadrangle()                 # Applying a quadrangle hypothesis
 131
 132               EdgeIDs = geompy.SubShapeAllSorted(RectFace,6)      # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
 133
 134               ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1])
 135
 136               Reference = [0,0,0]
 137               Vec = [(1,0,0),(0,1,0)]
 138               for Edge in EdgeIDs:
 139                               for i in range(0,2):
 140                                       if IsParallel(Edge,Vec[i]):
 141                                               if not Reference[i]:    # If this is the first found edge to be parallel to this direction, apply user preferences for meshing
 142                                                       Reference[i] = Edge
 143                                                       ApplyConstant1DMesh(MacObject.Mesh[0],Edge,int(round(ReducedRatio[i]*MacObject.MeshPar[0])))
 144                                                       break
 145                                               else:                   # If there already exists an edge parallel to this direction, then use a 1D projection
 146                                                       Apply1DProjMesh(MacObject.Mesh[0],Edge,Reference[i])
 147                                                       break
 148
 149               MacObject.Mesh[0].Compute()                     # Generates the mesh
 150
 151         MacObject.DirectionalMeshParams = [MacObject.MeshPar[0]*ReducedRatio[1],MacObject.MeshPar[0]*ReducedRatio[1],MacObject.MeshPar[0]*ReducedRatio[0],MacObject.MeshPar[0]*ReducedRatio[0]]
 152
 153         MacObject.status = 1
 154         Config.ListObj.append(MacObject)
 155         return MacObject
 156
 157 ##########################################################################################################
 158
 159 def CompBoxF (MacObject):
 160         if Config.debug : print "Generating composite box"
 161
 162         dummy1 = geompy.MakeScaleAlongAxes( ElemBox11 (), None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1)
 163         RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
 164
 165         MacObject.GeoChildren.append(RectFace)
 166         MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
 167
 168         if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
 169
 170         if Config.publish :
 171                MacObject.Mesh.append(smesh.Mesh(RectFace))                     # Creation of a new mesh
 172                Quad2D = MacObject.Mesh[0].Quadrangle()                 # Applying a quadrangle hypothesis
 173
 174                EdgeIDs = geompy.SubShapeAllSorted(RectFace,6)      # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
 175
 176                #ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1])
 177
 178                Reference = [0,0,0]
 179                Vec = [(1,0,0),(0,1,0)]
 180                for Edge in EdgeIDs:
 181                                for i in range(0,2):
 182                                        if IsParallel(Edge,Vec[i]):
 183                                                if not Reference[i]:    # If this is the first found edge to be parallel to this direction, apply user preferences for meshing
 184                                                        Reference[i] = Edge
 185                                                        ApplyConstant1DMesh(MacObject.Mesh[0],Edge,int(round(MacObject.MeshPar[0][i])))
 186                                                        break
 187                                                else:                   # If there already exists an edge parallel to this direction, then use a 1D projection
 188                                                        Apply1DProjMesh(MacObject.Mesh[0],Edge,Reference[i])
 189                                                        break
 190
 191                MacObject.Mesh[0].Compute()                 # Generates the mesh
 192
 193         MacObject.DirectionalMeshParams = [MacObject.MeshPar[0][1],MacObject.MeshPar[0][1],MacObject.MeshPar[0][0],MacObject.MeshPar[0][0]]
 194
 195         MacObject.status = 1
 196         Config.ListObj.append(MacObject)
 197         return MacObject
 198 ##########################################################################################################
 199
 200 def NonOrtho (MacObject):
 201         if Config.debug : print "Generating Non-orthogonal quadrangle"
 202
 203         RectFace = Quadrangler (MacObject.PtCoor)
 204
 205         MacObject.GeoChildren.append(RectFace)
 206         MacObject.GeoChildrenNames.append("Quad_"+ str(len(Config.ListObj)+1))
 207
 208
 209         if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
 210
 211         if Config.publish :
 212                MacObject.Mesh.append(smesh.Mesh(RectFace))                     # Creation of a new mesh
 213                Quad2D = MacObject.Mesh[0].Quadrangle()                 # Applying a quadrangle hypothesis
 214
 215                EdgeIDs = geompy.SubShapeAllSorted(RectFace,6)      # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
 216
 217                #ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1])
 218
 219                Vec = [MacObject.DirVectors(i) for i in range(4)]
 220                for Edge in EdgeIDs:
 221                         Dir = [IsParallel(Edge,Vec[j]) for j in range(4)].index(True)
 222                         DirConv = [0,0,1,1][Dir]
 223                         ApplyConstant1DMesh(MacObject.Mesh[0],Edge,int(round(MacObject.MeshPar[0][DirConv])))
 224
 225                MacObject.Mesh[0].Compute()                 # Generates the mesh
 226
 227         MacObject.DirectionalMeshParams = [MacObject.MeshPar[0][1],MacObject.MeshPar[0][1],MacObject.MeshPar[0][0],MacObject.MeshPar[0][0]]
 228
 229         MacObject.status = 1
 230         Config.ListObj.append(MacObject)
 231         return MacObject
 232
 233 ##########################################################################################################
 234
 235 def QuartCyl (MacObject):
 236         if Config.debug : print "Generating quarter cylinder"
 237         Z_Axis = geompy.MakeVectorDXDYDZ(0., 0., 1.)
 238         RotAngle = {'NE' : lambda : 0,
 239                     'NW' : lambda : math.pi/2,
 240                     'SW' : lambda : math.pi,
 241                     'SE' : lambda : -math.pi/2, }[MacObject.MeshPar[1]]()
 242         dummy0 = geompy.MakeRotation( ElemQuartCyl(MacObject.MeshPar[2]) , Z_Axis, RotAngle )
 243         dummy1 = geompy.MakeScaleAlongAxes( dummy0, None , MacObject.GeoPar[1][0]/10., MacObject.GeoPar[1][1]/10., 1)
 244         RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0)
 245
 246         MacObject.GeoChildren.append(RectFace)
 247         MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
 248
 249         if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
 250
 251         if Config.publish :
 252                MacObject.Mesh.append(smesh.Mesh(RectFace))                     # Creation of a new mesh
 253                Quad2D = MacObject.Mesh[0].Quadrangle()                 # Applying a quadrangle hypothesis
 254
 255                EdgeIDs = geompy.SubShapeAllSorted(RectFace,6)      # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
 256                Reg1D = MacObject.Mesh[0].Segment()
 257
 258                #if MacObject.MeshPar[0] == 2 and MacObject.MeshPar[2] <= 2.:
 259                #         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.")
 260                #         MacObject.MeshPar[0] = 3
 261
 262                Reg1D.NumberOfSegments(MacObject.MeshPar[0])
 263
 264                MacObject.Mesh[0].Compute()                     # Generates the mesh
 265
 266         MacObject.status = 1
 267
 268         x = MacObject.MeshPar[0]
 269         N = QuarCylParam(MacObject.MeshPar[2])+1
 270
 271         MacObject.DirectionalMeshParams = {'NE' : lambda : [2*x, N*x, 2*x, N*x],
 272                                            'NW' : lambda : [N*x, 2*x, 2*x, N*x],
 273                                            'SW' : lambda : [N*x, 2*x, N*x, 2*x],
 274                                            'SE' : lambda : [2*x, N*x, N*x, 2*x], }[MacObject.MeshPar[1]]()
 275
 276         Config.ListObj.append(MacObject)
 277         return MacObject
 278
 279 ##########################################################################################################
 280 # Below this are the elementary calculation/visualization functions
 281 ##########################################################################################################
 282
 283 def Publish (ObjToPublish,NamesToPublish):
 284         i = 0
 285         for GeoObj in ObjToPublish :
 286                 geompy.addToStudy(GeoObj,NamesToPublish[i])
 287                 i = i+1
 288
 289 def IsParallel (Edge, Vector):
 290         """
 291         Function checks whether a given edge object is parallel to a reference vector.
 292         Output can be 0 (not parallel) or 1 (parallel and same sense) or 2 (parallel and opposite sense).
 293         If the reference vector is null, the function returns 0
 294         """
 295         if Vector == (0,0,0) : return 0
 296         else :
 297                 P1 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,0))
 298                 P2 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,1))
 299                 V0 = [ P1[0] - P2[0], P1[1] - P2[1], P1[2] - P2[2] ]
 300                 if Distance2Pt((0,0,0),CrossProd(V0,Vector))<1e-7 and DotProd(V0,Vector) > 0 : return 1
 301                 elif Distance2Pt((0,0,0),CrossProd(V0,Vector))<1e-7 and DotProd(V0,Vector) < 0 : return 2
 302                 else : return 0
 303
 304 def IsOnCircle (Edge, Center, Radius):
 305         """
 306         Function checks whether a given edge object belong to the periphery of a circle defined by its
 307         center and radius.
 308         Output can be 0 (does not belong) or 1 (belongs).
 309         If the reference Radius is null, the function returns 0
 310         Note that this function is basic in the sense that it only checks if the two border points of a
 311         given edge belong to the arc of reference.
 312         """
 313         if Radius == 0 : return 0
 314         else :
 315                 P1 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,0))
 316                 P2 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,1))
 317                 if abs(Distance2Pt(Center,P1)-Radius) < 1e-6 and abs(Distance2Pt(Center,P2)-Radius) < 1e-6:
 318                         return 1
 319                 else :
 320                         return 0
 321
 322 def CrossProd(V1,V2):
 323         """
 324         Determines the cross product of two 3D vectors
 325         """
 326         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]])
 327
 328 def QuarCylParam(PitchRatio):
 329         R = float(PitchRatio)/(PitchRatio+1)
 330         Eps = 1. - R
 331         X = (R+Eps/2.)*math.sin(math.pi/4)+Eps/2.
 332         N = int(math.floor((math.pi*R/4.)/(Eps/2.)))
 333         return N
 334
 335 def DotProd(V1,V2):
 336         """
 337         Determines the dot product of two 3D vectors
 338         """
 339         if len(V1)==2 : V1.append(0)
 340         if len(V2)==2 : V2.append(0)
 341
 342         return (V1[0]*V2[0]+V1[1]*V2[1]+V1[2]*V2[2])
 343
 344 def Distance2Pt(P1,P2):
 345         """
 346         Returns the distance between two points
 347         """
 348         return (math.sqrt((P1[0]-P2[0])**2+(P1[1]-P2[1])**2+(P1[2]-P2[2])**2))
 349
 350 def ApplyConstant1DMesh (ParentMsh, Edge, Nseg):
 351         Reg1D = ParentMsh.Segment(geom=Edge)
 352         Len = Reg1D.NumberOfSegments(Nseg)
 353
 354 def Apply1DProjMesh (ParentMsh, Edge, Ref):
 355         Proj1D = ParentMsh.Projection1D(geom=Edge)
 356         SrcEdge = Proj1D.SourceEdge(Ref,None,None,None)
 357
 358 def EdgeLength (Edge):
 359         """
 360         This function returns the edge object length.
 361         """
 362         P1 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,0))
 363         P2 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,1))
 364         return Distance2Pt(P1,P2)
 365
 366
 367 def D2R (Angle):
 368         return Angle*math.pi/180
 369
 370 def R2D (Angle):
 371         return Angle*180/math.pi
 372
 373 def F2D (FloatNumber):
 374         return round(FloatNumber*100.)/100.
 375
 376 def BezierGen (PointA, PointB, AngleA, AngleB):
 377
 378         if AngleA == 0 and AngleB == 0 : return (geompy.MakeEdge(PointA, PointB))
 379         else :
 380                 A = geompy.PointCoordinates(PointA)
 381                 B = geompy.PointCoordinates(PointB)
 382                 dAB = Distance2Pt(A,B)
 383                 dAC = dAB * (math.tan(AngleA)*math.tan(AngleB)) / (math.sin(AngleA) * ( math.tan(AngleA)+math.tan(AngleB) ) )
 384                 AngleOX_AB = math.acos((B[0]-A[0])/dAB)
 385                 PointC = geompy.MakeVertex(A[0]+math.cos(AngleA+AngleOX_AB)*dAC,A[1]+math.sin(AngleA+AngleOX_AB)*dAC,0)
 386                 CurveACB = geompy.MakeBezier([PointA,PointC,PointB])
 387                 return CurveACB
 388
 389 def GetSideAngleForBezier (PointA , PointB):
 390         """
 391         This function takes for input two points A and B where the bezier line is needed. It calculates the incident
 392         angle needed at point A so that the final curve is either at 0 or 90 degrees from the x'Ox axis
 393         """
 394         A = geompy.PointCoordinates(PointA)
 395         B = geompy.PointCoordinates(PointB)
 396         ABx = B[0]-A[0]
 397         dAB = Distance2Pt(A,B)
 398         Alpha = math.acos(ABx/dAB)
 399         #print "New angle request"
 400         #print ABx, dAB, R2D(Alpha)
 401         if Alpha < math.pi/4 :
 402                 #print "returning", R2D(-Alpha)
 403                 return -Alpha
 404         elif Alpha < 3*math.pi/4 :
 405                 #print "returning", R2D(-(Alpha-math.pi/2))
 406                 return -(Alpha-math.pi/2)
 407         else :
 408                 #print "returning", R2D(-(Alpha-math.pi))
 409                 return -(Alpha-math.pi)
 410
 411 def VecDivRatio (Vec1, Vec2):
 412         """
 413         This function tries to find the ratio of Vec1 on Vec2 while neglecting any zero term in Vec1. This is used afterwards
 414         for determining the global mesh parameter from automatically detected directional mesh params. If no compatibility is
 415         possible, the function returns -1
 416         """
 417         Vec3 = []
 418         for i in range(len(Vec1)) :
 419                 Vec3.append(float(Vec1[i])/Vec2[i])
 420         Ratio=[]
 421         for i in Vec3 :
 422                 if not (abs(i)<1e-7) : Ratio.append(i)
 423         if Ratio :
 424                 if min(Ratio) == max(Ratio) and min(Ratio)==int(min(Ratio)) : return(min(Ratio))
 425                 else : return -1
 426         else :
 427                 return -2
 428
 429
 430 def ReduceRatio (dx, dy):
 431         """
 432         This function transforms a decimal ratio into a scale between two integers, for example : [0.2,0.05] --> [4,1] ;
 433         """
 434         Output = [0,0]
 435         ratio = float(dy)/dx
 436         if isinteger(ratio) : return [1,ratio]
 437         elif dx == 1 :                  # when this function is called recursively!
 438                 for i in range(1,20) :  # searches over 20 decimals
 439                         if isinteger(ratio * (10**i) ) :
 440                                 Output = GetScale((10**i),int(round(ratio * (10**i) ) ) )
 441                                 break
 442                         else :
 443                                 for n in range(0,i) :
 444                                         if isinteger(ratio * ( 10**(i)-10**(n) )) :
 445                                                 Output = GetScale( 10**(i)-10**(n)  ,  int(round(ratio * ( 10**(i)-10**(n) ) ) ) )
 446                                                 break
 447                                 if not (Output==[0,0]) : break
 448                 return Output
 449         else :
 450                 for i in range(1,10) :  # searches over 10 decimals
 451                         if isinteger(ratio * (10**i) ) :
 452                                 Output = GetScale((10**i),int(round(ratio * (10**i) ) ) )
 453                                 break
 454                         else :
 455                                 for n in range(0,i) :
 456                                         if isinteger(ratio * ( 10**(i)-10**(n) )) :
 457                                                 Output = GetScale( 10**(i)-10**(n)  ,  int(round(ratio * ( 10**(i)-10**(n) ) ) ) )
 458                                                 break
 459                                 if not (Output==[0,0]) : break
 460
 461                 if Output == [0,0] :
 462                         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..."
 463                         if dy > dx :
 464                                 A = ReduceRatio (dx, dy-dx)
 465                                 return ([A[0],A[1]+A[0]])
 466                         else :
 467                                 A = ReduceRatio (dy, dx-dy)
 468                                 return ([A[1]+A[0],A[0]])
 469
 470                 else : return Output
 471
 472 def GetScale (X,Y):
 473         """
 474         This function is called within ReduceRatio and aims to reduce down two integers X and Y by dividing them with their common divisors;
 475         Example: 25 and 5 ---> 5 and 1 / 63 and 12 ---> 21 and 4
 476         """
 477         MaxDiv = max(X,Y)
 478         Divisor = 2             # Initializing the divisor
 479         while MaxDiv >= Divisor :
 480                 X0 = 0
 481                 Y0 = 0
 482                 if not(X%Divisor) :
 483                         X0 = X/Divisor
 484                         MaxDiv = max(MaxDiv,X0)
 485                 if not(Y%Divisor) :
 486                         Y0 = Y/Divisor
 487                         MaxDiv = max(MaxDiv,Y0)
 488                 if (X0*Y0) :
 489                         X = X0
 490                         Y = Y0
 491                 else :
 492                         Divisor = Divisor + 1
 493         return [X,Y]
 494
 495 def isinteger (x) :
 496         """
 497         This functions applies a simple check if the entered value is an integer
 498         """
 499         x = float('%.5f' % (x))         #Truncate x to 5 digits after the decimal point
 500         if math.ceil(x) == math.floor(x) : return True
 501         else : return False
 502 ##########################################################################################
 503 # Below this are the functions that create the elementary forms for the macro objects
 504 ##########################################################################################
 505
 506 def ElemBox11 ():
 507         """
 508         This function returns a simple square face of 1 side length
 509         """
 510         RectFace = geompy.MakeFaceHW(1, 1, 1)
 511         return RectFace
 512
 513 def ElemBox42 ():
 514         """
 515         This function returns a square face of 1 side length, partitioned
 516         according to the elementary 4 to 2 reductor method
 517         """
 518         OrigRectFace = geompy.MakeFaceHW(1, 1, 1)
 519
 520         SouthPt1 = geompy.MakeVertex (-.25, -.5, 0)
 521         SouthPt2 = geompy.MakeVertex (0, -.5, 0)
 522         SouthPt3 = geompy.MakeVertex (.25, -.5, 0)
 523         WestPt1 = geompy.MakeVertex (-.5, -.5+1./3, 0)
 524         WestPt2 = geompy.MakeVertex (-.5, -.5+2./3, 0)
 525         EastPt1 = geompy.MakeVertex (.5, -.5+1./3, 0)
 526         EastPt2 = geompy.MakeVertex (.5, -.5+2./3, 0)
 527         NorthPt = geompy.MakeVertex (0, .5, 0)
 528         MidPt1 = geompy.MakeVertex (0, .05, 0)
 529         MidPt2 = geompy.MakeVertex (.2, -.18, 0)
 530         MidPt3 = geompy.MakeVertex (0, -.28, 0)
 531         MidPt4 = geompy.MakeVertex (-.2, -.18, 0)
 532
 533         Cutter = []
 534         Cutter.append(geompy.MakeEdge(SouthPt2, MidPt3))
 535         Cutter.append(geompy.MakeEdge(MidPt1, NorthPt))
 536         Cutter.append(BezierGen(SouthPt1, MidPt4, GetSideAngleForBezier(SouthPt1,MidPt4), D2R(15)))
 537         Cutter.append(BezierGen(SouthPt3, MidPt2, GetSideAngleForBezier(SouthPt3,MidPt2), D2R(-15)))
 538         Cutter.append(BezierGen(WestPt1, MidPt4, GetSideAngleForBezier(WestPt1,MidPt4), D2R(-10)))
 539         Cutter.append(BezierGen(EastPt1, MidPt2, GetSideAngleForBezier(EastPt1,MidPt2), D2R(10)))
 540         Cutter.append(BezierGen(WestPt2, MidPt1, GetSideAngleForBezier(WestPt2,MidPt1), D2R(-10)))
 541         Cutter.append(BezierGen(EastPt2, MidPt1, GetSideAngleForBezier(EastPt2,MidPt1), D2R(10)))
 542         Cutter.append(BezierGen(MidPt2, MidPt1, D2R(-15), D2R(-15)))
 543         Cutter.append(BezierGen(MidPt3, MidPt2, D2R(10), D2R(15)))
 544         Cutter.append(BezierGen(MidPt3, MidPt4, D2R(-10), D2R(-15)))
 545         Cutter.append(BezierGen(MidPt4, MidPt1, D2R(15), D2R(15)))
 546
 547         RectFace = geompy.MakePartition([OrigRectFace],Cutter, [], [],4, 0, [], 0)      #Creating the partition object
 548         #i=1
 549         #for SingleCut in Cutter :
 550         #       geompy.addToStudy(SingleCut,'Cutter'+str(i))
 551         #       i = i+1
 552         #geompy.addToStudy(RectFace,'RectFace')
 553         return RectFace
 554
 555 def ElemEdge32 ():
 556         """
 557         This function returns a square face of 1 side length, partitioned
 558         according to the elementary edge with 3 to 2 reductor
 559         """
 560         OrigRectFace = geompy.MakeFaceHW(1., 1., 1)
 561
 562         SouthPt1 = geompy.MakeVertex (-1./6, -0.5, 0.)
 563         SouthPt2 = geompy.MakeVertex ( 1./6, -0.5, 0.)
 564         WestPt1 = geompy.MakeVertex  (-0.5, -1./6, 0.)
 565         WestPt2 = geompy.MakeVertex  (-0.5,  1./6, 0.)
 566         EastPt = geompy.MakeVertex   ( 0.5, 0., 0.)
 567         NorthPt = geompy.MakeVertex  (0., 0.5, 0.)
 568
 569         MidPt1 = geompy.MakeVertex (-0.2, -0.2, 0.)
 570         MidPt2 = geompy.MakeVertex ( -0.02,  -0.02, 0.)
 571
 572         Cutter = []
 573         Cutter.append(BezierGen(SouthPt1,  MidPt1,  GetSideAngleForBezier(SouthPt1,MidPt1) , D2R(-5)))
 574         Cutter.append(BezierGen( WestPt1,  MidPt1,  GetSideAngleForBezier(WestPt1 ,MidPt1) , D2R(-5)))
 575         Cutter.append(BezierGen(SouthPt2,  MidPt2,  GetSideAngleForBezier(SouthPt2,MidPt2) , D2R(-10)))
 576         Cutter.append(BezierGen(  EastPt,  MidPt2,  GetSideAngleForBezier(EastPt  ,MidPt2) , D2R(5)))
 577         Cutter.append(BezierGen( WestPt2,  MidPt2,  GetSideAngleForBezier(WestPt2 ,MidPt2) , D2R(-10)))
 578         Cutter.append(BezierGen(  MidPt2, NorthPt,  GetSideAngleForBezier(NorthPt ,MidPt2) , D2R(-5)))
 579
 580         Cutter.append(geompy.MakeEdge(MidPt1, MidPt2))
 581
 582         RectFace = geompy.MakePartition([OrigRectFace],Cutter, [], [],4, 0, [], 0)      #Creating the partition object
 583         #i=1
 584         #for SingleCut in Cutter :
 585         #       geompy.addToStudy(SingleCut,'Cutter'+str(i))
 586         #       i = i+1
 587         #geompy.addToStudy(RectFace,'RectFace')
 588         return RectFace
 589
 590 def Quadrangler (Points):
 591         """
 592         This function returns a quadranglar face based on four points, non of which 3 are non-colinear.
 593         The points are defined by their 2D [(x1,y1),(x2,y2)..] coordinates.
 594         Note that the list of points is already arranged upon the creation in MacObject
 595         """
 596         Pt = []
 597         for Point in Points: Pt.append(geompy.MakeVertex(Point[0], Point[1], 0))
 598         # The first point is added at the end of the list in order to facilitate the line creation
 599         Pt.append(Pt[0])
 600         #Draw the lines in order to form the 4 side polygon
 601         Ln=[]
 602         for i in range(4) : Ln.append(geompy.MakeLineTwoPnt(Pt[i],Pt[i+1]))
 603         RectFace = geompy.MakeQuad (Ln[0],Ln[1],Ln[2],Ln[3])
 604         return RectFace
 605
 606 def ElemQuartCyl(K):
 607         """
 608         This function returns a quarter cylinder to box relay of 1 side length, partitioned
 609         with a pitch ratio of K, In other words the side of the box is R*(1+(1/K))
 610         """
 611         R = 10.*float(K)/(K+1)
 612         Eps = 10.- R
 613
 614         Config.theStudy.SetReal("R"  , R)
 615         Config.theStudy.SetReal("minusR"  , -R)
 616         Config.theStudy.SetReal("Eps", Eps)
 617
 618         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])
 619         CylFace = geompy.MakeFace(CylWire, 1)
 620
 621         SouthPt = geompy.MakeVertex (R+Eps/2., 0., 0)
 622         SouthWestPt = geompy.MakeVertex ( 0.,0., 0)   #The origin can be used for practical partionning objectifs
 623         WestPt = geompy.MakeVertex  (0., R+Eps/2., 0)
 624
 625         N = int(math.floor((math.pi*R/4.)/(Eps/2.)))
 626         X = 10.*(1.-1./(N+1))
 627
 628
 629         EastPt = geompy.MakeVertex  (10.0,  X, 0.)
 630         NorthPt = geompy.MakeVertex   ( X, 10.0, 0.)
 631
 632         DivFactor = 8./(F2D(math.log(K))-0.223)
 633         #MidPt = geompy.MakeVertex ((R+Eps)*math.cos(math.pi/4), (R+Eps)*math.sin(math.pi/4), 0.)
 634         MidPt = geompy.MakeVertex (X-Eps/DivFactor, X-Eps/DivFactor, 0.)
 635
 636         Cutter = []
 637         Cutter.append(BezierGen(SouthWestPt,  MidPt,  GetSideAngleForBezier(SouthWestPt,MidPt) , D2R(-5)))
 638         Cutter.append(BezierGen( EastPt, MidPt,  GetSideAngleForBezier(EastPt,MidPt) , D2R(5)))
 639         Cutter.append(BezierGen( MidPt, NorthPt, (-1)**((K<1.25)*1)*D2R(-5), GetSideAngleForBezier(NorthPt,MidPt)))
 640         SMBezier = BezierGen( SouthPt,  MidPt,  GetSideAngleForBezier(SouthPt ,MidPt) , D2R((K<1.25)*180-5))
 641         WMBezier = BezierGen( WestPt,  MidPt,  GetSideAngleForBezier(WestPt, MidPt) , D2R(-5))
 642         Cutter.append(WMBezier)
 643         Cutter.append(SMBezier)
 644
 645         for i in range(1,N) :
 646                 # Determining intermediate points on the bezier lines and then performing additional cuts
 647
 648                 TempAnglePlus = (math.pi/4)*(1+float(i)/N)
 649                 SectionResult = CutnGroup.Go(WMBezier, [(0,0,0,math.sin(TempAnglePlus),-math.cos(TempAnglePlus),0)], [1], ['Dummy'], 0)
 650                 TempPt1 = SectionResult[1][0]
 651                 TempPt11 = geompy.MakeVertex  ((N-i)*X/N, 10., 0)
 652
 653                 TempAngleMinus = (math.pi/4)*(1-float(i)/N)
 654                 SectionResult = CutnGroup.Go(SMBezier, [(0,0,0,math.sin(TempAngleMinus),-math.cos(TempAngleMinus),0)], [1], ['Dummy'], 0)
 655                 TempPt2 = SectionResult[1][0]
 656                 TempPt21 = geompy.MakeVertex  (10.,  (N-i)*X/N, 0)
 657
 658                 Cutter.append(geompy.MakeEdge(SouthWestPt, TempPt1))
 659                 Cutter.append(geompy.MakeEdge(SouthWestPt, TempPt2))
 660                 Cutter.append(geompy.MakeEdge(TempPt1, TempPt11))
 661                 Cutter.append(geompy.MakeEdge(TempPt2, TempPt21))
 662
 663         CylFace = geompy.MakePartition([CylFace],Cutter, [], [],4, 0, [], 0)    #Creating the partition object
 664         CylFace = geompy.MakeTranslation(CylFace, -5., -5., 0.0)
 665
 666         return CylFace
 667
 668 def CompatibilityTest(MacObject):
 669         Type = MacObject.Type
 670         if Type == 'Box11' :
 671                 BaseDirPar = [1,1,1,1]
 672                 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
 673         elif Type == 'Box42' :
 674                 BaseDirPar = {'SN' : lambda : [3, 3, 4, 2],
 675                               'NS' : lambda : [3, 3, 2, 4],
 676                               'EW' : lambda : [2, 4, 3, 3],
 677                               'WE' : lambda : [4, 2, 3, 3], }[MacObject.MeshPar[1]]()
 678                 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
 679         elif Type == 'BoxAng32' :
 680                 BaseDirPar = {'NE' : lambda : [3, 2, 3, 2],
 681                               'NW' : lambda : [2, 3, 3, 2],
 682                               'SW' : lambda : [2, 3, 2, 3],
 683                               'SE' : lambda : [3, 2, 2, 3], }[MacObject.MeshPar[1]]()
 684                 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
 685         elif Type == 'CompBox' :
 686                 #print "dx is: ", MacObject.GeoPar[1][1], ". dy is: ",MacObject.GeoPar[1][0]
 687                 ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0], MacObject.GeoPar[1][1])
 688                 #print ReducedRatio
 689                 BaseDirPar = [ReducedRatio[1], ReducedRatio[1], ReducedRatio[0], ReducedRatio[0]]
 690                 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
 691
 692         elif Type == 'QuartCyl' :
 693                 N = QuarCylParam(MacObject.MeshPar[2])+1
 694                 BaseDirPar = {'NE' : lambda : [2, N, 2, N],
 695                               'NW' : lambda : [N, 2, 2, N],
 696                               'SW' : lambda : [N, 2, N, 2],
 697                               'SE' : lambda : [2, N, N, 2], }[MacObject.MeshPar[1]]()
 698                 return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
 699         elif Type == 'CompBoxF' :
 700                 RealRatio = MacObject.GeoPar[1][1]/MacObject.GeoPar[1][0]
 701                 Xd = 0
 702                 Yd = 0
 703                 if MacObject.DirectionalMeshParams[2]+MacObject.DirectionalMeshParams[3] :
 704                         A = int(max(MacObject.DirectionalMeshParams[2:4]))
 705                         Xd = int(VecDivRatio([A,0,0,0], [1,1,1,1]))
 706                 if MacObject.DirectionalMeshParams[0]+MacObject.DirectionalMeshParams[1] :
 707                         A = int(max(MacObject.DirectionalMeshParams[0:2]))
 708                         Yd = int(VecDivRatio([0,0,A,0], [1,1,1,1]))
 709
 710                 if Xd == 0 and Yd : Xd = int(round(Yd/RealRatio))
 711                 elif Yd == 0 : Yd = int(round(RealRatio*Xd))
 712
 713                 return [Xd,Yd]
 714         elif Type == 'NonOrtho' :
 715                 MeanDX = 0.5*(IntLen(MacObject.DirBoundaries(0))+IntLen(MacObject.DirBoundaries(1)))
 716                 MeanDY = 0.5*(IntLen(MacObject.DirBoundaries(2))+IntLen(MacObject.DirBoundaries(3)))
 717                 RealRatio = MeanDY/MeanDX
 718                 Xd = 0
 719                 Yd = 0
 720                 if MacObject.DirectionalMeshParams[2]+MacObject.DirectionalMeshParams[3] :
 721                         A = int(max(MacObject.DirectionalMeshParams[2:4]))
 722                         Xd = int(VecDivRatio([A,0,0,0], [1,1,1,1]))
 723                 if MacObject.DirectionalMeshParams[0]+MacObject.DirectionalMeshParams[1] :
 724                         A = int(max(MacObject.DirectionalMeshParams[0:2]))
 725                         Yd = int(VecDivRatio([0,0,A,0], [1,1,1,1]))
 726
 727                 if Xd == 0 and Yd : Xd = int(round(Yd/RealRatio))
 728                 elif Yd == 0 : Yd = int(round(RealRatio*Xd))
 729
 730                 return [Xd,Yd]
 731
 732 def IntLen (Interval) :
 733         """
 734         This function returns the length of a given interval even if the latter is not sorted correctly.
 735         """
 736         return abs(Interval[1]-Interval[0])
 737
 738 def NextTo (RefBox, Direction, Extension):
 739         """
 740         This functions returns geometrical parameters for easy positioning of neighbouring objects.
 741         The input (RefBox) and output are in the form :  [(X0,Y0),(DX,DY)]
 742         """
 743         X0_0 = RefBox[0][0]
 744         Y0_0 = RefBox[0][1]
 745         DX_0 = RefBox[1][0]
 746         DY_0 = RefBox[1][1]
 747
 748         DirectionalCoef = {'Above' : lambda : [ 0, 1],
 749                            'Below' : lambda : [ 0,-1],
 750                            'Right' : lambda : [ 1, 0],
 751                            'Left ' : lambda : [-1, 0], }[Direction]()
 752
 753         X0_1 = X0_0+ DirectionalCoef[0] * (DX_0/2.+Extension/2.)
 754         DX_1 = abs(DirectionalCoef[0]) * (Extension) + abs(DirectionalCoef[1])*DX_0
 755         Y0_1 = Y0_0+ DirectionalCoef[1] * (DY_0/2.+Extension/2.)
 756         DY_1 = abs(DirectionalCoef[1]) * (Extension) + abs(DirectionalCoef[0])*DY_0
 757
 758         return [(X0_1,Y0_1),(DX_1,DY_1)]
 759
 760 def GeomMinMax (PtA, PtB):
 761         """
 762         This function returns geometrical parameters in the format  [(X0,Y0),(DX,DY)]. The input being
 763         the coordinates of two points (Xa,Ya), (Xb,Yb).
 764         """
 765         # First test that the vector relying the two points is oblique
 766         AB = [PtB[0]- PtA[0],PtB[1]- PtA[1]]
 767         if 0 in AB :
 768                 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")
 769                 return -1
 770         else:
 771                 X0 = 0.5*(PtA[0]+PtB[0])
 772                 Y0 = 0.5*(PtA[1]+PtB[1])
 773                 DX = abs(AB[0])
 774                 DY = abs(AB[1])
 775                 return [(X0,Y0),(DX,DY)]
 776
 777 def AddIfDifferent (List, Element):
 778         if not(Element in List):
 779                 List = List+(Element,)
 780         return List
 781
 782 def IndexMultiOcc (Array,Element) :
 783         """
 784         This functions returns the occurrences indices of Element in Array.
 785         As opposed to Array.index(Element) method, this allows determining
 786         multiple entries rather than just the first one!
 787         """
 788         Output = []
 789         try : Array.index(Element)
 790         except ValueError : print "No more occurrences"
 791         else : Output.append(Array.index(Element))
 792
 793         if not(Output == []) and len(Array) > 1 :
 794                 for index, ArrElem in enumerate(Array[Output[0]+1:]) :
 795                         if ArrElem == Element : Output.append(index+Output[0]+1)
 796
 797         return Output
 798
 799 def SortList (ValList, CritList):
 800         Output = []
 801         SortedCritList = copy.copy(CritList)
 802         SortedCritList.sort()
 803         for i in range(0,len(ValList)):
 804                 if i > 0 :
 805                         if not(SortedCritList[i]==SortedCritList[i-1]):
 806                                 index = IndexMultiOcc(CritList,SortedCritList[i])
 807                                 Output= Output + [ValList[j] for j in index]
 808                 else :
 809                         index = IndexMultiOcc(CritList,SortedCritList[i])
 810                         Output= Output + [ValList[j] for j in index]
 811
 812         return Output
 813
 814 def SortPoints(Points):
 815         """
 816         This function sorts a list of the coordinates of N points as to start at
 817         an origin that represents Xmin and Xmax and then proceed in a counter
 818         clock-wise sense
 819         """
 820         NbPts = len(Points)
 821         Xmin = min([Points[i][0] for i in range(NbPts)])
 822         Ymin = min([Points[i][1] for i in range(NbPts)])
 823         Xmax = max([Points[i][0] for i in range(NbPts)])
 824         Ymax = max([Points[i][1] for i in range(NbPts)])
 825         Crit = [(abs(Point[0]-Xmin)+0.1*(Xmax-Xmin))*(abs(Point[1]-Ymin)+0.1*(Ymax-Ymin)) for Point in Points]
 826         #print "Input Points      : ", Points
 827         #print "Sorting Criterion : ", Crit
 828         Order = SortList (range(NbPts), Crit)
 829         #print "Sorted Results    : ", Order
 830         Output = []
 831         Output.append(Points[Order[0]])
 832
 833         Point0 = Points[Order[0]]
 834         #print "Reference point :", Point0
 835
 836         V = [[Point1[0]-Point0[0],Point1[1]-Point0[1]] for Point1 in Points]
 837         Cosines = [-(vec[0]-1E-10)/(math.sqrt(DotProd(vec,vec)+1e-25)) for vec in V]
 838         #print "Cosines criterion :", Cosines
 839         Order = SortList(range(NbPts),Cosines)
 840         #print "Ordered points:", Order
 841         for PtIndex in Order[:-1]: Output.append(Points[PtIndex])
 842
 843         return Output
 844