src/Tools/MacMesh/MacMesh/CompositeBoxF.py

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