src/Tools/MacMesh/MacMesh/CompositeBoxF.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
  25 # INTRODUCTION HERE
  26
  27 import sys, math, copy, commands
  28 CWD = commands.getoutput('pwd')
  29 sys.path.append(CWD)
  30
  31 from MacObject import *
  32 import Config, GenFunctions
  33
  34 def CompositeBoxF (Pt1 , Pt2 , Pt3 , Pt4 , **args ) :
  35         [Pt1 , Pt2 , Pt3 , Pt4] = GenFunctions.SortPoints([Pt1 , Pt2 , Pt3 , Pt4])
  36         if args.__contains__('groups') :
  37                 GroupNames = args['groups']
  38         else : GroupNames = [None, None, None, None]
  39         # Create a full NonOrtho box just to inherit, globally, the mesh parameters of bounding objects
  40         dummy = MacObject('NonOrtho',[Pt1,Pt2,Pt3,Pt4],['auto'],publish=0)
  41         # Save the existing number of segments on each direction
  42         ExistingSeg0 = Config.ListObj[-1].DirectionalMeshParams
  43         Convention = [2,3,0,1]
  44         ExistingSegments = [ExistingSeg0[Convention[i]] for i in range(4)]
  45         # Save Boundary lengths on each direction
  46         BoundaryLengths = [IntLen(dummy.DirBoundaries(i)) for i in range(4) ]
  47         # Calculate global mesh element size on each direction
  48         GlobalDelta = [1.*BoundaryLengths[i]/ExistingSegments[i] for i in range(4) ]
  49         print "GlobalDelta :",GlobalDelta
  50         # Sort the connection list for the full Box
  51         [(X0,Y0),(DX,DY)] = dummy.GeoPar
  52         ObjIDLists = SortObjLists(Config.Connections[-1],X0 , Y0 , DX , DY )
  53         [Xmin,Xmax,Ymin,Ymax] = dummy.Boundaries() # Used for groups determination
  54         RemoveLastObj()
  55
  56         RealSegments = []
  57         Direction = []
  58         flag = 0
  59         if not(args.__contains__('recursive')) :
  60                 Config.Count = 0
  61
  62         Config.Criterion = GetCriterion(ObjIDLists)
  63         for index, ObjList in enumerate(ObjIDLists) :
  64                 if not (ObjList[0] == -1 or Config.Count >= Config.Criterion):
  65                         if not(args.__contains__('recursive')) :
  66                                 Config.DirIndex = index
  67                                 if index > 1 : Config.RefPts = [Pt2, Pt3]
  68                                 elif index == 0 : Config.RefPts = [Pt1, Pt2]
  69                                 else : Config.RefPts = [Pt4, Pt3]
  70
  71                         if len(ObjList)>1 : flag = 1
  72                         else : flag = 0
  73                         for ObjID in ObjList:
  74                            ToLook0 = [2,3,0,1][index]
  75                            ToLook1 = [3,2,1,0][index]
  76                            CommonSide =  FindCommonSide(Config.ListObj[ObjID].DirBoundaries(ToLook1),dummy.DirBoundaries(ToLook0))
  77                            ToLook2 = [1,0,3,2][index]
  78                            RealSegments = Config.ListObj[ObjID].DirectionalMeshParams[ToLook2]*IntLen(CommonSide)/IntLen(Config.ListObj[ObjID].DirBoundaries(ToLook1))
  79                            LocalDelta = 1.*IntLen(CommonSide)/RealSegments
  80                            print "Direction:", ["West","East","South","North"][ToLook2]
  81                            print "IntLen(CommonSide):",IntLen(CommonSide)
  82                            print "RealSegments:",RealSegments
  83                            print "LocalDelta:",LocalDelta
  84                            if flag and Config.Count < Config.Criterion:
  85                                 if index ==0 :
  86                                         if abs(CommonSide[0] - Ymin)<1e-7 : SouthGR = GroupNames[0]
  87                                         else : SouthGR = None
  88                                         if abs(CommonSide[1] - Ymax)<1e-7 : NorthGR = GroupNames[1]
  89                                         else : NorthGR = None
  90
  91                                         NDelta = Config.ListObj[ObjID].DirectionalMeshParams[ToLook2]* (LocalDelta-GlobalDelta[Convention[index]])
  92                                         [Pt1,Pt2] = Config.RefPts
  93                                         Coef = [1.,-1.][index]
  94                                         Vref1 = [Coef*(Pt2[0]-Pt1[0]),Coef*(Pt2[1]-Pt1[1])]
  95                                         Vref2 = NormalizeVector([Pt2[0]-Pt3[0],Pt2[1]-Pt3[1]])
  96                                         Ptref = Config.ListObj[ObjID].PtCoor[[2,3][index]]
  97                                         NewPt = ExtrapPoint (Ptref,Vref1,Vref2,NDelta)
  98                                         CompositeBoxF (Pt1, Pt2, NewPt, Ptref, recursive=1, groups = [SouthGR,NorthGR]+GroupNames[2:4])
  99                                 elif index == 1:
 100                                         if abs(CommonSide[0] - Ymin)<1e-7 : SouthGR = GroupNames[0]
 101                                         else : SouthGR = None
 102                                         if abs(CommonSide[1] - Ymax)<1e-7 : NorthGR = GroupNames[1]
 103                                         else : NorthGR = None
 104
 105                                         NDelta = Config.ListObj[ObjID].DirectionalMeshParams[ToLook2]* (LocalDelta-GlobalDelta[Convention[index]])
 106                                         [Pt4,Pt3] = Config.RefPts
 107                                         Coef = 1.
 108                                         Vref1 = [Coef*(Pt4[0]-Pt3[0]),Coef*(Pt4[1]-Pt3[1])]
 109                                         Vref2 = NormalizeVector([Pt1[0]-Pt4[0],Pt1[1]-Pt4[1]])
 110                                         Ptref = Config.ListObj[ObjID].PtCoor[0]
 111                                         NewPt = ExtrapPoint (Ptref,Vref1,Vref2,NDelta)
 112                                         CompositeBoxF (NewPt, Ptref, Pt3, Pt4, recursive=1, groups = [SouthGR,NorthGR]+GroupNames[2:4])
 113                                 else :
 114                                         if abs(CommonSide[0] - Xmin)<1e-7 : WestGR = GroupNames[2]
 115                                         else : WestGR = None
 116                                         if abs(CommonSide[1] - Xmax)<1e-7 : EastGR = GroupNames[3]
 117                                         else : EastGR = None
 118
 119                                         NDelta = Config.ListObj[ObjID].DirectionalMeshParams[ToLook2]* (LocalDelta-GlobalDelta[Convention[index]])
 120                                         [Pt2,Pt3] = Config.RefPts
 121                                         Coef = [1.,-1.][index-2]
 122                                         Vref1 = [Coef*(Pt3[0]-Pt2[0]),Coef*(Pt3[1]-Pt2[1])]
 123                                         Vref2 = NormalizeVector([Pt3[0]-Pt4[0],Pt3[1]-Pt4[1]])
 124                                         Ptref = Config.ListObj[ObjID].PtCoor[[3,0][index-2]]
 125                                         NewPt = ExtrapPoint (Ptref,Vref1,Vref2,NDelta)
 126                                         CompositeBoxF (Ptref, Pt2, Pt3, NewPt, recursive=1, groups = GroupNames[0:2] + [WestGR,EastGR])
 127
 128                            if Config.Count >= Config.Criterion :
 129                                 break
 130         if flag == 0 and Config.Count < Config.Criterion:
 131                 print "Creating NonOrtho object with the points:", Pt1,Pt2,Pt3,Pt4
 132                 MacObject('NonOrtho',[Pt1,Pt2,Pt3,Pt4] ,['auto'], groups = GroupNames)
 133
 134                 Config.Count += 1
 135                 if Config.DirIndex > 1 : Config.RefPts = [Pt1, Pt4]
 136                 elif Config.DirIndex==0 : Config.RefPts = [Pt4, Pt3]
 137                 else : Config.RefPts = [Pt1, Pt2]
 138
 139 def FindCommonSide (Int1, Int2) :
 140         if max(Int1[0],Int2[0])<min(Int1[1],Int2[1]): return [max(Int1[0],Int2[0]), min(Int1[1],Int2[1])]
 141         else :
 142                 print "Can not find interval intersection, returning [0,0]..."
 143                 return [0,0]
 144
 145 def IntLen (Interval) :
 146         return float(abs(Interval[1]-Interval[0]))
 147
 148 def RemoveLastObj() :
 149         Config.ListObj = Config.ListObj[:-1]
 150         Config.Connections = Config.Connections[:-1]
 151
 152 def NormalizeVector(V):
 153         Magnitude = math.sqrt(GenFunctions.DotProd(V,V))
 154         return [ V[i]/Magnitude for i in range(len(V))]
 155
 156 def GetCriterion (ObjListIDs):
 157         return max(Config.Criterion, max(len(ObjListIDs[0]),len(ObjListIDs[1]))*max(len(ObjListIDs[2]),len(ObjListIDs[3])))
 158
 159 def SortObjLists (List,X0,Y0,DX,DY) :
 160         """
 161         This function sorts the list of neighbouring objects on each side, according to their intersection
 162         with the object being created. From South to North and from East to West
 163         """
 164         Output = List
 165         # First find the directions where no neighbour exists
 166         # Important : Here we assume that exactly two directions have no neighbours !!!
 167         #             Should we change this to allow a more general case ????
 168         dummy = IndexMultiOcc(List,(-1,))
 169
 170         # dummy[0] is either 0, meaning there is no neighbour on X- (West)
 171         #                 or 1, meaning there is no neighbour on X+ (East)
 172         # Similarly dummy[1] can be either 2 or 3 (South and North respectively)
 173         # In order to get back to the formalism of groups (SNWE)
 174         #  => we do the following to define Sense of no neighbours and then the Direction list
 175         #       is calculated as to include uniquely the directions where we DO have neighbours
 176         if len(dummy) == 1 :
 177                 # This adds a second direction where neighbours are not regarded, it is either 0 or 2
 178                 dummy.append(2*(dummy[0]+2<4))
 179                 print("Careful, you have neighbours on 3 or more sides of the box, we will not check if on two parallel sides the boxes are compatible !!!")
 180         if len(dummy) == 2 or len(dummy) == 1 :
 181                 # Sense contains : Vertical then Horizontal
 182                 Sense = [dummy[1]%2,dummy[0]]
 183                 DirList = [[1,0][dummy[0]],[3,2][dummy[1]%2]]
 184                 for index,Direction in enumerate(DirList) :
 185                         ObjList = List[Direction]
 186                         RankMin = []
 187                         ToLook0 = [2,2,0,0][Direction]
 188                         ToLook1 = [3,2,1,0][Direction]
 189                         for index1,ObjID in enumerate(ObjList) :
 190                                 RankMin.append([-1.,1.][Sense[index]] * FindCommonSide(Config.ListObj[ObjID].DirBoundaries(ToLook1),[X0-DX/2.,X0+DX/2.,Y0-DY/2.,Y0+DY/2.][ToLook0:ToLook0+2])[Sense[index]])
 191                         Output[Direction] = SortList(ObjList,RankMin)
 192
 193         elif len(dummy) == 3 :
 194                 # We find the direction where we do have neighbours and then we sort the object list along it
 195                 Sense = dummy[0]%2
 196                 Direction = [ i not in dummy for i in range(4) ].index(True)
 197                 ObjList = List[Direction]
 198                 RankMin = []
 199                 ToLook0 = [2,2,0,0][Direction]
 200                 ToLook1 = [3,2,1,0][Direction]
 201                 for index1,ObjID in enumerate(ObjList) :
 202                         RankMin.append([-1.,1.][Sense] * FindCommonSide(Config.ListObj[ObjID].DirBoundaries(ToLook1),[X0-DX/2.,X0+DX/2.,Y0-DY/2.,Y0+DY/2.][ToLook0:ToLook0+2])[Sense])
 203                 Output[Direction] = SortList(ObjList,RankMin)
 204         else :
 205                 print ("Error : the composite box being created has no neighbours, how on earth do you want us to inherit its mesh parameters!!!")
 206
 207
 208         return Output
 209
 210 def IndexMultiOcc (Array,Element) :
 211         """
 212         This functions returns the occurrences indices of Element in Array.
 213         As opposed to Array.index(Element) method, this allows determining
 214         multiple entries rather than just the first one!
 215         """
 216         Output = []
 217         try : Array.index(Element)
 218         except ValueError : print "No more occurrences"
 219         else : Output.append(Array.index(Element))
 220
 221         if not(Output == []) and len(Array) > 1 :
 222                 for index, ArrElem in enumerate(Array[Output[0]+1:]) :
 223                         if ArrElem == Element : Output.append(index+Output[0]+1)
 224
 225         return Output
 226
 227 def SortList (ValList, CritList):
 228         Output = []
 229         SortedCritList = copy.copy(CritList)
 230         SortedCritList.sort()
 231         for i in range(0,len(ValList)):
 232                 index = CritList.index(SortedCritList[i])
 233                 Output.append(ValList[index])
 234         return Output
 235
 236 def ExtrapPoint (Ptref,Vref1,Vref2,Delta):
 237         """
 238         This function allows determining the absolute coordinates of an extrapolation point
 239         as shown in the following :
 240
 241
 242         ExtrapPoint x---Vref2->--------o
 243                    /       delta_glob  |Vref1
 244                   /                    |
 245            Ptref x---------------------+
 246                         delta_loc * Nseg
 247
 248              Delta = (delta_loc - delta_glob) * Nseg
 249         """
 250
 251         X = Ptref[0] + Vref1[0] + Delta*Vref2[0]
 252         Y = Ptref[1] + Vref1[1] + Delta*Vref2[1]
 253         return (X,Y,)
 254