src/Tools/MacMesh/MacMesh/GenFunctions.py

   1 # Copyright (C) 2014-2016  EDF R&D
   2 #
   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.
   7 #
   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.
  12 #
  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
  16 #
  17 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 #
  19
  20
  21
  22 # In this file are all the generation functions for manipulating the different created macro-objects
  23
  24 import math, copy
  25 import Config
  26 import CutnGroup
  27 import CompositeBox
  28
  29 from salome.geom import geomBuilder
  30 geompy = geomBuilder.New()
  31
  32 from salome.smesh import smeshBuilder
  33 smesh = smeshBuilder.New()
  34
  35 ##########################################################################################################
  36
  37 def Box11 (MacObject):
  38     if Config.debug : print("Generating regular box")
  39
  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)
  42
  43     MacObject.GeoChildren.append(RectFace)
  44     MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
  45
  46     if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
  47
  48     if Config.publish :
  49         MacObject.Mesh.append(smesh.Mesh(RectFace))                     # Creation of a new mesh
  50         Quad2D = MacObject.Mesh[0].Quadrangle()                 # Applying a quadrangle hypothesis
  51
  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])
  55
  56         MacObject.Mesh[0].Compute()                     # Generates the mesh
  57
  58     MacObject.DirectionalMeshParams = [MacObject.MeshPar[0],MacObject.MeshPar[0],MacObject.MeshPar[0],MacObject.MeshPar[0]]
  59
  60     MacObject.status = 1
  61     Config.ListObj.append(MacObject)
  62     return MacObject
  63
  64 ##########################################################################################################
  65
  66 def Box42 (MacObject):
  67     if Config.debug : print("Generating box 4-2 reducer")
  68
  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)
  77
  78     MacObject.GeoChildren.append(RectFace)
  79     MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
  80
  81     if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
  82
  83     if Config.publish :
  84         MacObject.Mesh.append(smesh.Mesh(RectFace))                     # Creation of a new mesh
  85         Quad2D = MacObject.Mesh[0].Quadrangle()                 # Applying a quadrangle hypothesis
  86
  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])
  90
  91         MacObject.Mesh[0].Compute()                     # Generates the mesh
  92
  93     MacObject.status = 1
  94
  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]]()
 100
 101     Config.ListObj.append(MacObject)
 102     return MacObject
 103
 104
 105 ##########################################################################################################
 106
 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)
 117
 118     MacObject.GeoChildren.append(RectFace)
 119     MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
 120
 121     if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
 122
 123     if Config.publish :
 124         MacObject.Mesh.append(smesh.Mesh(RectFace))                     # Creation of a new mesh
 125         Quad2D = MacObject.Mesh[0].Quadrangle()                 # Applying a quadrangle hypothesis
 126
 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])
 130
 131         MacObject.Mesh[0].Compute()                     # Generates the mesh
 132
 133     MacObject.status = 1
 134
 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]]()
 140
 141     Config.ListObj.append(MacObject)
 142     return MacObject
 143 ##########################################################################################################
 144 def CompBox (MacObject):
 145     if Config.debug : print("Generating composite box")
 146
 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)
 149
 150     MacObject.GeoChildren.append(RectFace)
 151     MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
 152
 153     if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
 154
 155     if Config.publish :
 156         MacObject.Mesh.append(smesh.Mesh(RectFace))                     # Creation of a new mesh
 157         Quad2D = MacObject.Mesh[0].Quadrangle()                 # Applying a quadrangle hypothesis
 158
 159         EdgeIDs = geompy.SubShapeAllSorted(RectFace,6)      # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
 160
 161         ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1])
 162
 163         Reference = [0,0,0]
 164         Vec = [(1,0,0),(0,1,0)]
 165         for Edge in EdgeIDs:
 166             for i in range(0,2):
 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
 169                         Reference[i] = Edge
 170                         ApplyConstant1DMesh(MacObject.Mesh[0],Edge,int(round(ReducedRatio[i]*MacObject.MeshPar[0])))
 171                         break
 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])
 174                         break
 175
 176         MacObject.Mesh[0].Compute()                     # Generates the mesh
 177
 178     MacObject.DirectionalMeshParams = [MacObject.MeshPar[0]*ReducedRatio[1],MacObject.MeshPar[0]*ReducedRatio[1],MacObject.MeshPar[0]*ReducedRatio[0],MacObject.MeshPar[0]*ReducedRatio[0]]
 179
 180     MacObject.status = 1
 181     Config.ListObj.append(MacObject)
 182     return MacObject
 183
 184 ##########################################################################################################
 185
 186 def CompBoxF (MacObject):
 187     if Config.debug : print("Generating composite box")
 188
 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)
 191
 192     MacObject.GeoChildren.append(RectFace)
 193     MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
 194
 195     if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
 196
 197     if Config.publish :
 198         MacObject.Mesh.append(smesh.Mesh(RectFace))                     # Creation of a new mesh
 199         Quad2D = MacObject.Mesh[0].Quadrangle()                 # Applying a quadrangle hypothesis
 200
 201         EdgeIDs = geompy.SubShapeAllSorted(RectFace,6)      # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
 202
 203         #ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1])
 204
 205         Reference = [0,0,0]
 206         Vec = [(1,0,0),(0,1,0)]
 207         for Edge in EdgeIDs:
 208             for i in range(0,2):
 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
 211                         Reference[i] = Edge
 212                         ApplyConstant1DMesh(MacObject.Mesh[0],Edge,int(round(MacObject.MeshPar[0][i])))
 213                         break
 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])
 216                         break
 217
 218         MacObject.Mesh[0].Compute()                 # Generates the mesh
 219
 220     MacObject.DirectionalMeshParams = [MacObject.MeshPar[0][1],MacObject.MeshPar[0][1],MacObject.MeshPar[0][0],MacObject.MeshPar[0][0]]
 221
 222     MacObject.status = 1
 223     Config.ListObj.append(MacObject)
 224     return MacObject
 225 ##########################################################################################################
 226
 227 def NonOrtho (MacObject):
 228     if Config.debug : print("Generating Non-orthogonal quadrangle")
 229
 230     RectFace = Quadrangler (MacObject.PtCoor)
 231
 232     MacObject.GeoChildren.append(RectFace)
 233     MacObject.GeoChildrenNames.append("Quad_"+ str(len(Config.ListObj)+1))
 234
 235
 236     if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
 237
 238     if Config.publish :
 239         MacObject.Mesh.append(smesh.Mesh(RectFace))                     # Creation of a new mesh
 240         Quad2D = MacObject.Mesh[0].Quadrangle()                 # Applying a quadrangle hypothesis
 241
 242         EdgeIDs = geompy.SubShapeAllSorted(RectFace,6)      # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
 243
 244         #ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1])
 245
 246         Vec = [MacObject.DirVectors(i) for i in range(4)]
 247         for Edge in EdgeIDs:
 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])))
 251
 252         MacObject.Mesh[0].Compute()                 # Generates the mesh
 253
 254     MacObject.DirectionalMeshParams = [MacObject.MeshPar[0][1],MacObject.MeshPar[0][1],MacObject.MeshPar[0][0],MacObject.MeshPar[0][0]]
 255
 256     MacObject.status = 1
 257     Config.ListObj.append(MacObject)
 258     return MacObject
 259
 260 ##########################################################################################################
 261
 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)
 272
 273     MacObject.GeoChildren.append(RectFace)
 274     MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1))
 275
 276     if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames)
 277
 278     if Config.publish :
 279         MacObject.Mesh.append(smesh.Mesh(RectFace))                     # Creation of a new mesh
 280         Quad2D = MacObject.Mesh[0].Quadrangle()                 # Applying a quadrangle hypothesis
 281
 282         EdgeIDs = geompy.SubShapeAllSorted(RectFace,6)      # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary
 283         Reg1D = MacObject.Mesh[0].Segment()
 284
 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
 288
 289         Reg1D.NumberOfSegments(MacObject.MeshPar[0])
 290
 291         MacObject.Mesh[0].Compute()                     # Generates the mesh
 292
 293     MacObject.status = 1
 294
 295     x = MacObject.MeshPar[0]
 296     N = QuarCylParam(MacObject.MeshPar[2])+1
 297
 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]]()
 302
 303     Config.ListObj.append(MacObject)
 304     return MacObject
 305
 306 ##########################################################################################################
 307 # Below this are the elementary calculation/visualization functions
 308 ##########################################################################################################
 309
 310 def Publish (ObjToPublish,NamesToPublish):
 311     i = 0
 312     for GeoObj in ObjToPublish :
 313         geompy.addToStudy(GeoObj,NamesToPublish[i])
 314         i = i+1
 315
 316 def IsParallel (Edge, Vector):
 317     """
 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
 321     """
 322     if Vector == (0,0,0) : return 0
 323     else :
 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
 329         else : return 0
 330
 331 def IsOnCircle (Edge, Center, Radius):
 332     """
 333     Function checks whether a given edge object belong to the periphery of a circle defined by its
 334     center and radius.
 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.
 339     """
 340     if Radius == 0 : return 0
 341     else :
 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:
 345             return 1
 346         else :
 347             return 0
 348
 349 def CrossProd(V1,V2):
 350     """
 351     Determines the cross product of two 3D vectors
 352     """
 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]])
 354
 355 def QuarCylParam(PitchRatio):
 356     R = float(PitchRatio)/(PitchRatio+1)
 357     Eps = 1. - R
 358     X = (R+Eps/2.)*math.sin(math.pi/4)+Eps/2.
 359     N = int(math.floor((math.pi*R/4.)/(Eps/2.)))
 360     return N
 361
 362 def DotProd(V1,V2):
 363     """
 364     Determines the dot product of two 3D vectors
 365     """
 366     if len(V1)==2 : V1.append(0)
 367     if len(V2)==2 : V2.append(0)
 368
 369     return (V1[0]*V2[0]+V1[1]*V2[1]+V1[2]*V2[2])
 370
 371 def Distance2Pt(P1,P2):
 372     """
 373     Returns the distance between two points
 374     """
 375     return (math.sqrt((P1[0]-P2[0])**2+(P1[1]-P2[1])**2+(P1[2]-P2[2])**2))
 376
 377 def ApplyConstant1DMesh (ParentMsh, Edge, Nseg):
 378     Reg1D = ParentMsh.Segment(geom=Edge)
 379     Len = Reg1D.NumberOfSegments(Nseg)
 380
 381 def Apply1DProjMesh (ParentMsh, Edge, Ref):
 382     Proj1D = ParentMsh.Projection1D(geom=Edge)
 383     SrcEdge = Proj1D.SourceEdge(Ref,None,None,None)
 384
 385 def EdgeLength (Edge):
 386     """
 387     This function returns the edge object length.
 388     """
 389     P1 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,0))
 390     P2 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,1))
 391     return Distance2Pt(P1,P2)
 392
 393
 394 def D2R (Angle):
 395     return Angle*math.pi/180
 396
 397 def R2D (Angle):
 398     return Angle*180/math.pi
 399
 400 def F2D (FloatNumber):
 401     return round(FloatNumber*100.)/100.
 402
 403 def BezierGen (PointA, PointB, AngleA, AngleB):
 404
 405     if AngleA == 0 and AngleB == 0 : return (geompy.MakeEdge(PointA, PointB))
 406     else :
 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])
 414         return CurveACB
 415
 416 def GetSideAngleForBezier (PointA , PointB):
 417     """
 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
 420     """
 421     A = geompy.PointCoordinates(PointA)
 422     B = geompy.PointCoordinates(PointB)
 423     ABx = B[0]-A[0]
 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)
 430         return -Alpha
 431     elif Alpha < 3*math.pi/4 :
 432         #print "returning", R2D(-(Alpha-math.pi/2))
 433         return -(Alpha-math.pi/2)
 434     else :
 435         #print "returning", R2D(-(Alpha-math.pi))
 436         return -(Alpha-math.pi)
 437
 438 def VecDivRatio (Vec1, Vec2):
 439     """
 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
 443     """
 444     Vec3 = []
 445     for i in range(len(Vec1)) :
 446         Vec3.append(float(Vec1[i])/Vec2[i])
 447     Ratio=[]
 448     for i in Vec3 :
 449         if not (abs(i)<1e-7) : Ratio.append(i)
 450     if Ratio :
 451         if min(Ratio) == max(Ratio) and min(Ratio)==int(min(Ratio)) : return(min(Ratio))
 452         else : return -1
 453     else :
 454         return -2
 455
 456
 457 def ReduceRatio (dx, dy):
 458     """
 459     This function transforms a decimal ratio into a scale between two integers, for example : [0.2,0.05] --> [4,1] ;
 460     """
 461     Output = [0,0]
 462     ratio = float(dy)/dx
 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) ) ) )
 468                 break
 469             else :
 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) ) ) ) )
 473                         break
 474                 if not (Output==[0,0]) : break
 475         return Output
 476     else :
 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) ) ) )
 480                 break
 481             else :
 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) ) ) ) )
 485                         break
 486                 if not (Output==[0,0]) : break
 487
 488         if Output == [0,0] :
 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...")
 490             if dy > dx :
 491                 A = ReduceRatio (dx, dy-dx)
 492                 return ([A[0],A[1]+A[0]])
 493             else :
 494                 A = ReduceRatio (dy, dx-dy)
 495                 return ([A[1]+A[0],A[0]])
 496
 497         else : return Output
 498
 499 def GetScale (X,Y):
 500     """
 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
 503     """
 504     MaxDiv = max(X,Y)
 505     Divisor = 2             # Initializing the divisor
 506     while MaxDiv >= Divisor :
 507         X0 = 0
 508         Y0 = 0
 509         if not(X%Divisor) :
 510             X0 = X/Divisor
 511             MaxDiv = max(MaxDiv,X0)
 512         if not(Y%Divisor) :
 513             Y0 = Y/Divisor
 514             MaxDiv = max(MaxDiv,Y0)
 515         if (X0*Y0) :
 516             X = X0
 517             Y = Y0
 518         else :
 519             Divisor = Divisor + 1
 520     return [X,Y]
 521
 522 def isinteger (x) :
 523     """
 524     This functions applies a simple check if the entered value is an integer
 525     """
 526     x = float('%.5f' % (x))         #Truncate x to 5 digits after the decimal point
 527     if math.ceil(x) == math.floor(x) : return True
 528     else : return False
 529 ##########################################################################################
 530 # Below this are the functions that create the elementary forms for the macro objects
 531 ##########################################################################################
 532
 533 def ElemBox11 ():
 534     """
 535     This function returns a simple square face of 1 side length
 536     """
 537     RectFace = geompy.MakeFaceHW(1, 1, 1)
 538     return RectFace
 539
 540 def ElemBox42 ():
 541     """
 542     This function returns a square face of 1 side length, partitioned
 543     according to the elementary 4 to 2 reductor method
 544     """
 545     OrigRectFace = geompy.MakeFaceHW(1, 1, 1)
 546
 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)
 559
 560     Cutter = []
 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)))
 573
 574     RectFace = geompy.MakePartition([OrigRectFace],Cutter, [], [],4, 0, [], 0)      #Creating the partition object
 575     #i=1
 576     #for SingleCut in Cutter :
 577     #       geompy.addToStudy(SingleCut,'Cutter'+str(i))
 578     #       i = i+1
 579     #geompy.addToStudy(RectFace,'RectFace')
 580     return RectFace
 581
 582 def ElemEdge32 ():
 583     """
 584     This function returns a square face of 1 side length, partitioned
 585     according to the elementary edge with 3 to 2 reductor
 586     """
 587     OrigRectFace = geompy.MakeFaceHW(1., 1., 1)
 588
 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.)
 595
 596     MidPt1 = geompy.MakeVertex (-0.2, -0.2, 0.)
 597     MidPt2 = geompy.MakeVertex ( -0.02,  -0.02, 0.)
 598
 599     Cutter = []
 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)))
 606
 607     Cutter.append(geompy.MakeEdge(MidPt1, MidPt2))
 608
 609     RectFace = geompy.MakePartition([OrigRectFace],Cutter, [], [],4, 0, [], 0)      #Creating the partition object
 610     #i=1
 611     #for SingleCut in Cutter :
 612     #       geompy.addToStudy(SingleCut,'Cutter'+str(i))
 613     #       i = i+1
 614     #geompy.addToStudy(RectFace,'RectFace')
 615     return RectFace
 616
 617 def Quadrangler (Points):
 618     """
 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
 622     """
 623     Pt = []
 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
 626     Pt.append(Pt[0])
 627     #Draw the lines in order to form the 4 side polygon
 628     Ln=[]
 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])
 631     return RectFace
 632
 633 def ElemQuartCyl(K):
 634     """
 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))
 637     """
 638     R = 10.*float(K)/(K+1)
 639     Eps = 10.- R
 640
 641     Config.theStudy.SetReal("R"  , R)
 642     Config.theStudy.SetReal("minusR"  , -R)
 643     Config.theStudy.SetReal("Eps", Eps)
 644
 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)
 647
 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)
 651
 652     N = int(math.floor((math.pi*R/4.)/(Eps/2.)))
 653     X = 10.*(1.-1./(N+1))
 654
 655
 656     EastPt = geompy.MakeVertex  (10.0,  X, 0.)
 657     NorthPt = geompy.MakeVertex   ( X, 10.0, 0.)
 658
 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.)
 662
 663     Cutter = []
 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)
 671
 672     for i in range(1,N) :
 673         # Determining intermediate points on the bezier lines and then performing additional cuts
 674
 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)
 679
 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)
 684
 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))
 689
 690     CylFace = geompy.MakePartition([CylFace],Cutter, [], [],4, 0, [], 0)    #Creating the partition object
 691     CylFace = geompy.MakeTranslation(CylFace, -5., -5., 0.0)
 692
 693     return CylFace
 694
 695 def CompatibilityTest(MacObject):
 696     Type = MacObject.Type
 697     if Type == 'Box11' :
 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])
 715         #print ReducedRatio
 716         BaseDirPar = [ReducedRatio[1], ReducedRatio[1], ReducedRatio[0], ReducedRatio[0]]
 717         return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar))
 718
 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]
 728         Xd = 0
 729         Yd = 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]))
 736
 737         if Xd == 0 and Yd : Xd = int(round(Yd/RealRatio))
 738         elif Yd == 0 : Yd = int(round(RealRatio*Xd))
 739
 740         return [Xd,Yd]
 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
 745         Xd = 0
 746         Yd = 0
 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]))
 753
 754         if Xd == 0 and Yd : Xd = int(round(Yd/RealRatio))
 755         elif Yd == 0 : Yd = int(round(RealRatio*Xd))
 756
 757         return [Xd,Yd]
 758
 759 def IntLen (Interval) :
 760     """
 761     This function returns the length of a given interval even if the latter is not sorted correctly.
 762     """
 763     return abs(Interval[1]-Interval[0])
 764
 765 def NextTo (RefBox, Direction, Extension):
 766     """
 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)]
 769     """
 770     X0_0 = RefBox[0][0]
 771     Y0_0 = RefBox[0][1]
 772     DX_0 = RefBox[1][0]
 773     DY_0 = RefBox[1][1]
 774
 775     DirectionalCoef = {'Above' : lambda : [ 0, 1],
 776                        'Below' : lambda : [ 0,-1],
 777                        'Right' : lambda : [ 1, 0],
 778                        'Left ' : lambda : [-1, 0], }[Direction]()
 779
 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
 784
 785     return [(X0_1,Y0_1),(DX_1,DY_1)]
 786
 787 def GeomMinMax (PtA, PtB):
 788     """
 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).
 791     """
 792     # First test that the vector relying the two points is oblique
 793     AB = [PtB[0]- PtA[0],PtB[1]- PtA[1]]
 794     if 0 in AB :
 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")
 796         return -1
 797     else:
 798         X0 = 0.5*(PtA[0]+PtB[0])
 799         Y0 = 0.5*(PtA[1]+PtB[1])
 800         DX = abs(AB[0])
 801         DY = abs(AB[1])
 802         return [(X0,Y0),(DX,DY)]
 803
 804 def AddIfDifferent (List, Element):
 805     if not(Element in List):
 806         List = List+(Element,)
 807     return List
 808
 809 def IndexMultiOcc (Array,Element) :
 810     """
 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!
 814     """
 815     Output = []
 816     try : Array.index(Element)
 817     except ValueError : print("No more occurrences")
 818     else : Output.append(Array.index(Element))
 819
 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)
 823
 824     return Output
 825
 826 def SortList (ValList, CritList):
 827     Output = []
 828     SortedCritList = sorted(copy.copy(CritList))
 829     for i in range(0,len(ValList)):
 830         if i > 0 :
 831             if not(SortedCritList[i]==SortedCritList[i-1]):
 832                 index = IndexMultiOcc(CritList,SortedCritList[i])
 833                 Output= Output + [ValList[j] for j in index]
 834         else :
 835             index = IndexMultiOcc(CritList,SortedCritList[i])
 836             Output= Output + [ValList[j] for j in index]
 837
 838     return Output
 839
 840 def SortPoints(Points):
 841     """
 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
 844     clock-wise sense
 845     """
 846     NbPts = len(Points)
 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
 856     Output = []
 857     Output.append(Points[Order[0]])
 858
 859     Point0 = Points[Order[0]]
 860     #print "Reference point :", Point0
 861
 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])
 868
 869     return Output