1 # Copyright (C) 2007-2014 CEA/DEN, EDF R&D, OPEN CASCADE
3 # Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 # CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
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.
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.
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
20 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
27 import sys, math, copy, commands
28 CWD = commands.getoutput('pwd')
31 from MacObject import *
32 import Config, GenFunctions
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
59 if not(args.__contains__('recursive')) :
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]
71 if len(ObjList)>1 : flag = 1
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:
86 if abs(CommonSide[0] - Ymin)<1e-7 : SouthGR = GroupNames[0]
88 if abs(CommonSide[1] - Ymax)<1e-7 : NorthGR = GroupNames[1]
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])
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
105 NDelta = Config.ListObj[ObjID].DirectionalMeshParams[ToLook2]* (LocalDelta-GlobalDelta[Convention[index]])
106 [Pt4,Pt3] = Config.RefPts
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])
114 if abs(CommonSide[0] - Xmin)<1e-7 : WestGR = GroupNames[2]
116 if abs(CommonSide[1] - Xmax)<1e-7 : EastGR = GroupNames[3]
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])
128 if Config.Count >= Config.Criterion :
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)
135 if Config.DirIndex > 1 : Config.RefPts = [Pt1, Pt4]
136 elif Config.DirIndex==0 : Config.RefPts = [Pt4, Pt3]
137 else : Config.RefPts = [Pt1, Pt2]
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])]
142 print "Can not find interval intersection, returning [0,0]..."
145 def IntLen (Interval) :
146 return float(abs(Interval[1]-Interval[0]))
148 def RemoveLastObj() :
149 Config.ListObj = Config.ListObj[:-1]
150 Config.Connections = Config.Connections[:-1]
152 def NormalizeVector(V):
153 Magnitude = math.sqrt(GenFunctions.DotProd(V,V))
154 return [ V[i]/Magnitude for i in range(len(V))]
156 def GetCriterion (ObjListIDs):
157 return max(Config.Criterion, max(len(ObjListIDs[0]),len(ObjListIDs[1]))*max(len(ObjListIDs[2]),len(ObjListIDs[3])))
159 def SortObjLists (List,X0,Y0,DX,DY) :
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
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,))
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
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]
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)
193 elif len(dummy) == 3 :
194 # We find the direction where we do have neighbours and then we sort the object list along it
196 Direction = [ i not in dummy for i in range(4) ].index(True)
197 ObjList = List[Direction]
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)
205 print ("Error : the composite box being created has no neighbours, how on earth do you want us to inherit its mesh parameters!!!")
210 def IndexMultiOcc (Array,Element) :
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!
217 try : Array.index(Element)
218 except ValueError : print "No more occurrences"
219 else : Output.append(Array.index(Element))
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)
227 def SortList (ValList, CritList):
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])
236 def ExtrapPoint (Ptref,Vref1,Vref2,Delta):
238 This function allows determining the absolute coordinates of an extrapolation point
239 as shown in the following :
242 ExtrapPoint x---Vref2->--------o
245 Ptref x---------------------+
248 Delta = (delta_loc - delta_glob) * Nseg
251 X = Ptref[0] + Vref1[0] + Delta*Vref2[0]
252 Y = Ptref[1] + Vref1[1] + Delta*Vref2[1]