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