1 # Copyright (C) 2014-2022 EDF R&D
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.
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.
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
17 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 import sys, math, copy, subprocess
24 CWD = subprocess.getoutput('pwd')
27 from MacObject import *
28 import Config, GenFunctions
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
55 if not(args.__contains__('recursive')) :
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]
67 if len(ObjList)>1 : flag = 1
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:
82 if abs(CommonSide[0] - Ymin)<1e-7 : SouthGR = GroupNames[0]
84 if abs(CommonSide[1] - Ymax)<1e-7 : NorthGR = GroupNames[1]
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])
96 if abs(CommonSide[0] - Ymin)<1e-7 : SouthGR = GroupNames[0]
98 if abs(CommonSide[1] - Ymax)<1e-7 : NorthGR = GroupNames[1]
101 NDelta = Config.ListObj[ObjID].DirectionalMeshParams[ToLook2]* (LocalDelta-GlobalDelta[Convention[index]])
102 [Pt4,Pt3] = Config.RefPts
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])
110 if abs(CommonSide[0] - Xmin)<1e-7 : WestGR = GroupNames[2]
112 if abs(CommonSide[1] - Xmax)<1e-7 : EastGR = GroupNames[3]
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])
124 if Config.Count >= Config.Criterion :
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)
131 if Config.DirIndex > 1 : Config.RefPts = [Pt1, Pt4]
132 elif Config.DirIndex==0 : Config.RefPts = [Pt4, Pt3]
133 else : Config.RefPts = [Pt1, Pt2]
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])]
138 print("Can not find interval intersection, returning [0,0]...")
141 def IntLen (Interval) :
142 return float(abs(Interval[1]-Interval[0]))
144 def RemoveLastObj() :
145 Config.ListObj = Config.ListObj[:-1]
146 Config.Connections = Config.Connections[:-1]
148 def NormalizeVector(V):
149 Magnitude = math.sqrt(GenFunctions.DotProd(V,V))
150 return [ V[i]/Magnitude for i in range(len(V))]
152 def GetCriterion (ObjListIDs):
153 return max(Config.Criterion, max(len(ObjListIDs[0]),len(ObjListIDs[1]))*max(len(ObjListIDs[2]),len(ObjListIDs[3])))
155 def SortObjLists (List,X0,Y0,DX,DY) :
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
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,))
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
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]
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)
189 elif len(dummy) == 3 :
190 # We find the direction where we do have neighbours and then we sort the object list along it
192 Direction = [ i not in dummy for i in range(4) ].index(True)
193 ObjList = List[Direction]
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)
201 print ("Error : the composite box being created has no neighbours, how on earth do you want us to inherit its mesh parameters!!!")
206 def IndexMultiOcc (Array,Element) :
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!
213 try : Array.index(Element)
214 except ValueError : print("No more occurrences")
215 else : Output.append(Array.index(Element))
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)
223 def SortList (ValList, CritList):
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])
231 def ExtrapPoint (Ptref,Vref1,Vref2,Delta):
233 This function allows determining the absolute coordinates of an extrapolation point
234 as shown in the following :
237 ExtrapPoint x---Vref2->--------o
240 Ptref x---------------------+
243 Delta = (delta_loc - delta_glob) * Nseg
246 X = Ptref[0] + Vref1[0] + Delta*Vref2[0]
247 Y = Ptref[1] + Vref1[1] + Delta*Vref2[1]