src/Tools/MacMesh/MacMesh/CompositeBoxF.py

   1 # Copyright (C) 2014-2016  EDF R&D
   2 #
   3 # This library is free software; you can redistribute it and/or
   4 # modify it under the terms of the GNU Lesser General Public
   5 # License as published by the Free Software Foundation; either
   6 # version 2.1 of the License, or (at your option) any later version.
   7 #
   8 # This library is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 # Lesser General Public License for more details.
  12 #
  13 # You should have received a copy of the GNU Lesser General Public
  14 # License along with this library; if not, write to the Free Software
  15 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 #
  17 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 #
  19
  20
  21 # INTRODUCTION HERE
  22
  23 import sys, math, copy, subprocess
  24 CWD = subprocess.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 = sorted(copy.copy(CritList))
 226     for i in range(0,len(ValList)):
 227         index = CritList.index(SortedCritList[i])
 228         Output.append(ValList[index])
 229     return Output
 230
 231 def ExtrapPoint (Ptref,Vref1,Vref2,Delta):
 232     """
 233     This function allows determining the absolute coordinates of an extrapolation point
 234     as shown in the following :
 235
 236
 237     ExtrapPoint x---Vref2->--------o
 238                /       delta_glob  |Vref1
 239               /                    |
 240        Ptref x---------------------+
 241                     delta_loc * Nseg
 242
 243          Delta = (delta_loc - delta_glob) * Nseg
 244     """
 245
 246     X = Ptref[0] + Vref1[0] + Delta*Vref2[0]
 247     Y = Ptref[1] + Vref1[1] + Delta*Vref2[1]
 248     return (X,Y,)