source: issm/trunk/src/py3/parameterization/setflowequation.py@ 20500

Last change on this file since 20500 was 20500, checked in by Mathieu Morlighem, 9 years ago

merged trunk-jpl and trunk for revision 20497
File size: 13.5 KB
Line 
1import numpy
2from model import model
3from pairoptions import pairoptions
4import MatlabFuncs as m
5import PythonFuncs as p
6from FlagElements import FlagElements
7
8def setflowequation(md,**kwargs):
9 """
10 SETFLOWEQUATION - associate a solution type to each element
11
12 This routine works like plotmodel: it works with an even number of inputs
13 'SIA','SSA','HO','L1L2','FS' and 'fill' are the possible options
14 that must be followed by the corresponding exp file or flags list
15 It can either be a domain file (argus type, .exp extension), or an array of element flags.
16 If user wants every element outside the domain to be
17 setflowequationd, add '~' to the name of the domain file (ex: '~HO.exp');
18 an empty string '' will be considered as an empty domain
19 a string 'all' will be considered as the entire domain
20 You can specify the type of coupling, 'penalties' or 'tiling', to use with the input 'coupling'
21
22 Usage:
23 md=setflowequation(md,varargin)
24
25 Example:
26 md=setflowequation(md,'HO','HO.exp',fill','SIA','coupling','tiling');
27 """
28
29 #some checks on list of arguments
30 if not isinstance(md,model) or not len(kwargs):
31 raise TypeError("setflowequation error message")
32
33 #process options
34 options=pairoptions(**kwargs)
35 print(options)
36# options=deleteduplicates(options,1);
37
38 #Find_out what kind of coupling to use
39 coupling_method=options.getfieldvalue('coupling','tiling')
40 if coupling_method is not 'tiling' or not 'penalties':
41 raise TypeError("coupling type can only be: tiling or penalties")
42
43 #recover elements distribution
44 SIAflag = FlagElements(md,options.getfieldvalue('SIA',''))
45 SSAflag = FlagElements(md,options.getfieldvalue('SSA',''))
46 HOflag = FlagElements(md,options.getfieldvalue('HO',''))
47 L1L2flag = FlagElements(md,options.getfieldvalue('L1L2',''))
48 FSflag = FlagElements(md,options.getfieldvalue('FS',''))
49 filltype = options.getfieldvalue('fill','none')
50
51 #Flag the elements that have not been flagged as filltype
52 if filltype is 'SIA':
53 SIAflag[numpy.nonzero(numpy.logical_not(p.logical_or_n(SSAflag,HOflag)))]=True
54 elif filltype is 'SSA':
55 SSAflag[numpy.nonzero(numpy.logical_not(p.logical_or_n(SIAflag,HOflag,FSflag)))]=True
56 elif filltype is 'HO':
57 HOflag[numpy.nonzero(numpy.logical_not(p.logical_or_n(SIAflag,SSAflag,FSflag)))]=True
58
59 #check that each element has at least one flag
60 if not any(SIAflag+SSAflag+L1L2flag+HOflag+FSflag):
61 raise TypeError("elements type not assigned, supported models are 'SIA','SSA','HO' and 'FS'")
62
63 #check that each element has only one flag
64 if any(SIAflag+SSAflag+L1L2flag+HOflag+FSflag>1):
65 print("setflowequation warning message: some elements have several types, higher order type is used for them")
66 SIAflag[numpy.nonzero(numpy.logical_and(SIAflag,SSAflag))]=False
67 SIAflag[numpy.nonzero(numpy.logical_and(SIAflag,HOflag))]=False
68 SSAflag[numpy.nonzero(numpy.logical_and(SSAflag,HOflag))]=False
69
70 #FS can only be used alone for now:
71 if any(FSflag) and any(SIAflag):
72 raise TypeError("FS cannot be used with any other model for now, put FS everywhere")
73
74 #Initialize node fields
75 nodeonSIA=numpy.zeros(md.mesh.numberofvertices,bool)
76 nodeonSIA[md.mesh.elements[numpy.nonzero(SIAflag),:]-1]=True
77 nodeonSSA=numpy.zeros(md.mesh.numberofvertices,bool)
78 nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True
79 nodeonL1L2=numpy.zeros(md.mesh.numberofvertices,bool)
80 nodeonL1L2[md.mesh.elements[numpy.nonzero(L1L2flag),:]-1]=True
81 nodeonHO=numpy.zeros(md.mesh.numberofvertices,bool)
82 nodeonHO[md.mesh.elements[numpy.nonzero(HOflag),:]-1]=True
83 nodeonFS=numpy.zeros(md.mesh.numberofvertices,bool)
84 noneflag=numpy.zeros(md.mesh.numberofelements,bool)
85
86 #First modify FSflag to get rid of elements contrained everywhere (spc + border with HO or SSA)
87 if any(FSflag):
88# fullspcnodes=double((~isnan(md.stressbalance.spcvx)+~isnan(md.stressbalance.spcvy)+~isnan(md.stressbalance.spcvz))==3 | (nodeonHO & nodeonFS)); %find all the nodes on the boundary of the domain without icefront
89 fullspcnodes=numpy.logical_or(numpy.logical_not(numpy.isnan(md.stressbalance.spcvx)).astype(int)+ \
90 numpy.logical_not(numpy.isnan(md.stressbalance.spcvy)).astype(int)+ \
91 numpy.logical_not(numpy.isnan(md.stressbalance.spcvz)).astype(int)==3, \
92 numpy.logical_and(nodeonHO,nodeonFS)).astype(int) #find all the nodes on the boundary of the domain without icefront
93# fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront
94 fullspcelems=(numpy.sum(fullspcnodes[md.mesh.elements-1],axis=1)==6).astype(int) #find all the nodes on the boundary of the domain without icefront
95 FSflag[numpy.nonzero(fullspcelems.reshape(-1))]=False
96 nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True
97
98 #Then complete with NoneApproximation or the other model used if there is no FS
99 if any(FSflag):
100 if any(HOflag): #fill with HO
101 HOflag[numpy.logical_not(FSflag)]=True
102 nodeonHO[md.mesh.elements[numpy.nonzero(HOflag),:]-1]=True
103 elif any(SSAflag): #fill with SSA
104 SSAflag[numpy.logical_not(FSflag)]=True
105 nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True
106 else: #fill with none
107 noneflag[numpy.nonzero(numpy.logical_not(FSflag))]=True
108
109 #Now take care of the coupling between SSA and HO
110 md.stressbalance.vertex_pairing=numpy.array([])
111 nodeonSSAHO=numpy.zeros(md.mesh.numberofvertices,bool)
112 nodeonHOFS=numpy.zeros(md.mesh.numberofvertices,bool)
113 nodeonSSAFS=numpy.zeros(md.mesh.numberofvertices,bool)
114 SSAHOflag=numpy.zeros(md.mesh.numberofelements,bool)
115 SSAFSflag=numpy.zeros(md.mesh.numberofelements,bool)
116 HOFSflag=numpy.zeros(md.mesh.numberofelements,bool)
117 if coupling_method is 'penalties':
118 #Create the border nodes between HO and SSA and extrude them
119 numnodes2d=md.mesh.numberofvertices2d
120 numlayers=md.mesh.numberoflayers
121 bordernodes2d=numpy.nonzero(numpy.logical_and(nodeonHO[0:numnodes2d],nodeonSSA[0:numnodes2d]))[0]+1 #Nodes connected to two different types of elements
122
123 #initialize and fill in penalties structure
124 if numpy.all(numpy.logical_not(numpy.isnan(bordernodes2d))):
125 penalties=numpy.zeros((0,2))
126 for i in range(1,numlayers):
127 penalties=numpy.vstack((penalties,numpy.hstack((bordernodes2d.reshape(-1,1),bordernodes2d.reshape(-1,1)+md.mesh.numberofvertices2d*(i)))))
128 md.stressbalance.vertex_pairing=penalties
129
130 elif coupling_method is 'tiling':
131 if any(SSAflag) and any(HOflag): #coupling SSA HO
132 #Find node at the border
133 nodeonSSAHO[numpy.nonzero(numpy.logical_and(nodeonSSA,nodeonHO))]=True
134 #SSA elements in contact with this layer become SSAHO elements
135 matrixelements=m.ismember(md.mesh.elements-1,numpy.nonzero(nodeonSSAHO)[0])
136 commonelements=numpy.sum(matrixelements,axis=1)!=0
137 commonelements[numpy.nonzero(HOflag)]=False #only one layer: the elements previously in SSA
138 SSAflag[numpy.nonzero(commonelements)]=False #these elements are now SSAHOelements
139 SSAHOflag[numpy.nonzero(commonelements)]=True
140 nodeonSSA[:]=False
141 nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True
142
143 #rule out elements that don't touch the 2 boundaries
144 pos=numpy.nonzero(SSAHOflag)[0]
145 elist=numpy.zeros(numpy.size(pos),dtype=int)
146 elist = elist + numpy.sum(nodeonSSA[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
147 elist = elist - numpy.sum(nodeonHO[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool)
148 pos1=numpy.nonzero(elist==1)[0]
149 SSAflag[pos[pos1]]=True
150 SSAHOflag[pos[pos1]]=False
151 pos2=numpy.nonzero(elist==-1)[0]
152 HOflag[pos[pos2]]=True
153 SSAHOflag[pos[pos2]]=False
154
155 #Recompute nodes associated to these elements
156 nodeonSSA[:]=False
157 nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True
158 nodeonHO[:]=False
159 nodeonHO[md.mesh.elements[numpy.nonzero(HOflag),:]-1]=True
160 nodeonSSAHO[:]=False
161 nodeonSSAHO[md.mesh.elements[numpy.nonzero(SSAHOflag),:]-1]=True
162
163 elif any(HOflag) and any(FSflag): #coupling HO FS
164 #Find node at the border
165 nodeonHOFS[numpy.nonzero(numpy.logical_and(nodeonHO,nodeonFS))]=True
166 #FS elements in contact with this layer become HOFS elements
167 matrixelements=m.ismember(md.mesh.elements-1,numpy.nonzero(nodeonHOFS)[0])
168 commonelements=numpy.sum(matrixelements,axis=1)!=0
169 commonelements[numpy.nonzero(HOflag)]=False #only one layer: the elements previously in SSA
170 FSflag[numpy.nonzero(commonelements)]=False #these elements are now SSAHOelements
171 HOFSflag[numpy.nonzero(commonelements)]=True
172 nodeonFS=numpy.zeros(md.mesh.numberofvertices,bool)
173 nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True
174
175 #rule out elements that don't touch the 2 boundaries
176 pos=numpy.nonzero(HOFSflag)[0]
177 elist=numpy.zeros(numpy.size(pos),dtype=int)
178 elist = elist + numpy.sum(nodeonFS[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
179 elist = elist - numpy.sum(nodeonHO[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
180 pos1=numpy.nonzero(elist==1)[0]
181 FSflag[pos[pos1]]=True
182 HOFSflag[pos[pos1]]=False
183 pos2=numpy.nonzero(elist==-1)[0]
184 HOflag[pos[pos2]]=True
185 HOFSflag[pos[pos2]]=False
186
187 #Recompute nodes associated to these elements
188 nodeonFS[:]=False
189 nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True
190 nodeonHO[:]=False
191 nodeonHO[md.mesh.elements[numpy.nonzero(HOflag),:]-1]=True
192 nodeonHOFS[:]=False
193 nodeonHOFS[md.mesh.elements[numpy.nonzero(HOFSflag),:]-1]=True
194
195 elif any(FSflag) and any(SSAflag):
196 #Find node at the border
197 nodeonSSAFS[numpy.nonzero(numpy.logical_and(nodeonSSA,nodeonFS))]=True
198 #FS elements in contact with this layer become SSAFS elements
199 matrixelements=m.ismember(md.mesh.elements-1,numpy.nonzero(nodeonSSAFS)[0])
200 commonelements=numpy.sum(matrixelements,axis=1)!=0
201 commonelements[numpy.nonzero(SSAflag)]=False #only one layer: the elements previously in SSA
202 FSflag[numpy.nonzero(commonelements)]=False #these elements are now SSASSAelements
203 SSAFSflag[numpy.nonzero(commonelements)]=True
204 nodeonFS=numpy.zeros(md.mesh.numberofvertices,bool)
205 nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True
206
207 #rule out elements that don't touch the 2 boundaries
208 pos=numpy.nonzero(SSAFSflag)[0]
209 elist=numpy.zeros(numpy.size(pos),dtype=int)
210 elist = elist + numpy.sum(nodeonSSA[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
211 elist = elist - numpy.sum(nodeonFS[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool)
212 pos1=numpy.nonzero(elist==1)[0]
213 SSAflag[pos[pos1]]=True
214 SSAFSflag[pos[pos1]]=False
215 pos2=numpy.nonzero(elist==-1)[0]
216 FSflag[pos[pos2]]=True
217 SSAFSflag[pos[pos2]]=False
218
219 #Recompute nodes associated to these elements
220 nodeonSSA[:]=False
221 nodeonSSA[md.mesh.elements[numpy.nonzero(SSAflag),:]-1]=True
222 nodeonFS[:]=False
223 nodeonFS[md.mesh.elements[numpy.nonzero(FSflag),:]-1]=True
224 nodeonSSAFS[:]=False
225 nodeonSSAFS[md.mesh.elements[numpy.nonzero(SSAFSflag),:]-1]=True
226
227 elif any(FSflag) and any(SIAflag):
228 raise TypeError("type of coupling not supported yet")
229
230 #Create SSAHOApproximation where needed
231 md.flowequation.element_equation=numpy.zeros(md.mesh.numberofelements,int)
232 md.flowequation.element_equation[numpy.nonzero(noneflag)]=0
233 md.flowequation.element_equation[numpy.nonzero(SIAflag)]=1
234 md.flowequation.element_equation[numpy.nonzero(SSAflag)]=2
235 md.flowequation.element_equation[numpy.nonzero(L1L2flag)]=3
236 md.flowequation.element_equation[numpy.nonzero(HOflag)]=4
237 md.flowequation.element_equation[numpy.nonzero(FSflag)]=5
238 md.flowequation.element_equation[numpy.nonzero(SSAHOflag)]=6
239 md.flowequation.element_equation[numpy.nonzero(SSAFSflag)]=7
240 md.flowequation.element_equation[numpy.nonzero(HOFSflag)]=8
241
242 #border
243 md.flowequation.borderHO=nodeonHO
244 md.flowequation.borderSSA=nodeonSSA
245 md.flowequation.borderFS=nodeonFS
246
247 #Create vertices_type
248 md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices,int)
249 pos=numpy.nonzero(nodeonSSA)
250 md.flowequation.vertex_equation[pos]=2
251 pos=numpy.nonzero(nodeonL1L2)
252 md.flowequation.vertex_equation[pos]=3
253 pos=numpy.nonzero(nodeonHO)
254 md.flowequation.vertex_equation[pos]=4
255 pos=numpy.nonzero(nodeonFS)
256 md.flowequation.vertex_equation[pos]=5
257 #DO SIA LAST! Otherwise spcs might not be set up correctly (SIA should have priority)
258 pos=numpy.nonzero(nodeonSIA)
259 md.flowequation.vertex_equation[pos]=1
260 if any(FSflag):
261 pos=numpy.nonzero(numpy.logical_not(nodeonFS))
262 if not (any(HOflag) or any(SSAflag)):
263 md.flowequation.vertex_equation[pos]=0
264 pos=numpy.nonzero(nodeonSSAHO)
265 md.flowequation.vertex_equation[pos]=6
266 pos=numpy.nonzero(nodeonHOFS)
267 md.flowequation.vertex_equation[pos]=7
268 pos=numpy.nonzero(nodeonSSAFS)
269 md.flowequation.vertex_equation[pos]=8
270
271 #figure out solution types
272 md.flowequation.isSIA=any(md.flowequation.element_equation==1)
273 md.flowequation.isSSA=any(md.flowequation.element_equation==2)
274 md.flowequation.isL1L2=any(md.flowequation.element_equation==3)
275 md.flowequation.isHO=any(md.flowequation.element_equation==4)
276 md.flowequation.isFS=any(md.flowequation.element_equation==5)
277
278 return md
279
280 #Check that tiling can work:
281 if any(md.flowequation.borderSSA) and any(md.flowequation.borderHO) and any(md.flowequation.borderHO + md.flowequation.borderSSA !=1):
282 raise TypeError("error coupling domain too irregular")
283 if any(md.flowequation.borderSSA) and any(md.flowequation.borderFS) and any(md.flowequation.borderFS + md.flowequation.borderSSA !=1):
284 raise TypeError("error coupling domain too irregular")
285 if any(md.flowequation.borderFS) and any(md.flowequation.borderHO) and any(md.flowequation.borderHO + md.flowequation.borderFS !=1):
286 raise TypeError("error coupling domain too irregular")
287
288 return md
289
Note: See TracBrowser for help on using the repository browser.