src/Tools/MacMesh/MacMesh/GenFunctions.py

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