source: issm/oecreview/Archive/25834-26739/ISSM-26554-26555.diff@ 26740

Last change on this file since 26740 was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739
File size: 41.5 KB
RevLine 
[26740]1Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26554)
4+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26555)
5@@ -20,7 +20,6 @@
6 //#define FSANALYTICAL 10
7 //#define LATERALFRICTION 1
8 //#define DISCSLOPE 1 //testing for SSA
9-#define NOSPCSHEARVEL 1 //MLHO
10
11 /*Model processing*/
12 void StressbalanceAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
13@@ -195,16 +194,10 @@
14 }
15 }
16 else{//MLHO
17- #ifdef NOSPCSHEARVEL
18- IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0);
19- IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2);
20- #else
21- /*Default: apply spcs to shear vx and shear vy*/
22- IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0);
23- IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,1);
24- IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2);
25- IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,3);
26- #endif
27+ IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx_base",StressbalanceAnalysisEnum,finiteelement,0);
28+ IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx_shear",StressbalanceAnalysisEnum,finiteelement,1);
29+ IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy_base",StressbalanceAnalysisEnum,finiteelement,2);
30+ IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy_shear",StressbalanceAnalysisEnum,finiteelement,3);
31 }
32 }
33
34Index: ../trunk-jpl/src/m/boundaryconditions/SetMLHOBC.m
35===================================================================
36--- ../trunk-jpl/src/m/boundaryconditions/SetMLHOBC.m (nonexistent)
37+++ ../trunk-jpl/src/m/boundaryconditions/SetMLHOBC.m (revision 26555)
38@@ -0,0 +1,16 @@
39+function md=SetMLHOBC(md)
40+%SETMLHOBC - Create the boundary conditions for stressbalance for MLHO: VxBase, VyBase, VxShear, VyShear
41+%
42+% Usage:
43+% md=SetMLHOBC(md)
44+%
45+
46+
47+%node on Dirichlet
48+if md.flowequation.isMLHO
49+ md.stressbalance.spcvx_base=md.stressbalance.spcvx;
50+ md.stressbalance.spcvy_base=md.stressbalance.spcvy;
51+
52+ md.stressbalance.spcvx_shear=NaN*ones(size(md.stressbalance.spcvx_base));
53+ md.stressbalance.spcvy_shear=NaN*ones(size(md.stressbalance.spcvy_base));
54+end
55Index: ../trunk-jpl/src/m/classes/flowequation.m
56===================================================================
57--- ../trunk-jpl/src/m/classes/flowequation.m (revision 26554)
58+++ ../trunk-jpl/src/m/classes/flowequation.m (revision 26555)
59@@ -61,10 +61,9 @@
60 if isfield(objstruct,'borderstokes'), self.borderFS = objstruct.borderstokes; end;
61 end
62
63- %Oct 8 2020
64- %Warning if MLHO
65+ %Nov 6 2021
66 if any(self.vertex_equation==4)
67- warning(['Monolayer Higher-Order (MLHO) detected in md.flowequation, this is probably a mistake, you need to run setflowequation again']);
68+ disp(['Monolayer Higher-Order (MLHO) detected in md.flowequation, this is still under development. Please double check your settings.']);
69 end
70
71 end% }}}
72Index: ../trunk-jpl/src/m/classes/stressbalance.m
73===================================================================
74--- ../trunk-jpl/src/m/classes/stressbalance.m (revision 26554)
75+++ ../trunk-jpl/src/m/classes/stressbalance.m (revision 26555)
76@@ -7,6 +7,10 @@
77 properties (SetAccess=public)
78 spcvx = NaN;
79 spcvy = NaN;
80+ spcvx_base = NaN;
81+ spcvy_base = NaN;
82+ spcvx_shear = NaN;
83+ spcvy_shear = NaN;
84 spcvz = NaN;
85 restol = 0;
86 reltol = 0;
87@@ -31,6 +35,13 @@
88 self.referential=project3d(md,'vector',self.referential,'type','node');
89 self.loadingforce=project3d(md,'vector',self.loadingforce,'type','node');
90
91+ % for MLHO
92+ if md.flowequation.isMLHO
93+ self.spcvx_base=project3d(md,'vector',self.spcvx_base,'type','node');
94+ self.spcvy_base=project3d(md,'vector',self.spcvy_base,'type','node');
95+ self.spcvx_shear=project3d(md,'vector',self.spcvx_shear,'type','poly','degree',4);
96+ self.spcvy_shear=project3d(md,'vector',self.spcvy_shear,'type','poly','degree',4);
97+ end
98 end % }}}
99 function self = stressbalance(varargin) % {{{
100 switch nargin
101@@ -122,6 +133,13 @@
102 end
103 md = checkfield(md,'fieldname','stressbalance.FSreconditioning','>',0);
104 end
105+ % CHECK THIS ONLY WORKS FOR MLHO
106+ if md.flowequation.isMLHO
107+ md = checkfield(md,'fieldname','stressbalance.spcvx_base','Inf',1,'timeseries',1);
108+ md = checkfield(md,'fieldname','stressbalance.spcvy_base','Inf',1,'timeseries',1);
109+ md = checkfield(md,'fieldname','stressbalance.spcvx_shear','Inf',1,'timeseries',1);
110+ md = checkfield(md,'fieldname','stressbalance.spcvy_shear','Inf',1,'timeseries',1);
111+ end
112 end % }}}
113 function list=defaultoutputs(self,md) % {{{
114
115@@ -150,6 +168,12 @@
116 fielddisplay(self,'spcvy','y-axis velocity constraint (NaN means no constraint) [m/yr]');
117 fielddisplay(self,'spcvz','z-axis velocity constraint (NaN means no constraint) [m/yr]');
118
119+ disp(sprintf('\n %s','MLHO boundary conditions:'));
120+ fielddisplay(self,'spcvx_base','x-axis basal velocity constraint (NaN means no constraint) [m/yr]');
121+ fielddisplay(self,'spcvy_base','y-axis basal velocity constraint (NaN means no constraint) [m/yr]');
122+ fielddisplay(self,'spcvx_shear','x-axis shear velocity constraint (NaN means no constraint) [m/yr]');
123+ fielddisplay(self,'spcvy_shear','y-axis shear velocity constraint (NaN means no constraint) [m/yr]');
124+
125 disp(sprintf('\n %s','Rift options:'));
126 fielddisplay(self,'rift_penalty_threshold','threshold for instability of mechanical constraints');
127 fielddisplay(self,'rift_penalty_lock','number of iterations before rift penalties are locked');
128@@ -201,6 +225,13 @@
129 outputs = [outputs defaultoutputs(self,md)]; %add defaults
130 end
131 WriteData(fid,prefix,'data',outputs,'name','md.stressbalance.requested_outputs','format','StringArray');
132+ % for MLHO
133+ if (md.flowequation.isMLHO)
134+ WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','spcvx_base','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
135+ WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','spcvy_base','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
136+ WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','spcvx_shear','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
137+ WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','spcvy_shear','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
138+ end
139 end % }}}
140 function savemodeljs(self,fid,modelname) % {{{
141
142@@ -222,6 +253,10 @@
143 writejs2Darray(fid,[modelname '.stressbalance.loadingforce'],self.loadingforce);
144 writejscellstring(fid,[modelname '.stressbalance.requested_outputs'],self.requested_outputs);
145
146+ writejs1Darray(fid,[modelname '.stressbalance.spcvx_base'],self.spcvx_shear);
147+ writejs1Darray(fid,[modelname '.stressbalance.spcvy_base'],self.spcvy_shear);
148+ writejs1Darray(fid,[modelname '.stressbalance.spcvx_shear'],self.spcvx_shear);
149+ writejs1Darray(fid,[modelname '.stressbalance.spcvy_shear'],self.spcvy_shear);
150 end % }}}
151 end
152 end
153Index: ../trunk-jpl/src/m/classes/stressbalance.py
154===================================================================
155--- ../trunk-jpl/src/m/classes/stressbalance.py (revision 26554)
156+++ ../trunk-jpl/src/m/classes/stressbalance.py (revision 26555)
157@@ -19,6 +19,10 @@
158 def __init__(self): # {{{
159 self.spcvx = np.nan
160 self.spcvy = np.nan
161+ self.spcvx_base = np.nan
162+ self.spcvy_base = np.nan
163+ self.spcvx_shear = np.nan
164+ self.spcvy_shear = np.nan
165 self.spcvz = np.nan
166 self.restol = 0
167 self.reltol = 0
168@@ -54,6 +58,11 @@
169 s += '{}\n'.format(fielddisplay(self, 'spcvy', 'y-axis velocity constraint (NaN means no constraint) [m / yr]'))
170 s += '{}\n'.format(fielddisplay(self, 'spcvz', 'z-axis velocity constraint (NaN means no constraint) [m / yr]'))
171 s += '{}\n'.format(fielddisplay(self, 'icefront', 'segments on ice front list (last column 0: Air, 1: Water, 2: Ice'))
172+ s += ' MLHO boundary conditions:\n'
173+ s += '{}\n'.format(fielddisplay(self, 'spcvx_base', 'x-axis basal velocity constraint (NaN means no constraint) [m / yr]'))
174+ s += '{}\n'.format(fielddisplay(self, 'spcvy_base', 'y-axis basal velocity constraint (NaN means no constraint) [m / yr]'))
175+ s += '{}\n'.format(fielddisplay(self, 'spcvx_shear', 'x-axis shear velocity constraint (NaN means no constraint) [m / yr]'))
176+ s += '{}\n'.format(fielddisplay(self, 'spcvy_shear', 'y-axis shear velocity constraint (NaN means no constraint) [m / yr]'))
177 s += ' Rift options:\n'
178 s += '{}\n'.format(fielddisplay(self, 'rift_penalty_threshold', 'threshold for instability of mechanical constraints'))
179 s += '{}\n'.format(fielddisplay(self, 'rift_penalty_lock', 'number of iterations before rift penalties are locked'))
180@@ -76,6 +85,12 @@
181 self.referential = project3d(md, 'vector', self.referential, 'type', 'node')
182 self.loadingforce = project3d(md, 'vector', self.loadingforce, 'type', 'node')
183
184+ if md.flowequation.isMLHO:
185+ self.spcvx_base = project3d(md, 'vector', self.spcvx_base, 'type', 'node')
186+ self.spcvy_base = project3d(md, 'vector', self.spcvy_base, 'type', 'node')
187+ self.spcvx_shear = project3d(md, 'vector', self.spcvx_shear, 'type', 'poly', 'degree', 4)
188+ self.spcvy_shear = project3d(md, 'vector', self.spcvy_shear, 'type', 'poly', 'degree', 4)
189+
190 return self
191 # }}}
192
193@@ -160,6 +175,11 @@
194 pos = np.nonzero(np.logical_and(md.mask.ocean_levelset, md.mesh.vertexonbase))
195 if np.any(np.logical_not(np.isnan(md.stressbalance.referential[pos, :]))):
196 md.checkmessage("no referential should be specified for basal vertices of grounded ice")
197+ if md.flowequation.isMLHO:
198+ md = checkfield(md, 'fieldname', 'stressbalance.spcvx_base', 'Inf', 1, 'timeseries', 1)
199+ md = checkfield(md, 'fieldname', 'stressbalance.spcvy_base', 'Inf', 1, 'timeseries', 1)
200+ md = checkfield(md, 'fieldname', 'stressbalance.spcvx_shear', 'Inf', 1, 'timeseries', 1)
201+ md = checkfield(md, 'fieldname', 'stressbalance.spcvy_shear', 'Inf', 1, 'timeseries', 1)
202 return md
203 # }}}
204
205@@ -194,4 +214,10 @@
206 outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
207 outputs = outputscopy
208 WriteData(fid, prefix, 'data', outputs, 'name', 'md.stressbalance.requested_outputs', 'format', 'StringArray')
209+ # MLHO
210+ if md.flowequation.isMLHO:
211+ WriteData(fid, prefix, 'object', self, 'class', 'stressbalance', 'fieldname', 'spcvx_base', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
212+ WriteData(fid, prefix, 'object', self, 'class', 'stressbalance', 'fieldname', 'spcvy_base', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
213+ WriteData(fid, prefix, 'object', self, 'class', 'stressbalance', 'fieldname', 'spcvx_shear', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
214+ WriteData(fid, prefix, 'object', self, 'class', 'stressbalance', 'fieldname', 'spcvy_shear', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
215 # }}}
216Index: ../trunk-jpl/src/m/extrusion/project3d.m
217===================================================================
218--- ../trunk-jpl/src/m/extrusion/project3d.m (revision 26554)
219+++ ../trunk-jpl/src/m/extrusion/project3d.m (revision 26555)
220@@ -6,11 +6,12 @@
221 % element vector of size (md.mesh.numberofelements2d,N/A).
222 % arguments:
223 % 'vector': 2d vector
224-% 'type': 'element' or 'node'.
225+% 'type': 'element' or 'node' or 'poly'
226 % options:
227 % 'layer' a layer number where vector should keep its values. If not specified, all layers adopt the
228 % value of the 2d vector.
229 % 'padding': default to 0 (value adopted by other 3d layers not being projected0
230+% 'degree': degree of polynomials when extrude from bottom to the top
231 %
232 % Egs:
233 % extruded_vector=project3d(md,'vector',vector2d,'type','node','layer',1,'padding',NaN);
234@@ -32,6 +33,7 @@
235 type = getfieldvalue(options,'type'); %mandatory
236 layer = getfieldvalue(options,'layer',0); %optional (do all layers otherwise)
237 paddingvalue = getfieldvalue(options,'padding',0); %0 by default
238+polyexponent = getfieldvalue(options,'degree',0); %0 by default, 0-degree polynomial
239
240 if length(vector2d)==1,
241 projected_vector=vector2d;
242@@ -79,6 +81,29 @@
243 projected_vector( ((layer-1)*md.mesh.numberofelements2d+1):(layer*md.mesh.numberofelements2d),:)=vector2d;
244 end
245
246+elseif strcmpi(type,'poly'), % interpolate values from 0 to 1 with a polynomial degree n
247+ %Initialize 3d vector
248+ if size(vector2d,1)==md.mesh.numberofvertices2d
249+ projected_vector=paddingvalue*ones(md.mesh.numberofvertices, size(vector2d,2));
250+ elseif size(vector2d,1)==md.mesh.numberofvertices2d+1
251+ projected_vector=paddingvalue*ones(md.mesh.numberofvertices+1,size(vector2d,2));
252+ projected_vector(end,:)=vector2d(end,:);
253+ vector2d=vector2d(1:end-1,:);
254+ else
255+ error('vector length not supported')
256+ end
257+
258+ polycoeff = [0:1./(md.mesh.numberoflayers-1):1];
259+
260+ %Fill in
261+ if layer==0,
262+ for i=1:md.mesh.numberoflayers,
263+ projected_vector(((i-1)*md.mesh.numberofvertices2d+1):(i*md.mesh.numberofvertices2d),:)=vector2d*(1-(1-polycoeff(i)).^polyexponent);
264+ end
265+ else
266+ projected_vector(((layer-1)*md.mesh.numberofvertices2d+1):(layer*md.mesh.numberofvertices2d),:)=vector2d*(1-(1-polycoeff(layer)).^polyexponent);
267+ end
268+
269 else
270 error('project3d error message: unknown projection type');
271 end
272Index: ../trunk-jpl/test/NightlyRun/test127.m
273===================================================================
274--- ../trunk-jpl/test/NightlyRun/test127.m (revision 26554)
275+++ ../trunk-jpl/test/NightlyRun/test127.m (revision 26555)
276@@ -17,7 +17,7 @@
277 massfluxatgate('name','MassFlux5','profilename',['../Exp/MassFlux5.exp'],'definitionstring','Outputdefinition5'),...
278 massfluxatgate('name','MassFlux6','profilename',['../Exp/MassFlux6.exp'],'definitionstring','Outputdefinition6')...
279 };
280-
281+md=SetMLHOBC(md);
282 md=solve(md,'Stressbalance');
283
284 %Fields and tolerances to track changes
285Index: ../trunk-jpl/test/NightlyRun/test128.m
286===================================================================
287--- ../trunk-jpl/test/NightlyRun/test128.m (revision 26554)
288+++ ../trunk-jpl/test/NightlyRun/test128.m (revision 26555)
289@@ -6,6 +6,7 @@
290 md.cluster=generic('name',oshostname(),'np',3);
291 md.transient.requested_outputs={'IceVolume','VxShear','VyShear','VxBase','VyBase','VxSurface','VySurface'};
292
293+md=SetMLHOBC(md);
294 md=solve(md,'Transient');
295
296 %Fields and tolerances to track changes
297Index: ../trunk-jpl/test/NightlyRun/test128.py
298===================================================================
299--- ../trunk-jpl/test/NightlyRun/test128.py (revision 26554)
300+++ ../trunk-jpl/test/NightlyRun/test128.py (revision 26555)
301@@ -15,6 +15,7 @@
302 md.cluster = generic('name', gethostname(), 'np', 3)
303 md.transient.requested_outputs = ['IceVolume','VxSurface','VySurface','VxShear','VyShear','VxBase','VyBase']
304
305+md = SetMLHOBC(md);
306 md = solve(md, 'Transient')
307
308 #Fields and tolerances to track changes
309Index: ../trunk-jpl/test/NightlyRun/test129.m
310===================================================================
311--- ../trunk-jpl/test/NightlyRun/test129.m (revision 26554)
312+++ ../trunk-jpl/test/NightlyRun/test129.m (revision 26555)
313@@ -13,6 +13,7 @@
314 md.timestepping.final_time=19;
315 md.settings.output_frequency=2;
316
317+md=SetMLHOBC(md);
318 md=solve(md,'Transient');
319 md2=solve(md,'Transient','restart',1);
320
321Index: ../trunk-jpl/test/NightlyRun/test129.py
322===================================================================
323--- ../trunk-jpl/test/NightlyRun/test129.py (revision 26554)
324+++ ../trunk-jpl/test/NightlyRun/test129.py (revision 26555)
325@@ -22,6 +22,7 @@
326 # time steps and resolution
327 md.timestepping.final_time = 19
328 md.settings.output_frequency = 2
329+md = SetMLHOBC(md);
330
331 md = solve(md, 'Transient')
332 md2 = copy.deepcopy(md)
333Index: ../trunk-jpl/test/NightlyRun/test248.m
334===================================================================
335--- ../trunk-jpl/test/NightlyRun/test248.m (revision 26554)
336+++ ../trunk-jpl/test/NightlyRun/test248.m (revision 26555)
337@@ -5,6 +5,7 @@
338 md=setflowequation(md,'MLHO','all');
339 md.cluster=generic('name',oshostname(),'np',3);
340 md.stressbalance.requested_outputs={'default','VxSurface','VySurface','VxShear','VyShear','VxBase','VyBase'};
341+md=SetMLHOBC(md);
342 md=solve(md,'Stressbalance');
343
344 %Fields and tolerances to track changes
345Index: ../trunk-jpl/test/NightlyRun/test248.py
346===================================================================
347--- ../trunk-jpl/test/NightlyRun/test248.py (revision 26554)
348+++ ../trunk-jpl/test/NightlyRun/test248.py (revision 26555)
349@@ -15,6 +15,7 @@
350 md = setflowequation(md, 'MLHO', 'all')
351 md.cluster = generic('name', gethostname(), 'np', 3)
352 md.stressbalance.requested_outputs = ['default', 'VxSurface', 'VySurface', 'VxShear', 'VyShear', 'VxBase', 'VyBase']
353+md = SetMLHOBC(md);
354 md = solve(md, 'Stressbalance')
355
356 #Fields and tolerances to track changes
357Index: ../trunk-jpl/test/NightlyRun/test249.m
358===================================================================
359--- ../trunk-jpl/test/NightlyRun/test249.m (revision 26554)
360+++ ../trunk-jpl/test/NightlyRun/test249.m (revision 26555)
361@@ -6,6 +6,7 @@
362 md.cluster=generic('name',oshostname(),'np',3);
363 md.transient.requested_outputs={'default','FloatingArea','GroundedArea','TotalGroundedBmb','TotalFloatingBmb'};
364 md.basalforcings.floatingice_melting_rate(:)=1;
365+md=SetMLHOBC(md);
366 md=solve(md,'Transient');
367
368 %Fields and tolerances to track changes
369Index: ../trunk-jpl/test/NightlyRun/test249.py
370===================================================================
371--- ../trunk-jpl/test/NightlyRun/test249.py (revision 26554)
372+++ ../trunk-jpl/test/NightlyRun/test249.py (revision 26555)
373@@ -16,6 +16,7 @@
374 md = setflowequation(md, 'MLHO', 'all')
375 md.cluster = generic('name', gethostname(), 'np', 3)
376 md.transient.requested_outputs = ['default', 'FloatingArea', 'GroundedArea', 'TotalFloatingBmb', 'TotalGroundedBmb']
377+md = SetMLHOBC(md);
378 md = solve(md, 'Transient')
379
380
381Index: ../trunk-jpl/test/NightlyRun/test254.m
382===================================================================
383--- ../trunk-jpl/test/NightlyRun/test254.m (revision 26554)
384+++ ../trunk-jpl/test/NightlyRun/test254.m (revision 26555)
385@@ -57,6 +57,7 @@
386 md.mask.ice_levelset=-1+nodeonicefront;
387
388 md.stressbalance.requested_outputs={'default','VySurface','VyShear','VyBase'};
389+md=SetMLHOBC(md);
390 md=solve(md,'Stressbalance');
391
392 %create analytical solution: strain rate is constant = ((rho_ice*g*h)/4B)^3 (Paterson, 4th Edition, page 292.
393Index: ../trunk-jpl/test/NightlyRun/test254.py
394===================================================================
395--- ../trunk-jpl/test/NightlyRun/test254.py (revision 26554)
396+++ ../trunk-jpl/test/NightlyRun/test254.py (revision 26555)
397@@ -69,6 +69,7 @@
398 md.mask.ice_levelset = -1 + nodeonicefront
399
400 md.stressbalance.requested_outputs = ['default', 'VySurface', 'VyShear', 'VyBase']
401+md = SetMLHOBC(md);
402 md = solve(md, 'Stressbalance')
403
404 # create analytical solution: strain rate is constant = ((rho_ice * g * h) / 4B)^3 (Paterson, 4th Edition, page 292.
405Index: ../trunk-jpl/test/NightlyRun/test255.m
406===================================================================
407--- ../trunk-jpl/test/NightlyRun/test255.m (revision 26554)
408+++ ../trunk-jpl/test/NightlyRun/test255.m (revision 26555)
409@@ -5,6 +5,7 @@
410 md=setflowequation(md,'MLHO','all');
411 md.cluster=generic('name',oshostname(),'np',3);
412 md.masstransport.hydrostatic_adjustment='Incremental';
413+md=SetMLHOBC(md);
414 md=solve(md,'Transient');
415
416 %Fields and tolerances to track changes
417Index: ../trunk-jpl/test/NightlyRun/test255.py
418===================================================================
419--- ../trunk-jpl/test/NightlyRun/test255.py (revision 26554)
420+++ ../trunk-jpl/test/NightlyRun/test255.py (revision 26555)
421@@ -14,6 +14,7 @@
422 md = setflowequation(md, 'MLHO', 'all')
423 md.cluster = generic('name', gethostname(), 'np', 3)
424 md.masstransport.hydrostatic_adjustment = 'Incremental'
425+md = SetMLHOBC(md);
426 md = solve(md, 'Transient')
427
428 #Fields and tolerances to track changes
429Index: ../trunk-jpl/test/NightlyRun/test256.m
430===================================================================
431--- ../trunk-jpl/test/NightlyRun/test256.m (revision 26554)
432+++ ../trunk-jpl/test/NightlyRun/test256.m (revision 26555)
433@@ -6,6 +6,7 @@
434 md.geometry.base=md.geometry.base+50.; md.geometry.surface=md.geometry.surface+50.;
435 md.cluster=generic('name',oshostname(),'np',1);
436 md.masstransport.hydrostatic_adjustment='Incremental';
437+md=SetMLHOBC(md);
438 md=solve(md,'Transient');
439
440 %Fields and tolerances to track changes
441Index: ../trunk-jpl/test/NightlyRun/test256.py
442===================================================================
443--- ../trunk-jpl/test/NightlyRun/test256.py (revision 26554)
444+++ ../trunk-jpl/test/NightlyRun/test256.py (revision 26555)
445@@ -16,6 +16,7 @@
446 md.geometry.surface = md.geometry.surface + 50.
447 md.cluster = generic('name', gethostname(), 'np', 1)
448 md.masstransport.hydrostatic_adjustment = 'Incremental'
449+md = SetMLHOBC(md);
450 md = solve(md, 'Transient')
451
452 #Fields and tolerances to track changes
453Index: ../trunk-jpl/test/NightlyRun/test330.m
454===================================================================
455--- ../trunk-jpl/test/NightlyRun/test330.m (revision 26554)
456+++ ../trunk-jpl/test/NightlyRun/test330.m (revision 26555)
457@@ -5,6 +5,7 @@
458 md=setflowequation(md,'MLHO','all');
459 md.cluster=generic('name',oshostname(),'np',3);
460 md.stressbalance.requested_outputs={'default','VxSurface','VySurface','VxShear','VyShear','VxBase','VyBase'};
461+md=SetMLHOBC(md);
462 md=solve(md,'Stressbalance');
463
464 %Fields and tolerances to track changes
465Index: ../trunk-jpl/test/NightlyRun/test330.py
466===================================================================
467--- ../trunk-jpl/test/NightlyRun/test330.py (revision 26554)
468+++ ../trunk-jpl/test/NightlyRun/test330.py (revision 26555)
469@@ -14,6 +14,7 @@
470 md = setflowequation(md, 'MLHO', 'all')
471 md.cluster = generic('name', gethostname(), 'np', 3)
472 md.stressbalance.requested_outputs = ['default', 'VxSurface', 'VySurface', 'VxShear', 'VyShear', 'VxBase', 'VyBase']
473+md = SetMLHOBC(md);
474 md = solve(md, 'Stressbalance')
475
476 #Fields and tolerances to track changes
477Index: ../trunk-jpl/test/NightlyRun/test332.m
478===================================================================
479--- ../trunk-jpl/test/NightlyRun/test332.m (revision 26554)
480+++ ../trunk-jpl/test/NightlyRun/test332.m (revision 26555)
481@@ -4,6 +4,7 @@
482 md=parameterize(md,'../Par/SquareSheetConstrained.par');
483 md=setflowequation(md,'MLHO','all');
484 md.cluster=generic('name',oshostname(),'np',3);
485+md=SetMLHOBC(md);
486 md=solve(md,'Transient');
487
488 %Fields and tolerances to track changes
489Index: ../trunk-jpl/test/NightlyRun/test332.py
490===================================================================
491--- ../trunk-jpl/test/NightlyRun/test332.py (revision 26554)
492+++ ../trunk-jpl/test/NightlyRun/test332.py (revision 26555)
493@@ -13,6 +13,7 @@
494 md = parameterize(md, '../Par/SquareSheetConstrained.py')
495 md = setflowequation(md, 'MLHO', 'all')
496 md.cluster = generic('name', gethostname(), 'np', 3)
497+md = SetMLHOBC(md);
498 md = solve(md, 'Transient')
499
500 #Fields and tolerances to track changes
501Index: ../trunk-jpl/test/NightlyRun/test333.m
502===================================================================
503--- ../trunk-jpl/test/NightlyRun/test333.m (revision 26554)
504+++ ../trunk-jpl/test/NightlyRun/test333.m (revision 26555)
505@@ -18,6 +18,7 @@
506 md.inversion.vx_obs=md.initialization.vx; md.inversion.vy_obs=md.initialization.vy;
507
508 md.cluster=generic('name',oshostname(),'np',3);
509+md=SetMLHOBC(md);
510 md=solve(md,'Stressbalance');
511
512 %Fields and tolerances to track changes
513Index: ../trunk-jpl/test/NightlyRun/test333.py
514===================================================================
515--- ../trunk-jpl/test/NightlyRun/test333.py (revision 26554)
516+++ ../trunk-jpl/test/NightlyRun/test333.py (revision 26555)
517@@ -31,6 +31,7 @@
518 md.inversion.vy_obs = md.initialization.vy
519
520 md.cluster = generic('name', gethostname(), 'np', 3)
521+md = SetMLHOBC(md);
522 md = solve(md, 'Stressbalance')
523
524 #Fields and tolerances to track changes
525Index: ../trunk-jpl/test/NightlyRun/test334.m
526===================================================================
527--- ../trunk-jpl/test/NightlyRun/test334.m (revision 26554)
528+++ ../trunk-jpl/test/NightlyRun/test334.m (revision 26555)
529@@ -19,6 +19,7 @@
530 md.verbose.control=true;
531
532 md.cluster=generic('name',oshostname(),'np',3);
533+md=SetMLHOBC(md);
534 md=solve(md,'Stressbalance');
535
536 %Fields and tolerances to track changes
537Index: ../trunk-jpl/test/NightlyRun/test334.py
538===================================================================
539--- ../trunk-jpl/test/NightlyRun/test334.py (revision 26554)
540+++ ../trunk-jpl/test/NightlyRun/test334.py (revision 26555)
541@@ -31,6 +31,7 @@
542
543
544 md.cluster = generic('name', gethostname(), 'np', 3)
545+md = SetMLHOBC(md);
546 md = solve(md, 'Stressbalance')
547
548
549Index: ../trunk-jpl/test/NightlyRun/test335.m
550===================================================================
551--- ../trunk-jpl/test/NightlyRun/test335.m (revision 26554)
552+++ ../trunk-jpl/test/NightlyRun/test335.m (revision 26555)
553@@ -6,6 +6,7 @@
554 md = extrude(md, 5, 1);
555 md.cluster=generic('name',oshostname(),'np',3);
556 md.stressbalance.requested_outputs={'default','VxSurface','VySurface','VxShear','VyShear','VxBase','VyBase'};
557+md=SetMLHOBC(md);
558 md=solve(md,'Stressbalance');
559
560 %Fields and tolerances to track changes
561Index: ../trunk-jpl/test/NightlyRun/test335.py
562===================================================================
563--- ../trunk-jpl/test/NightlyRun/test335.py (revision 26554)
564+++ ../trunk-jpl/test/NightlyRun/test335.py (revision 26555)
565@@ -12,6 +12,7 @@
566 md = setmask(md, '', '')
567 md = parameterize(md, '../Par/SquareSheetConstrained.py')
568 md = setflowequation(md, 'MLHO', 'all')
569+md = SetMLHOBC(md);
570 md.extrude(5, 1.)
571 md.cluster = generic('name', gethostname(), 'np', 3)
572 md.stressbalance.requested_outputs = ['default', 'VxSurface', 'VySurface', 'VxShear', 'VyShear', 'VxBase', 'VyBase']
573Index: ../trunk-jpl/test/NightlyRun/test446.m
574===================================================================
575--- ../trunk-jpl/test/NightlyRun/test446.m (revision 26554)
576+++ ../trunk-jpl/test/NightlyRun/test446.m (revision 26555)
577@@ -5,6 +5,7 @@
578 md=setflowequation(md,'MLHO','all');
579 md.cluster=generic('name',oshostname(),'np',3);
580 md.stressbalance.requested_outputs={'default','VxSurface','VySurface','VxShear','VyShear','VxBase','VyBase'};
581+md=SetMLHOBC(md);
582 md=solve(md,'Stressbalance');
583
584 %Fields and tolerances to track changes
585Index: ../trunk-jpl/test/NightlyRun/test446.py
586===================================================================
587--- ../trunk-jpl/test/NightlyRun/test446.py (revision 26554)
588+++ ../trunk-jpl/test/NightlyRun/test446.py (revision 26555)
589@@ -14,6 +14,7 @@
590 md = setflowequation(md, 'MLHO', 'all')
591 md.cluster = generic('name', gethostname(), 'np', 3)
592 md.stressbalance.requested_outputs = ['default', 'VxSurface', 'VySurface', 'VxShear', 'VyShear', 'VxBase', 'VyBase']
593+md = SetMLHOBC(md);
594 md = solve(md, 'Stressbalance')
595
596 #Fields and tolerances to track changes
597Index: ../trunk-jpl/test/NightlyRun/test447.m
598===================================================================
599--- ../trunk-jpl/test/NightlyRun/test447.m (revision 26554)
600+++ ../trunk-jpl/test/NightlyRun/test447.m (revision 26555)
601@@ -15,6 +15,7 @@
602 md=setflowequation(md,'MLHO','all');
603 md.cluster=generic('name',oshostname(),'np',3);
604 md.transient.requested_outputs={'default','GroundedArea','FloatingArea','TotalFloatingBmb','TotalGroundedBmb','TotalSmb'};
605+md=SetMLHOBC(md);
606 md=solve(md,'Transient');
607
608 %Fields and tolerances to track changes
609Index: ../trunk-jpl/test/NightlyRun/test447.py
610===================================================================
611--- ../trunk-jpl/test/NightlyRun/test447.py (revision 26554)
612+++ ../trunk-jpl/test/NightlyRun/test447.py (revision 26555)
613@@ -26,6 +26,7 @@
614 md = setflowequation(md, 'MLHO', 'all')
615 md.cluster = generic('name', gethostname(), 'np', 3)
616 md.transient.requested_outputs = ['default', 'GroundedArea', 'FloatingArea', 'TotalFloatingBmb', 'TotalGroundedBmb', 'TotalSmb']
617+md = SetMLHOBC(md);
618 md = solve(md, 'Transient')
619
620 #Fields and tolerances to track changes
621Index: ../trunk-jpl/test/NightlyRun/test448.m
622===================================================================
623--- ../trunk-jpl/test/NightlyRun/test448.m (revision 26554)
624+++ ../trunk-jpl/test/NightlyRun/test448.m (revision 26555)
625@@ -24,6 +24,7 @@
626 md.transient.isstressbalance=1;
627 md.transient.isgroundingline=1;
628
629+md=SetMLHOBC(md);
630 %test different grounding line dynamics.
631 md.groundingline.migration='AggressiveMigration';
632 md=solve(md,'Transient');
633Index: ../trunk-jpl/test/NightlyRun/test448.py
634===================================================================
635--- ../trunk-jpl/test/NightlyRun/test448.py (revision 26554)
636+++ ../trunk-jpl/test/NightlyRun/test448.py (revision 26555)
637@@ -37,6 +37,7 @@
638
639 #test different grounding line dynamics.
640 md.groundingline.migration = 'AggressiveMigration'
641+md = SetMLHOBC(md);
642 md = solve(md, 'Transient')
643 element_on_iceshelf_agressive = md.results.TransientSolution[0].MaskOceanLevelset
644 vel_agressive = md.results.TransientSolution[0].Vel
645Index: ../trunk-jpl/test/NightlyRun/test449.m
646===================================================================
647--- ../trunk-jpl/test/NightlyRun/test449.m (revision 26554)
648+++ ../trunk-jpl/test/NightlyRun/test449.m (revision 26555)
649@@ -44,6 +44,7 @@
650 md.timestepping.time_step=10;
651
652 md.cluster=generic('name',oshostname(),'np',3);
653+md=SetMLHOBC(md);
654 md=solve(md,'Transient');
655
656 %Fields and tolerances to track changes
657Index: ../trunk-jpl/test/NightlyRun/test449.py
658===================================================================
659--- ../trunk-jpl/test/NightlyRun/test449.py (revision 26554)
660+++ ../trunk-jpl/test/NightlyRun/test449.py (revision 26555)
661@@ -55,6 +55,7 @@
662 md.timestepping.time_step = 10
663
664 md.cluster = generic('name', gethostname(), 'np', 3)
665+md = SetMLHOBC(md);
666 md = solve(md, 'Transient')
667 #print md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate
668 #print md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate
669Index: ../trunk-jpl/test/NightlyRun/test518.m
670===================================================================
671--- ../trunk-jpl/test/NightlyRun/test518.m (revision 26554)
672+++ ../trunk-jpl/test/NightlyRun/test518.m (revision 26555)
673@@ -5,6 +5,7 @@
674 md=setflowequation(md,'MLHO','all');
675 md.cluster=generic('name',oshostname(),'np',3);
676 md.stressbalance.requested_outputs={'default','VxSurface','VySurface','VxShear','VyShear','VxBase','VyBase'};
677+md=SetMLHOBC(md);
678 md=solve(md,'Stressbalance');
679
680 %Fields and tolerances to track changes
681Index: ../trunk-jpl/test/NightlyRun/test518.py
682===================================================================
683--- ../trunk-jpl/test/NightlyRun/test518.py (revision 26554)
684+++ ../trunk-jpl/test/NightlyRun/test518.py (revision 26555)
685@@ -14,6 +14,7 @@
686 md = setflowequation(md, 'MLHO', 'all')
687 md.cluster = generic('name', gethostname(), 'np', 3)
688 md.stressbalance.requested_outputs = ['default', 'VxSurface', 'VySurface', 'VxShear', 'VyShear', 'VxBase', 'VyBase']
689+md = SetMLHOBC(md);
690 md = solve(md, 'Stressbalance')
691
692 #Fields and tolerances to track changes
693Index: ../trunk-jpl/test/NightlyRun/test519.m
694===================================================================
695--- ../trunk-jpl/test/NightlyRun/test519.m (revision 26554)
696+++ ../trunk-jpl/test/NightlyRun/test519.m (revision 26555)
697@@ -8,6 +8,7 @@
698 [md.mesh.lat,md.mesh.long] = xy2ll(md.mesh.x,md.mesh.y,-1);
699 md.mesh.scale_factor=0.9*ones(md.mesh.numberofvertices,1);
700 md.transient.requested_outputs={'default','IceVolume','IceVolumeScaled','GroundedArea','GroundedAreaScaled','FloatingArea','FloatingAreaScaled','TotalSmb','TotalSmbScaled','TotalFloatingBmb','TotalFloatingBmbScaled'};
701+md=SetMLHOBC(md);
702 md=solve(md,'Transient');
703
704 %Fields and tolerances to track changes
705Index: ../trunk-jpl/test/NightlyRun/test519.py
706===================================================================
707--- ../trunk-jpl/test/NightlyRun/test519.py (revision 26554)
708+++ ../trunk-jpl/test/NightlyRun/test519.py (revision 26555)
709@@ -15,6 +15,7 @@
710 md.mesh.scale_factor = 0.9 * np.ones((md.mesh.numberofvertices))
711 md.transient.requested_outputs = ['default', 'IceVolume', 'IceVolumeScaled', 'GroundedArea', 'GroundedAreaScaled', 'FloatingArea', 'FloatingAreaScaled', 'TotalSmb', 'TotalSmbScaled', 'TotalFloatingBmb', 'TotalFloatingBmbScaled']
712 md.cluster = generic('name', gethostname(), 'np', 3)
713+md = SetMLHOBC(md);
714 md = solve(md, 'Transient')
715
716 # Fields and tolerances to track changes
717Index: ../trunk-jpl/test/NightlyRun/test810.m
718===================================================================
719--- ../trunk-jpl/test/NightlyRun/test810.m (revision 26554)
720+++ ../trunk-jpl/test/NightlyRun/test810.m (revision 26555)
721@@ -14,6 +14,7 @@
722 md.transient.isthermal=0;
723 md.transient.isgroundingline=1;
724
725+md=SetMLHOBC(md);
726 md=solve(md,'Transient');
727
728 %Fields and tolerances to track changes
729Index: ../trunk-jpl/test/NightlyRun/test810.py
730===================================================================
731--- ../trunk-jpl/test/NightlyRun/test810.py (revision 26554)
732+++ ../trunk-jpl/test/NightlyRun/test810.py (revision 26555)
733@@ -21,6 +21,7 @@
734 md.transient.issmb = True
735 md.transient.isthermal = False
736 md.transient.isgroundingline = True
737+md = SetMLHOBC(md);
738
739 md = solve(md, 'Transient')
740
741Index: ../trunk-jpl/test/NightlyRun/test811.m
742===================================================================
743--- ../trunk-jpl/test/NightlyRun/test811.m (revision 26554)
744+++ ../trunk-jpl/test/NightlyRun/test811.m (revision 26555)
745@@ -17,6 +17,7 @@
746 md.frontalforcings.meltingrate=zeros(md.mesh.numberofvertices,1);
747 md.levelset.migration_max = 1e10;
748
749+md=SetMLHOBC(md);
750 md=solve(md,'Transient');
751
752 %Fields and tolerances to track changes
753Index: ../trunk-jpl/test/NightlyRun/test811.py
754===================================================================
755--- ../trunk-jpl/test/NightlyRun/test811.py (revision 26554)
756+++ ../trunk-jpl/test/NightlyRun/test811.py (revision 26555)
757@@ -27,6 +27,7 @@
758 md.frontalforcings.meltingrate = np.zeros((md.mesh.numberofvertices))
759 md.levelset.migration_max = 1e10
760
761+md = SetMLHOBC(md);
762 md = solve(md, 'Transient')
763
764 #Fields and tolerances to track changes
765Index: ../trunk-jpl/test/NightlyRun/test812.m
766===================================================================
767--- ../trunk-jpl/test/NightlyRun/test812.m (revision 26554)
768+++ ../trunk-jpl/test/NightlyRun/test812.m (revision 26555)
769@@ -34,6 +34,7 @@
770
771 md.transient.requested_outputs={'default','StrainRateparallel','StrainRateperpendicular','Calvingratex','Calvingratey','CalvingCalvingrate'};
772
773+md=SetMLHOBC(md);
774 md=solve(md,'Transient');
775
776 %Fields and tolerances to track changes
777Index: ../trunk-jpl/test/NightlyRun/test812.py
778===================================================================
779--- ../trunk-jpl/test/NightlyRun/test812.py (revision 26554)
780+++ ../trunk-jpl/test/NightlyRun/test812.py (revision 26555)
781@@ -43,6 +43,7 @@
782
783 md.transient.requested_outputs = ['default', 'StrainRateparallel', 'StrainRateperpendicular', 'Calvingratex', 'Calvingratey', 'CalvingCalvingrate']
784
785+md = SetMLHOBC(md);
786 md = solve(md, 'Transient')
787
788 #Fields and tolerances to track changes
789Index: ../trunk-jpl/src/m/boundaryconditions/SetMLHOBC.py
790===================================================================
791--- ../trunk-jpl/src/m/boundaryconditions/SetMLHOBC.py (nonexistent)
792+++ ../trunk-jpl/src/m/boundaryconditions/SetMLHOBC.py (revision 26555)
793@@ -0,0 +1,21 @@
794+import os
795+import MatlabFuncs as m
796+
797+
798+def SetMLHOBC(md):
799+ """
800+ SETMLHOBC - Create the boundary conditions for stressbalance for MLHO: VxBase, VyBase, VxShear, VyShear
801+
802+ Usage:
803+ md = SetIceShelfBC(md, varargin)
804+
805+ Example:
806+ md = SetIceShelfBC(md)
807+
808+ """
809+
810+ #node on Dirichlet (boundary and ~icefront)
811+ md.stressbalance.spcvx_base = md.stressbalance.spcvx
812+ md.stressbalance.spcvy_base = md.stressbalance.spcvy
813+ md.stressbalance.spcvx_shear = float('nan') * md.stressbalance.spcvx
814+ md.stressbalance.spcvy_shear = float('nan') * md.stressbalance.spcvy
815Index: ../trunk-jpl/src/m/extrusion/project3d.py
816===================================================================
817--- ../trunk-jpl/src/m/extrusion/project3d.py (revision 26554)
818+++ ../trunk-jpl/src/m/extrusion/project3d.py (revision 26555)
819@@ -14,7 +14,7 @@
820
821 arguments:
822 'vector': 2d vector
823- 'type': 'element' or 'node'
824+ 'type': 'element' or 'node' or 'poly'
825
826 options:
827 'layer' a layer number where vector should keep its values. If
828@@ -22,6 +22,7 @@
829 vector.
830 'padding': default to 0 (value adopted by other 3d layers not
831 being projected.
832+ 'degree': degree of polynomials when extrude from bottom to the top
833
834 Examples:
835 extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node', 'layer', 1, 'padding', NaN)
836@@ -41,6 +42,7 @@
837 vectype = options.getfieldvalue('type') #mandatory
838 layer = options.getfieldvalue('layer', 0) #optional (do all layers otherwise)
839 paddingvalue = options.getfieldvalue('padding', 0) #0 by default
840+ polyexponent = options.getfieldvalue('degree', 0) #0 by default, 0-degree polynomial
841
842 #Handle special case where vector2d is single element (differs from representation in MATLAB)
843 if isinstance(vector2d, (bool, int, float)):
844@@ -114,6 +116,37 @@
845 projected_vector[(i * md.mesh.numberofelements2d):((i + 1) * md.mesh.numberofelements2d), :] = vector2d
846 else:
847 projected_vector[((layer - 1) * md.mesh.numberofelements2d):(layer * md.mesh.numberofelements2d), :] = vector2d
848+ elif vectype.lower() == 'poly':
849+ #Initialize 3d vector
850+ if np.ndim(vector2d) == 1:
851+ if vector2d.shape[0] == md.mesh.numberofelements2d:
852+ projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements))).astype(vector2d.dtype)
853+ elif vector2d.shape[0] == md.mesh.numberofelements2d + 1:
854+ projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements + 1))).astype(vector2d.dtype)
855+ projected_vector[-1] = vector2d[-1]
856+ vector2d = vector2d[:-1]
857+ else:
858+ raise TypeError("vector length not supported")
859+ #Fill in
860+ if layer == 0:
861+ for i in range(md.mesh.numberoflayers - 1):
862+ projected_vector[(i * md.mesh.numberofelements2d):((i + 1) * md.mesh.numberofelements2d)] = vector2d*(1.0-(1.0-i/(md.mesh.numberoflayers - 1.0))**polyexponent)
863+ else:
864+ projected_vector[((layer - 1) * md.mesh.numberofelements2d):(layer * md.mesh.numberofelements2d)] = vector2d*(1.0-(1.0-layer/(md.mesh.numberoflayers - 1.0))**polyexponent)
865+ else:
866+ if vector2d.shape[0] == md.mesh.numberofelements2d:
867+ projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements, np.size(vector2d, axis=1)))).astype(vector2d.dtype)
868+ elif vector2d.shape[0] == md.mesh.numberofelements2d + 1:
869+ projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements + 1, np.size(vector2d, axis=1)))).astype(vector2d.dtype)
870+ projected_vector[-1, :] = vector2d[-1, :]
871+ vector2d = vector2d[:-1, :]
872+ else:
873+ raise TypeError("vector length not supported")
874+ #Fill in
875+ if layer == 0:
876+ for i in range(md.mesh.numberoflayers - 1):
877+ projected_vector[(i * md.mesh.numberofelements2d):((i + 1) * md.mesh.numberofelements2d), :] = vector2d*(1.0-(1.0-i/(md.mesh.numberoflayers - 1.0))**polyexponent)
878+ else:
879
880 else:
881 raise TypeError("project3d error message: unknown projection type")
882Index: ../trunk-jpl/test/NightlyRun/test127.py
883===================================================================
884--- ../trunk-jpl/test/NightlyRun/test127.py (revision 26554)
885+++ ../trunk-jpl/test/NightlyRun/test127.py (revision 26555)
886@@ -24,7 +24,7 @@
887 massfluxatgate('name', 'MassFlux4', 'profilename', '../Exp/MassFlux4.exp', 'definitionstring', 'Outputdefinition4'),
888 massfluxatgate('name', 'MassFlux5', 'profilename', '../Exp/MassFlux5.exp', 'definitionstring', 'Outputdefinition5'),
889 massfluxatgate('name', 'MassFlux6', 'profilename', '../Exp/MassFlux6.exp', 'definitionstring', 'Outputdefinition6')]
890-
891+md = SetMLHOBC(md);
892 md = solve(md, 'Stressbalance')
893
894 #Fields and tolerances to track changes
Note: See TracBrowser for help on using the repository browser.