src/Tools/MacMesh/MacMesh/GenFunctions.py

   1 # Copyright (C) 2014-2015  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( Config.theStudy )
  31
  32 from salome.smesh import smeshBuilder
  33 smesh = smeshBuilder.New( Config.theStudy )
  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 = copy.copy(CritList)
 829         SortedCritList.sort()
 830         for i in range(0,len(ValList)):
 831                 if i > 0 :
 832                         if not(SortedCritList[i]==SortedCritList[i-1]):
 833                                 index = IndexMultiOcc(CritList,SortedCritList[i])
 834                                 Output= Output + [ValList[j] for j in index]
 835                 else :
 836                         index = IndexMultiOcc(CritList,SortedCritList[i])
 837                         Output= Output + [ValList[j] for j in index]
 838
 839         return Output
 840
 841 def SortPoints(Points):
 842         """
 843         This function sorts a list of the coordinates of N points as to start at
 844         an origin that represents Xmin and Xmax and then proceed in a counter
 845         clock-wise sense
 846         """
 847         NbPts = len(Points)
 848         Xmin = min([Points[i][0] for i in range(NbPts)])
 849         Ymin = min([Points[i][1] for i in range(NbPts)])
 850         Xmax = max([Points[i][0] for i in range(NbPts)])
 851         Ymax = max([Points[i][1] for i in range(NbPts)])
 852         Crit = [(abs(Point[0]-Xmin)+0.1*(Xmax-Xmin))*(abs(Point[1]-Ymin)+0.1*(Ymax-Ymin)) for Point in Points]
 853         #print "Input Points      : ", Points
 854         #print "Sorting Criterion : ", Crit
 855         Order = SortList (range(NbPts), Crit)
 856         #print "Sorted Results    : ", Order
 857         Output = []
 858         Output.append(Points[Order[0]])
 859
 860         Point0 = Points[Order[0]]
 861         #print "Reference point :", Point0
 862
 863         V = [[Point1[0]-Point0[0],Point1[1]-Point0[1]] for Point1 in Points]
 864         Cosines = [-(vec[0]-1E-10)/(math.sqrt(DotProd(vec,vec)+1e-25)) for vec in V]
 865         #print "Cosines criterion :", Cosines
 866         Order = SortList(range(NbPts),Cosines)
 867         #print "Ordered points:", Order
 868         for PtIndex in Order[:-1]: Output.append(Points[PtIndex])
 869
 870         return Output
 871