Context Navigation

model.m@ 19527

Last change on this file since 19527 was 19527, checked in by Eric.Larour, 10 years ago

CHG: moved md.surfaceforcings to md.smb.
By doing so, had to rename the SMB class to SMBforcing class (it's just that, a mass_balance forcing inside
a SMB class, hence the name).
We also now have an smb_core solution, taken out of the mass transport core. Makes more sense long term.
Synced all enums according to the new changes, and operated the adjustments in all the test decks.

In addition, progressing in terms of GEMB integration into ISSM, specifically at the SMBgemb level (which
is spurring all the changes described above). Brought the class up to the level of the GEMB.m call in Alex's
code. Starting the C integration now.

File size: 61.7 KB

Line
1	%MODEL class definition
2	%
3	% Usage:
4	% md = model(varargin)
5
6	classdef model
7	properties (SetAccess=public) %Model fields
8	% {{{
9	%Careful here: no other class should be used as default value this is a bug of matlab
10	mesh = 0;
11	mask = 0;
12
13	geometry = 0;
14	constants = 0;
15	smb = 0;
16	basalforcings = 0;
17	materials = 0;
18	damage = 0;
19	friction = 0;
20	flowequation = 0;
21	timestepping = 0;
22	initialization = 0;
23	rifts = 0;
24
25	debug = 0;
26	verbose = 0;
27	settings = 0;
28	toolkits = 0;
29	cluster = 0;
30
31	balancethickness = 0;
32	stressbalance = 0;
33	groundingline = 0;
34	hydrology = 0;
35	masstransport = 0;
36	thermal = 0;
37	steadystate = 0;
38	transient = 0;
39	calving = 0;
40	gia = 0;
41
42	autodiff = 0;
43	flaim = 0;
44	inversion = 0;
45	qmu = 0;
46
47	results = 0;
48	outputdefinition = 0;
49	radaroverlay = 0;
50	miscellaneous = 0;
51	private = 0;
52
53	%}}}
54	end
55	methods (Static)
56	function md = loadobj(md) % {{{
57	% This function is directly called by matlab when a model object is
58	% loaded. If the input is a struct it is an old version of model and
59	% old fields must be recovered (make sure they are in the deprecated
60	% model properties)
61
62	if verLessThan('matlab','7.9'),
63	disp('Warning: your matlab version is old and there is a risk that load does not work correctly');
64	disp(' if the model is not loaded correctly, rename temporarily loadobj so that matlab does not use it');
65
66	% This is a Matlab bug: all the fields of md have their default value
67	% Example of error message:
68	% Warning: Error loading an object of class 'model':
69	% Undefined function or method 'exist' for input arguments of type 'cell'
70	%
71	% This has been fixed in MATLAB 7.9 (R2009b) and later versions
72	end
73
74	if isstruct(md)
75	disp('Recovering model object from a previous version');
76	md = structtomodel(model,md);
77	end
78
79	%2012 August 4th
80	if isa(md.materials,'materials'),
81	disp('Recovering old materials');
82	if numel(md.materials.rheology_Z)==1 & isnan(md.materials.rheology_Z),
83	md.materials=matice(md.materials);
84	else
85	md.materials=matdamageice(md.materials);
86	end
87	end
88	%2013 April 12
89	if numel(md.stressbalance.loadingforce==1)
90	md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
91	end
92	%2013 April 17
93	if isa(md.hydrology,'hydrology'),
94	disp('Recovering old hydrology class');
95	md.hydrology=hydrologyshreve(md.materials);
96	end
97	%2013 October 9
98	if ~isa(md.damage,'damage'),
99	md.damage=damage();
100	md.damage.D=zeros(md.mesh.numberofvertices,1);
101	md.damage.spcdamage=NaN*ones(md.mesh.numberofvertices,1);
102	end
103	%2013 November 18
104	if ~isa(md.outputdefinition,'outputdefinition'),
105	md.outputdefinition=outputdefinition();
106	end
107	%2014 January 9th
108	if isa(md.smb,'smb'),
109	disp('Recovering old smb class');
110	mass_balance=md.smb.mass_balance;
111	md.smb=SMBforcing();
112	md.smb.mass_balance=mass_balance;
113	end
114	%2014 March 26th
115	if isa(md.mesh,'mesh'),
116	disp('Recovering old mesh class');
117	if isprop(md.mesh,'dimension'),
118	if md.mesh.dimension==2,
119	md.mesh=mesh2d(md.mesh);
120	else
121	md.mesh=mesh3dprisms(md.mesh);
122	end
123	else
124	md.mesh=mesh2dvertical(md.mesh);
125	end
126	end
127	%2014 November 12
128	if isa(md.calving,'double');
129	md.calving=calving();
130	end
131	end% }}}
132	end
133	methods
134	function md = model(varargin) % {{{
135
136	switch nargin
137	case 0
138	md=setdefaultparameters(md);
139	case 1
140	error('model constructor not supported yet');
141
142	otherwise
143	error('model constructor error message: 0 of 1 argument only in input.');
144	end
145	end
146	%}}}
147	function md = checkmessage(md,string) % {{{
148	if(nargout~=1) error('wrong usage, model must be an output'); end
149	disp(['model not consistent: ' string]);
150	md.private.isconsistent=false;
151	end
152	%}}}
153	function md = collapse(md)% {{{
154	%COLLAPSE - collapses a 3d mesh into a 2d mesh
155	%
156	% This routine collapses a 3d model into a 2d model
157	% and collapses all the fileds of the 3d model by
158	% taking their depth-averaged values
159	%
160	% Usage:
161	% md=collapse(md)
162	%
163	% See also: EXTRUDE, MODELEXTRACT
164
165	%Check that the model is really a 3d model
166	if ~strcmp(md.mesh.elementtype(),'Penta'),
167	error('collapse error message: only 3d mesh can be collapsed')
168	end
169
170	%Start with changing all the fields from the 3d mesh
171
172	%dealing with the friciton law
173	%drag is limited to nodes that are on the bedrock.
174	if isa(md.friction,'friction'),
175	md.friction.coefficient=project2d(md,md.friction.coefficient,1);
176	md.friction.p=project2d(md,md.friction.p,1);
177	md.friction.q=project2d(md,md.friction.q,1);
178	elseif isa(md.friction,'frictionhydro'),
179	md.friction.q=project2d(md,md.friction.q,1);
180	md.friction.C=project2d(md,md.friction.C,1);
181	md.friction.As=project2d(md,md.friction.As,1);
182	md.friction.effective_pressure=project2d(md,md.friction.effective_pressure,1);
183	elseif isa(md.friction,'frictionwaterlayer'),
184	md.friction.coefficient=project2d(md,md.friction.coefficient,1);
185	md.friction.p=project2d(md,md.friction.p,1);
186	md.friction.q=project2d(md,md.friction.q,1);
187	md.friction.water_layer=project2d(md,md.friction.water_layer,1);
188	elseif isa(md.friction,'frictionweertman'),
189	md.friction.C=project2d(md,md.friction.C,1);
190	md.friction.m=project2d(md,md.friction.m,1);
191	end
192
193	%observations
194	if ~isnan(md.inversion.vx_obs), md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); end;
195	if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); end;
196	if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end;
197	if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end;
198	if numel(md.inversion.min_parameters)>1, md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end;
199	if numel(md.inversion.max_parameters)>1, md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end;
200	if isa(md.smb,'SMBforcing') & ~isnan(md.smb.mass_balance),
201	md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers);
202	elseif isa(md.smb,'SMBhenning') & ~isnan(md.smb.smbref),
203	md.smb.smbref=project2d(md,md.smb.smbref,md.mesh.numberoflayers);
204	end;
205	if ~isnan(md.balancethickness.thickening_rate), md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.mesh.numberoflayers); end;
206
207	%results
208	if ~isnan(md.initialization.vx),md.initialization.vx=DepthAverage(md,md.initialization.vx);end;
209	if ~isnan(md.initialization.vy),md.initialization.vy=DepthAverage(md,md.initialization.vy);end;
210	if ~isnan(md.initialization.vz),md.initialization.vz=DepthAverage(md,md.initialization.vz);end;
211	if ~isnan(md.initialization.vel),md.initialization.vel=DepthAverage(md,md.initialization.vel);end;
212	if ~isnan(md.initialization.temperature),md.initialization.temperature=DepthAverage(md,md.initialization.temperature);end;
213	if ~isnan(md.initialization.pressure),md.initialization.pressure=project2d(md,md.initialization.pressure,1);end;
214	if ~isnan(md.initialization.sediment_head),md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1);end;
215	if ~isnan(md.initialization.epl_head),md.initialization.epl_head=project2d(md,md.initialization.epl_head,1);end;
216	if ~isnan(md.initialization.epl_thickness),md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1);end;
217	%gia
218	if ~isnan(md.gia.mantle_viscosity), md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1); end
219	if ~isnan(md.gia.lithosphere_thickness), md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1); end
220
221	%elementstype
222	if ~isnan(md.flowequation.element_equation)
223	md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1);
224	md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1);
225	md.flowequation.borderSSA=project2d(md,md.flowequation.borderSSA,1);
226	md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1);
227	md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1);
228	end
229
230	%boundary conditions
231	md.stressbalance.spcvx=project2d(md,md.stressbalance.spcvx,md.mesh.numberoflayers);
232	md.stressbalance.spcvy=project2d(md,md.stressbalance.spcvy,md.mesh.numberoflayers);
233	md.stressbalance.spcvz=project2d(md,md.stressbalance.spcvz,md.mesh.numberoflayers);
234	md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers);
235	md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers);
236	md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers);
237	if ~isnan(md.damage.spcdamage), md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers); end
238	md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);
239
240	% Hydrologydc variables
241	if isa(md.hydrology,'hydrologydc');
242	md.hydrology.spcsediment_head=project2d(md,md.hydrology.spcsediment_head,1);
243	md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1);
244	md.hydrology.sediment_transmitivity=project2d(md,md.hydrology.sediment_transmitivity,1);
245	md.hydrology.basal_moulin_input=project2d(md,md.hydrology.basal_moulin_input,1);
246	if(md.hydrology.isefficientlayer==1)
247	md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1);
248	end
249	end
250
251	%materials
252	md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B);
253	md.materials.rheology_n=project2d(md,md.materials.rheology_n,1);
254
255	%damage:
256	if md.damage.isdamage,
257	md.damage.D=DepthAverage(md,md.damage.D);
258	end
259
260	%special for thermal modeling:
261	if ~isnan(md.basalforcings.groundedice_melting_rate),
262	md.basalforcings.groundedice_melting_rate=project2d(md,md.basalforcings.groundedice_melting_rate,1);
263	end
264	if ~isnan(md.basalforcings.floatingice_melting_rate),
265	md.basalforcings.floatingice_melting_rate=project2d(md,md.basalforcings.floatingice_melting_rate,1);
266	end
267	md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); %bedrock only gets geothermal flux
268
269	%update of connectivity matrix
270	md.mesh.average_vertex_connectivity=25;
271
272	%Collapse the mesh
273	nodes2d=md.mesh.numberofvertices2d;
274	elements2d=md.mesh.numberofelements2d;
275
276	%parameters
277	md.geometry.surface=project2d(md,md.geometry.surface,1);
278	md.geometry.thickness=project2d(md,md.geometry.thickness,1);
279	md.geometry.base=project2d(md,md.geometry.base,1);
280	if ~isnan(md.geometry.bed),
281	md.geometry.bed=project2d(md,md.geometry.bed,1);
282	end
283
284	if ~isnan(md.mask.groundedice_levelset),
285	md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1);
286	end
287	if ~isnan(md.mask.ice_levelset),
288	md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1);
289	end
290
291	%lat long
292	if numel(md.mesh.lat) ==md.mesh.numberofvertices, md.mesh.lat=project2d(md,md.mesh.lat,1); end
293	if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end
294
295	%Initialize with the 2d mesh
296	mesh=mesh2d();
297	mesh.x=md.mesh.x2d;
298	mesh.y=md.mesh.y2d;
299	mesh.numberofvertices=md.mesh.numberofvertices2d;
300	mesh.numberofelements=md.mesh.numberofelements2d;
301	mesh.elements=md.mesh.elements2d;
302	if ~isnan(md.mesh.vertexonboundary), mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1); end
303	if ~isnan(md.mesh.elementconnectivity), mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1); end
304	md.mesh=mesh;
305	md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
306	md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
307	md.mesh.segments=contourenvelope(md);
308
309	end % }}}
310	function md2 = extract(md,area) % {{{
311	%extract - extract a model according to an Argus contour or flag list
312	%
313	% This routine extracts a submodel from a bigger model with respect to a given contour
314	% md must be followed by the corresponding exp file or flags list
315	% It can either be a domain file (argus type, .exp extension), or an array of element flags.
316	% If user wants every element outside the domain to be
317	% extract2d, add '~' to the name of the domain file (ex: '~HO.exp');
318	% an empty string '' will be considered as an empty domain
319	% a string 'all' will be considered as the entire domain
320	%
321	% Usage:
322	% md2=extract(md,area);
323	%
324	% Examples:
325	% md2=extract(md,'Domain.exp');
326	%
327	% See also: EXTRUDE, COLLAPSE
328
329	%copy model
330	md1=md;
331
332	%some checks
333	if ((nargin~=2) \| (nargout~=1)),
334	help extract
335	error('extract error message: bad usage');
336	end
337
338	%get elements that are inside area
339	flag_elem=FlagElements(md1,area);
340	if ~any(flag_elem),
341	error('extracted model is empty');
342	end
343
344	%kick out all elements with 3 dirichlets
345	spc_elem=find(~flag_elem);
346	spc_node=sort(unique(md1.mesh.elements(spc_elem,:)));
347	flag=ones(md1.mesh.numberofvertices,1);
348	flag(spc_node)=0;
349	pos=find(sum(flag(md1.mesh.elements),2)==0);
350	flag_elem(pos)=0;
351
352	%extracted elements and nodes lists
353	pos_elem=find(flag_elem);
354	pos_node=sort(unique(md1.mesh.elements(pos_elem,:)));
355
356	%keep track of some fields
357	numberofvertices1=md1.mesh.numberofvertices;
358	numberofelements1=md1.mesh.numberofelements;
359	numberofvertices2=length(pos_node);
360	numberofelements2=length(pos_elem);
361	flag_node=zeros(numberofvertices1,1);
362	flag_node(pos_node)=1;
363
364	%Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
365	Pelem=zeros(numberofelements1,1);
366	Pelem(pos_elem)=[1:numberofelements2]';
367	Pnode=zeros(numberofvertices1,1);
368	Pnode(pos_node)=[1:numberofvertices2]';
369
370	%renumber the elements (some nodes won't exist anymore)
371	elements_1=md1.mesh.elements;
372	elements_2=elements_1(pos_elem,:);
373	elements_2(:,1)=Pnode(elements_2(:,1));
374	elements_2(:,2)=Pnode(elements_2(:,2));
375	elements_2(:,3)=Pnode(elements_2(:,3));
376	if isa(md1.mesh,'mesh3dprisms'),
377	elements_2(:,4)=Pnode(elements_2(:,4));
378	elements_2(:,5)=Pnode(elements_2(:,5));
379	elements_2(:,6)=Pnode(elements_2(:,6));
380	end
381
382	%OK, now create the new model!
383
384	%take every field from model
385	md2=md1;
386
387	%automatically modify fields
388
389	%loop over model fields
390	model_fields=fields(md1);
391	for i=1:length(model_fields),
392	%get field
393	field=md1.(model_fields{i});
394	fieldsize=size(field);
395	if isobject(field), %recursive call
396	object_fields=fields(md1.(model_fields{i}));
397	for j=1:length(object_fields),
398	%get field
399	field=md1.(model_fields{i}).(object_fields{j});
400	fieldsize=size(field);
401	%size = number of nodes * n
402	if fieldsize(1)==numberofvertices1
403	md2.(model_fields{i}).(object_fields{j})=field(pos_node,:);
404	elseif (fieldsize(1)==numberofvertices1+1)
405	md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)];
406	%size = number of elements * n
407	elseif fieldsize(1)==numberofelements1
408	md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:);
409	end
410	end
411	else
412	%size = number of nodes * n
413	if fieldsize(1)==numberofvertices1
414	md2.(model_fields{i})=field(pos_node,:);
415	elseif (fieldsize(1)==numberofvertices1+1)
416	md2.(model_fields{i})=[field(pos_node,:); field(end,:)];
417	%size = number of elements * n
418	elseif fieldsize(1)==numberofelements1
419	md2.(model_fields{i})=field(pos_elem,:);
420	end
421	end
422	end
423
424	%modify some specific fields
425
426	%Mesh
427	md2.mesh.numberofelements=numberofelements2;
428	md2.mesh.numberofvertices=numberofvertices2;
429	md2.mesh.elements=elements_2;
430
431	%mesh.uppervertex mesh.lowervertex
432	if isa(md1.mesh,'mesh3dprisms'),
433	md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node);
434	pos=find(~isnan(md2.mesh.uppervertex));
435	md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos));
436
437	md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node);
438	pos=find(~isnan(md2.mesh.lowervertex));
439	md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos));
440
441	md2.mesh.upperelements=md1.mesh.upperelements(pos_elem);
442	pos=find(~isnan(md2.mesh.upperelements));
443	md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos));
444
445	md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem);
446	pos=find(~isnan(md2.mesh.lowerelements));
447	md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos));
448	end
449
450	%Initial 2d mesh
451	if isa(md1.mesh,'mesh3dprisms'),
452	flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d);
453	pos_elem_2d=find(flag_elem_2d);
454	flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d);
455	pos_node_2d=find(flag_node_2d);
456
457	md2.mesh.numberofelements2d=length(pos_elem_2d);
458	md2.mesh.numberofvertices2d=length(pos_node_2d);
459	md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:);
460	md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1));
461	md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2));
462	md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3));
463
464	md2.mesh.x2d=md1.mesh.x(pos_node_2d);
465	md2.mesh.y2d=md1.mesh.y(pos_node_2d);
466	end
467
468	%Edges
469	if(dimension(md.mesh)==2),
470	if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs...
471	%renumber first two columns
472	pos=find(md2.mesh.edges(:,4)~=-1);
473	md2.mesh.edges(: ,1)=Pnode(md2.mesh.edges(:,1));
474	md2.mesh.edges(: ,2)=Pnode(md2.mesh.edges(:,2));
475	md2.mesh.edges(: ,3)=Pelem(md2.mesh.edges(:,3));
476	md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4));
477	%remove edges when the 2 vertices are not in the domain.
478	md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:);
479	%Replace all zeros by -1 in the last two columns
480	pos=find(md2.mesh.edges(:,3)==0);
481	md2.mesh.edges(pos,3)=-1;
482	pos=find(md2.mesh.edges(:,4)==0);
483	md2.mesh.edges(pos,4)=-1;
484	%Invert -1 on the third column with last column (Also invert first two columns!!)
485	pos=find(md2.mesh.edges(:,3)==-1);
486	md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4);
487	md2.mesh.edges(pos,4)=-1;
488	values=md2.mesh.edges(pos,2);
489	md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1);
490	md2.mesh.edges(pos,1)=values;
491	%Finally remove edges that do not belong to any element
492	pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1);
493	md2.mesh.edges(pos,:)=[];
494	end
495	end
496
497	%Penalties
498	if ~isnan(md2.stressbalance.vertex_pairing),
499	for i=1:size(md1.stressbalance.vertex_pairing,1);
500	md2.stressbalance.vertex_pairing(i,:)=Pnode(md1.stressbalance.vertex_pairing(i,:));
501	end
502	md2.stressbalance.vertex_pairing=md2.stressbalance.vertex_pairing(find(md2.stressbalance.vertex_pairing(:,1)),:);
503	end
504	if ~isnan(md2.masstransport.vertex_pairing),
505	for i=1:size(md1.masstransport.vertex_pairing,1);
506	md2.masstransport.vertex_pairing(i,:)=Pnode(md1.masstransport.vertex_pairing(i,:));
507	end
508	md2.masstransport.vertex_pairing=md2.masstransport.vertex_pairing(find(md2.masstransport.vertex_pairing(:,1)),:);
509	end
510
511	%recreate segments
512	if isa(md1.mesh,'mesh2d'),
513	md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices);
514	md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
515	md2.mesh.segments=contourenvelope(md2);
516	md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
517	else
518	%First do the connectivity for the contourenvelope in 2d
519	md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d);
520	md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity);
521	segments=contourenvelope(md2);
522	md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(segments(:,1:2))=1;
523	md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1);
524	%Then do it for 3d as usual
525	md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices);
526	md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
527	end
528
529	%Boundary conditions: Dirichlets on new boundary
530	%Catch the elements that have not been extracted
531	orphans_elem=find(~flag_elem);
532	orphans_node=unique(md1.mesh.elements(orphans_elem,:))';
533	%Figure out which node are on the boundary between md2 and md1
534	nodestoflag1=intersect(orphans_node,pos_node);
535	nodestoflag2=Pnode(nodestoflag1);
536	if numel(md1.stressbalance.spcvx)>1 & numel(md1.stressbalance.spcvy)>2 & numel(md1.stressbalance.spcvz)>2,
537	if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1
538	md2.stressbalance.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2);
539	md2.stressbalance.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2);
540	else
541	md2.stressbalance.spcvx(nodestoflag2)=NaN;
542	md2.stressbalance.spcvy(nodestoflag2)=NaN;
543	disp(' ')
544	disp('!! extract warning: spc values should be checked !!')
545	disp(' ')
546	end
547	%put 0 for vz
548	md2.stressbalance.spcvz(nodestoflag2)=0;
549	end
550	if ~isnan(md1.thermal.spctemperature),
551	md2.thermal.spctemperature(nodestoflag2,1)=1;
552	end
553
554	%Results fields
555	if isstruct(md1.results),
556	md2.results=struct();
557	solutionfields=fields(md1.results);
558	for i=1:length(solutionfields),
559	if isstruct(md1.results.(solutionfields{i}))
560	%get subfields
561	solutionsubfields=fields(md1.results.(solutionfields{i}));
562	for j=1:length(solutionsubfields),
563	field=md1.results.(solutionfields{i}).(solutionsubfields{j});
564	if length(field)==numberofvertices1,
565	md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_node);
566	elseif length(field)==numberofelements1,
567	md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_elem);
568	else
569	md2.results.(solutionfields{i}).(solutionsubfields{j})=field;
570	end
571	end
572	else
573	field=md1.results.(solutionfields{i});
574	if length(field)==numberofvertices1,
575	md2.results.(solutionfields{i})=field(pos_node);
576	elseif length(field)==numberofelements1,
577	md2.results.(solutionfields{i})=field(pos_elem);
578	else
579	md2.results.(solutionfields{i})=field;
580	end
581	end
582	end
583	end
584
585	%Keep track of pos_node and pos_elem
586	md2.mesh.extractedvertices=pos_node;
587	md2.mesh.extractedelements=pos_elem;
588	end % }}}
589	function md = extrude(md,varargin) % {{{
590	%EXTRUDE - vertically extrude a 2d mesh
591	%
592	% vertically extrude a 2d mesh and create corresponding 3d mesh.
593	% The vertical distribution can:
594	% - follow a polynomial law
595	% - follow two polynomial laws, one for the lower part and one for the upper part of the mesh
596	% - be discribed by a list of coefficients (between 0 and 1)
597	%
598	%
599	% Usage:
600	% md=extrude(md,numlayers,extrusionexponent);
601	% md=extrude(md,numlayers,lowerexponent,upperexponent);
602	% md=extrude(md,listofcoefficients);
603	%
604	% Example:
605	% md=extrude(md,15,1.3);
606	% md=extrude(md,15,1.3,1.2);
607	% md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]);
608	%
609	% See also: MODELEXTRACT, COLLAPSE
610
611	%some checks on list of arguments
612	if ((nargin>4) \| (nargin<2) \| (nargout~=1)),
613	help extrude;
614	error('extrude error message');
615	end
616
617	%Extrude the mesh
618	if nargin==2, %list of coefficients
619	clist=varargin{1};
620	if any(clist<0) \| any(clist>1),
621	error('extrusioncoefficients must be between 0 and 1');
622	end
623	extrusionlist=sort(unique([clist(:);0;1]));
624	numlayers=length(extrusionlist);
625	elseif nargin==3, %one polynomial law
626	if varargin{2}<=0,
627	help extrude;
628	error('extrusionexponent must be >=0');
629	end
630	numlayers=varargin{1};
631	extrusionlist=((0:1:numlayers-1)/(numlayers-1)).^varargin{2};
632	elseif nargin==4, %two polynomial laws
633	numlayers=varargin{1};
634	lowerexp=varargin{2};
635	upperexp=varargin{3};
636
637	if varargin{2}<=0 \| varargin{3}<=0,
638	help extrude;
639	error('lower and upper extrusionexponents must be >=0');
640	end
641
642	lowerextrusionlist=[(0:2/(numlayers-1):1).^lowerexp]/2;
643	upperextrusionlist=[(0:2/(numlayers-1):1).^upperexp]/2;
644	extrusionlist=sort(unique([lowerextrusionlist 1-upperextrusionlist]));
645
646	end
647
648	if numlayers<2,
649	error('number of layers should be at least 2');
650	end
651	if strcmp(md.mesh.domaintype(),'3D')
652	error('Cannot extrude a 3d mesh (extrude cannot be called more than once)');
653	end
654
655	%Initialize with the 2d mesh
656	mesh2d = md.mesh;
657	md.mesh=mesh3dprisms();
658	md.mesh.x = mesh2d.x;
659	md.mesh.y = mesh2d.y;
660	md.mesh.elements = mesh2d.elements;
661	md.mesh.numberofelements = mesh2d.numberofelements;
662	md.mesh.numberofvertices = mesh2d.numberofvertices;
663
664	md.mesh.lat = mesh2d.lat;
665	md.mesh.long = mesh2d.long;
666	md.mesh.epsg = mesh2d.epsg;
667
668	md.mesh.vertexonboundary = mesh2d.vertexonboundary;
669	md.mesh.vertexconnectivity = mesh2d.vertexconnectivity;
670	md.mesh.elementconnectivity = mesh2d.elementconnectivity;
671	md.mesh.average_vertex_connectivity = mesh2d.average_vertex_connectivity;
672
673	md.mesh.extractedvertices = mesh2d.extractedvertices;
674	md.mesh.extractedelements = mesh2d.extractedelements;
675
676	x3d=[];
677	y3d=[];
678	z3d=[]; %the lower node is on the bed
679	thickness3d=md.geometry.thickness; %thickness and bed for these nodes
680	bed3d=md.geometry.base;
681
682	%Create the new layers
683	for i=1:numlayers,
684	x3d=[x3d; md.mesh.x];
685	y3d=[y3d; md.mesh.y];
686	%nodes are distributed between bed and surface accordingly to the given exponent
687	z3d=[z3d; bed3d+thickness3d*extrusionlist(i)];
688	end
689	number_nodes3d=size(x3d,1); %number of 3d nodes for the non extruded part of the mesh
690
691	%Extrude elements
692	elements3d=[];
693	for i=1:numlayers-1,
694	elements3d=[elements3d;[md.mesh.elements+(i-1)md.mesh.numberofvertices md.mesh.elements+imd.mesh.numberofvertices]]; %Create the elements of the 3d mesh for the non extruded part
695	end
696	number_el3d=size(elements3d,1); %number of 3d nodes for the non extruded part of the mesh
697
698	%Keep a trace of lower and upper nodes
699	lowervertex=NaN*ones(number_nodes3d,1);
700	uppervertex=NaN*ones(number_nodes3d,1);
701	lowervertex(md.mesh.numberofvertices+1:end)=1:(numlayers-1)*md.mesh.numberofvertices;
702	uppervertex(1:(numlayers-1)*md.mesh.numberofvertices)=md.mesh.numberofvertices+1:number_nodes3d;
703	md.mesh.lowervertex=lowervertex;
704	md.mesh.uppervertex=uppervertex;
705
706	%same for lower and upper elements
707	lowerelements=NaN*ones(number_el3d,1);
708	upperelements=NaN*ones(number_el3d,1);
709	lowerelements(md.mesh.numberofelements+1:end)=1:(numlayers-2)*md.mesh.numberofelements;
710	upperelements(1:(numlayers-2)md.mesh.numberofelements)=md.mesh.numberofelements+1:(numlayers-1)md.mesh.numberofelements;
711	md.mesh.lowerelements=lowerelements;
712	md.mesh.upperelements=upperelements;
713
714	%Save old mesh
715	md.mesh.x2d=md.mesh.x;
716	md.mesh.y2d=md.mesh.y;
717	md.mesh.elements2d=md.mesh.elements;
718	md.mesh.numberofelements2d=md.mesh.numberofelements;
719	md.mesh.numberofvertices2d=md.mesh.numberofvertices;
720
721	%Build global 3d mesh
722	md.mesh.elements=elements3d;
723	md.mesh.x=x3d;
724	md.mesh.y=y3d;
725	md.mesh.z=z3d;
726	md.mesh.numberofelements=number_el3d;
727	md.mesh.numberofvertices=number_nodes3d;
728	md.mesh.numberoflayers=numlayers;
729
730	%Ok, now deal with the other fields from the 2d mesh:
731
732	%bedinfo and surface info
733	md.mesh.vertexonbase=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1);
734	md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers);
735	md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node');
736
737	%lat long
738	md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node');
739	md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node');
740
741	md.geometry=extrude(md.geometry,md);
742	md.friction = extrude(md.friction,md);
743	md.inversion = extrude(md.inversion,md);
744	md.smb = extrude(md.smb,md);
745	md.initialization = extrude(md.initialization,md);
746
747	md.flowequation=md.flowequation.extrude(md);
748	md.stressbalance=extrude(md.stressbalance,md);
749	md.thermal=md.thermal.extrude(md);
750	md.masstransport=md.masstransport.extrude(md);
751	md.calving=extrude(md.calving,md);
752	md.hydrology = extrude(md.hydrology,md);
753
754	%connectivity
755	if ~isnan(md.mesh.elementconnectivity)
756	md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1);
757	md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN;
758	for i=2:numlayers-1,
759	md.mesh.elementconnectivity((i-1)md.mesh.numberofelements2d+1:(i)md.mesh.numberofelements2d,:)...
760	=md.mesh.elementconnectivity((i-1)md.mesh.numberofelements2d+1:(i)md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d;
761	end
762	md.mesh.elementconnectivity(find(isnan(md.mesh.elementconnectivity)))=0;
763	end
764
765	md.materials=extrude(md.materials,md);
766	md.damage=extrude(md.damage,md);
767	md.mask=extrude(md.mask,md);
768	md.qmu=extrude(md.qmu,md);
769	md.basalforcings=extrude(md.basalforcings,md);
770
771	%increase connectivity if less than 25:
772	if md.mesh.average_vertex_connectivity<=25,
773	md.mesh.average_vertex_connectivity=100;
774	end
775	end % }}}
776	function md = structtomodel(md,structmd) % {{{
777
778	if ~isstruct(structmd) error('input model is not a structure'); end
779
780	%loaded model is a struct, initialize output and recover all fields
781	md = structtoobj(model,structmd);
782
783	%Old field now classes
784	if (isfield(structmd,'timestepping') & isnumeric(md.timestepping)), md.timestepping=timestepping(); end
785	if (isfield(structmd,'mask') & isnumeric(md.mask)),md.mask=mask(); end
786
787	%Field name change
788	if isfield(structmd,'drag'), md.friction.coefficient=structmd.drag; end
789	if isfield(structmd,'p'), md.friction.p=structmd.p; end
790	if isfield(structmd,'q'), md.friction.q=structmd.p; end
791	if isfield(structmd,'melting'), md.basalforcings.floatingice_melting_rate=structmd.melting; end
792	if isfield(structmd,'melting_rate'), md.basalforcings.floatingice_melting_rate=structmd.melting_rate; end
793	if isfield(structmd,'melting_rate'), md.basalforcings.groundedice_melting_rate=structmd.melting_rate; end
794	if isfield(structmd,'accumulation'), md.smb.mass_balance=structmd.accumulation; end
795	if isfield(structmd,'numberofgrids'), md.mesh.numberofvertices=structmd.numberofgrids; end
796	if isfield(structmd,'numberofgrids2d'), md.mesh.numberofvertices2d=structmd.numberofgrids2d; end
797	if isfield(structmd,'uppergrids'), md.mesh.uppervertex=structmd.uppergrids; end
798	if isfield(structmd,'lowergrids'), md.mesh.lowervertex=structmd.lowergrids; end
799	if isfield(structmd,'gridonbase'), md.mesh.vertexonbase=structmd.gridonbase; end
800	if isfield(structmd,'gridonsurface'), md.mesh.vertexonsurface=structmd.gridonsurface; end
801	if isfield(structmd,'extractedgrids'), md.mesh.extractedvertices=structmd.extractedgrids; end
802	if isfield(structmd,'gridonboundary'), md.mesh.vertexonboundary=structmd.gridonboundary; end
803	if isfield(structmd,'petscoptions') & ~isempty(structmd.petscoptions), md.toolkits=structmd.petscoptions; end
804	if isfield(structmd,'g'), md.constants.g=structmd.g; end
805	if isfield(structmd,'yts'), md.constants.yts=structmd.yts; end
806	if isfield(structmd,'surface_mass_balance'), md.smb.mass_balance=structmd.surface_mass_balance; end
807	if isfield(structmd,'basal_melting_rate'), md.basalforcings.floatingice_melting_rate=structmd.basal_melting_rate; end
808	if isfield(structmd,'geothermalflux'), md.basalforcings.geothermalflux=structmd.geothermalflux; end
809	if isfield(structmd,'drag'), md.friction.coefficient=structmd.drag; end
810	if isfield(structmd,'drag_coefficient'), md.friction.coefficient=structmd.drag_coefficient; end
811	if isfield(structmd,'drag_p'), md.friction.p=structmd.drag_p; end
812	if isfield(structmd,'drag_q'), md.friction.q=structmd.drag_q; end
813	if isfield(structmd,'riftproperties'), %old implementation
814	md.rifts=rifts();
815	md.rifts.riftproperties=structmd.riftproperties;
816	md.rifts.riftstruct=structmd.rifts;
817	md.rifts.riftproperties=structmd.riftinfo;
818	end
819	if isfield(structmd,'bamg'), md.private.bamg=structmd.bamg; end
820	if isfield(structmd,'lowmem'), md.settings.lowmem=structmd.lowmem; end
821	if isfield(structmd,'io_gather'), md.settings.io_gather=structmd.io_gather; end
822	if isfield(structmd,'spcwatercolumn'), md.hydrology.spcwatercolumn=structmd.spcwatercolumn; end
823	if isfield(structmd,'hydro_n'), md.hydrology.n=structmd.hydro_n; end
824	if isfield(structmd,'hydro_p'), md.hydrology.p=structmd.hydro_p; end
825	if isfield(structmd,'hydro_q'), md.hydrology.q=structmd.hydro_q; end
826	if isfield(structmd,'hydro_CR'), md.hydrology.CR=structmd.hydro_CR; end
827	if isfield(structmd,'hydro_kn'), md.hydrology.kn=structmd.hydro_kn; end
828	if isfield(structmd,'spctemperature'), md.thermal.spctemperature=structmd.spctemperature; end
829	if isfield(structmd,'min_thermal_constraints'), md.thermal.penalty_threshold=structmd.min_thermal_constraints; end
830	if isfield(structmd,'artificial_diffusivity'), md.thermal.stabilization=structmd.artificial_diffusivity; end
831	if isfield(structmd,'max_nonlinear_iterations'), md.thermal.maxiter=structmd.max_nonlinear_iterations; end
832	if isfield(structmd,'stabilize_constraints'), md.thermal.penalty_lock=structmd.stabilize_constraints; end
833	if isfield(structmd,'penalty_offset'), md.thermal.penalty_factor=structmd.penalty_offset; end
834	if isfield(structmd,'name'), md.miscellaneous.name=structmd.name; end
835	if isfield(structmd,'notes'), md.miscellaneous.notes=structmd.notes; end
836	if isfield(structmd,'dummy'), md.miscellaneous.dummy=structmd.dummy; end
837	if isfield(structmd,'dt'), md.timestepping.time_step=structmd.dt; end
838	if isfield(structmd,'ndt'), md.timestepping.final_time=structmd.ndt; end
839	if isfield(structmd,'time_adapt'), md.timestepping.time_adapt=structmd.time_adapt; end
840	if isfield(structmd,'cfl_coefficient'), md.timestepping.cfl_coefficient=structmd.cfl_coefficient; end
841	if isfield(structmd,'spcthickness'), md.masstransport.spcthickness=structmd.spcthickness; end
842	if isfield(structmd,'artificial_diffusivity'), md.masstransport.stabilization=structmd.artificial_diffusivity; end
843	if isfield(structmd,'hydrostatic_adjustment'), md.masstransport.hydrostatic_adjustment=structmd.hydrostatic_adjustment; end
844	if isfield(structmd,'penalties'), md.masstransport.vertex_pairing=structmd.penalties; end
845	if isfield(structmd,'penalty_offset'), md.masstransport.penalty_factor=structmd.penalty_offset; end
846	if isfield(structmd,'B'), md.materials.rheology_B=structmd.B; end
847	if isfield(structmd,'n'), md.materials.rheology_n=structmd.n; end
848	if isfield(structmd,'rheology_B'), md.materials.rheology_B=structmd.rheology_B; end
849	if isfield(structmd,'rheology_n'), md.materials.rheology_n=structmd.rheology_n; end
850	if isfield(structmd,'rheology_Z'), md.damage.D=(1-structmd.rheology_Z); end
851	if isfield(structmd,'spcthickness'), md.balancethickness.spcthickness=structmd.spcthickness; end
852	if isfield(structmd,'artificial_diffusivity'), md.balancethickness.stabilization=structmd.artificial_diffusivity; end
853	if isfield(structmd,'dhdt'), md.balancethickness.thickening_rate=structmd.dhdt; end
854	if isfield(structmd,'isSIA'), md.flowequation.isSIA=structmd.isSIA; end
855	if isfield(structmd,'isFS'), md.flowequation.isFS=structmd.isFS; end
856	if isfield(structmd,'elements_type'), md.flowequation.element_equation=structmd.elements_type; end
857	if isfield(structmd,'vertices_type'), md.flowequation.vertex_equation=structmd.vertices_type; end
858	if isfield(structmd,'eps_rel'), md.steadystate.reltol=structmd.eps_rel; end
859	if isfield(structmd,'max_steadystate_iterations'), md.steadystate.maxiter=structmd.max_steadystate_iterations; end
860	if isfield(structmd,'isdiagnostic'), md.transient.isstressbalance=structmd.isdiagnostic; end
861	if isfield(structmd,'isprognostic'), md.transient.ismasstransport=structmd.isprognostic; end
862	if isfield(structmd,'isthermal'), md.transient.isthermal=structmd.isthermal; end
863	if isfield(structmd,'control_analysis'), md.inversion.iscontrol=structmd.control_analysis; end
864	if isfield(structmd,'weights'), md.inversion.cost_functions_coefficients=structmd.weights; end
865	if isfield(structmd,'nsteps'), md.inversion.nsteps=structmd.nsteps; end
866	if isfield(structmd,'maxiter_per_step'), md.inversion.maxiter_per_step=structmd.maxiter_per_step; end
867	if isfield(structmd,'cm_min'), md.inversion.min_parameters=structmd.cm_min; end
868	if isfield(structmd,'cm_max'), md.inversion.max_parameters=structmd.cm_max; end
869	if isfield(structmd,'vx_obs'), md.inversion.vx_obs=structmd.vx_obs; end
870	if isfield(structmd,'vy_obs'), md.inversion.vy_obs=structmd.vy_obs; end
871	if isfield(structmd,'vel_obs'), md.inversion.vel_obs=structmd.vel_obs; end
872	if isfield(structmd,'thickness_obs'), md.inversion.thickness_obs=structmd.thickness_obs; end
873	if isfield(structmd,'vx'), md.initialization.vx=structmd.vx; end
874	if isfield(structmd,'vy'), md.initialization.vy=structmd.vy; end
875	if isfield(structmd,'vz'), md.initialization.vz=structmd.vz; end
876	if isfield(structmd,'vel'), md.initialization.vel=structmd.vel; end
877	if isfield(structmd,'pressure'), md.initialization.pressure=structmd.pressure; end
878	if isfield(structmd,'temperature'), md.initialization.temperature=structmd.temperature; end
879	if isfield(structmd,'waterfraction'), md.initialization.waterfraction=structmd.waterfraction; end
880	if isfield(structmd,'watercolumn'), md.initialization.watercolumn=structmd.watercolumn; end
881	if isfield(structmd,'surface'), md.geometry.surface=structmd.surface; end
882	if isfield(structmd,'bed'), md.geometry.base=structmd.bed; end
883	if isfield(structmd,'thickness'), md.geometry.thickness=structmd.thickness; end
884	if isfield(structmd,'bathymetry'), md.geometry.bed=structmd.bathymetry; end
885	if isfield(structmd,'thickness_coeff'), md.geometry.hydrostatic_ratio=structmd.thickness_coeff; end
886	if isfield(structmd,'connectivity'), md.mesh.average_vertex_connectivity=structmd.connectivity; end
887	if isfield(structmd,'extractednodes'), md.mesh.extractedvertices=structmd.extractednodes; end
888	if isfield(structmd,'extractedelements'), md.mesh.extractedelements=structmd.extractedelements; end
889	if isfield(structmd,'nodeonboundary'), md.mesh.vertexonboundary=structmd.nodeonboundary; end
890	if isfield(structmd,'lat'), md.mesh.lat=structmd.lat; end
891	if isfield(structmd,'long'), md.mesh.long=structmd.long; end
892	if isfield(structmd,'segments'), md.mesh.segments=structmd.segments; end
893	if isfield(structmd,'segmentmarkers'), md.mesh.segmentmarkers=structmd.segmentmarkers; end
894	if isfield(structmd,'numlayers'), md.mesh.numberoflayers=structmd.numlayers; end
895	if isfield(structmd,'numberofelements'), md.mesh.numberofelements=structmd.numberofelements; end
896	if isfield(structmd,'numberofvertices'), md.mesh.numberofvertices=structmd.numberofvertices; end
897	if isfield(structmd,'numberofnodes'), md.mesh.numberofvertices=structmd.numberofnodes; end
898	if isfield(structmd,'numberofedges'), md.mesh.numberofedges=structmd.numberofedges; end
899	if isfield(structmd,'numberofelements2d'), md.mesh.numberofelements2d=structmd.numberofelements2d; end
900	if isfield(structmd,'numberofnodes2d'), md.mesh.numberofvertices2d=structmd.numberofnodes2d; end
901	if isfield(structmd,'nodeconnectivity'), md.mesh.vertexconnectivity=structmd.nodeconnectivity; end
902	if isfield(structmd,'elementconnectivity'), md.mesh.elementconnectivity=structmd.elementconnectivity; end
903	if isfield(structmd,'uppernodes'), md.mesh.uppervertex=structmd.uppernodes; end
904	if isfield(structmd,'lowernodes'), md.mesh.lowervertex=structmd.lowernodes; end
905	if isfield(structmd,'upperelements'), md.mesh.upperelements=structmd.upperelements; end
906	if isfield(structmd,'lowerelements'), md.mesh.lowerelements=structmd.lowerelements; end
907	if isfield(structmd,'nodeonsurface'), md.mesh.vertexonsurface=structmd.nodeonsurface; end
908	if isfield(structmd,'nodeonbase'), md.mesh.vertexonbase=structmd.nodeonbase; end
909	if isfield(structmd,'elements2d'), md.mesh.elements2d=structmd.elements2d; end
910	if isfield(structmd,'y2d'), md.mesh.y2d=structmd.y2d; end
911	if isfield(structmd,'x2d'), md.mesh.x2d=structmd.x2d; end
912	if isfield(structmd,'elements'), md.mesh.elements=structmd.elements; end
913	if isfield(structmd,'edges'),
914	md.mesh.edges=structmd.edges;
915	md.mesh.edges(isnan(md.mesh.edges))=-1;
916	end
917	if isfield(structmd,'y'), md.mesh.y=structmd.y; end
918	if isfield(structmd,'x'), md.mesh.x=structmd.x; end
919	if isfield(structmd,'z'), md.mesh.z=structmd.z; end
920	if isfield(structmd,'mask'), md.flaim.criterion=structmd.mask; end
921	if isfield(structmd,'diagnostic_ref'), md.stressbalance.referential=structmd.diagnostic_ref; end
922	if isfield(structmd,'npart'); md.qmu.numberofpartitions=structmd.npart; end
923	if isfield(structmd,'part'); md.qmu.partition=structmd.part; end
924
925	if isnumeric(md.verbose),
926	md.verbose=verbose;
927	end
928
929	if isfield(structmd,'spcvelocity'),
930	md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
931	md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
932	md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
933	pos=find(structmd.spcvelocity(:,1)); md.stressbalance.spcvx(pos)=structmd.spcvelocity(pos,4);
934	pos=find(structmd.spcvelocity(:,2)); md.stressbalance.spcvy(pos)=structmd.spcvelocity(pos,5);
935	pos=find(structmd.spcvelocity(:,3)); md.stressbalance.spcvz(pos)=structmd.spcvelocity(pos,6);
936	end
937	if isfield(structmd,'spcvx'),
938	md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
939	pos=find(~isnan(structmd.spcvx)); md.stressbalance.spcvx(pos)=structmd.spcvx(pos);
940	end
941	if isfield(structmd,'spcvy'),
942	md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
943	pos=find(~isnan(structmd.spcvy)); md.stressbalance.spcvy(pos)=structmd.spcvy(pos);
944	end
945	if isfield(structmd,'spcvz'),
946	md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
947	pos=find(~isnan(structmd.spcvz)); md.stressbalance.spcvz(pos)=structmd.spcvz(pos);
948	end
949	if isfield(structmd,'pressureload'),
950	if ~isempty(structmd.pressureload) & ismember(structmd.pressureload(end,end),[118 119 120]),
951	pos=find(structmd.pressureload(:,end)==120); md.stressbalance.icefront(pos,end)=0;
952	pos=find(structmd.pressureload(:,end)==118); md.stressbalance.icefront(pos,end)=1;
953	pos=find(structmd.pressureload(:,end)==119); md.stressbalance.icefront(pos,end)=2;
954	end
955	end
956	if isfield(structmd,'elements_type') & structmd.elements_type(end,end)>50,
957	pos=find(structmd.elements_type==59); md.flowequation.element_equation(pos,end)=0;
958	pos=find(structmd.elements_type==55); md.flowequation.element_equation(pos,end)=1;
959	pos=find(structmd.elements_type==56); md.flowequation.element_equation(pos,end)=2;
960	pos=find(structmd.elements_type==60); md.flowequation.element_equation(pos,end)=3;
961	pos=find(structmd.elements_type==62); md.flowequation.element_equation(pos,end)=4;
962	pos=find(structmd.elements_type==57); md.flowequation.element_equation(pos,end)=5;
963	pos=find(structmd.elements_type==58); md.flowequation.element_equation(pos,end)=6;
964	pos=find(structmd.elements_type==61); md.flowequation.element_equation(pos,end)=7;
965	end
966	if isfield(structmd,'vertices_type') & structmd.vertices_type(end,end)>50,
967	pos=find(structmd.vertices_type==59); md.flowequation.vertex_equation(pos,end)=0;
968	pos=find(structmd.vertices_type==55); md.flowequation.vertex_equation(pos,end)=1;
969	pos=find(structmd.vertices_type==56); md.flowequation.vertex_equation(pos,end)=2;
970	pos=find(structmd.vertices_type==60); md.flowequation.vertex_equation(pos,end)=3;
971	pos=find(structmd.vertices_type==62); md.flowequation.vertex_equation(pos,end)=4;
972	pos=find(structmd.vertices_type==57); md.flowequation.vertex_equation(pos,end)=5;
973	pos=find(structmd.vertices_type==58); md.flowequation.vertex_equation(pos,end)=6;
974	pos=find(structmd.vertices_type==61); md.flowequation.vertex_equation(pos,end)=7;
975	end
976	if isfield(structmd,'rheology_law') & isnumeric(structmd.rheology_law),
977	if (structmd.rheology_law==272), md.materials.rheology_law='None'; end
978	if (structmd.rheology_law==368), md.materials.rheology_law='Paterson'; end
979	if (structmd.rheology_law==369), md.materials.rheology_law='Arrhenius'; end
980	end
981	if isfield(structmd,'groundingline_migration') & isnumeric(structmd.groundingline_migration),
982	if (structmd.groundingline_migration==272), md.groundingline.migration='None'; end
983	if (structmd.groundingline_migration==273), md.groundingline.migration='AggressiveMigration'; end
984	if (structmd.groundingline_migration==274), md.groundingline.migration='SoftMigration'; end
985	end
986	if isfield(structmd,'control_type') & isnumeric(structmd.control_type),
987	if (structmd.control_type==143), md.inversion.control_parameters={'FrictionCoefficient'}; end
988	if (structmd.control_type==190), md.inversion.control_parameters={'RheologyBbar'}; end
989	if (structmd.control_type==147), md.inversion.control_parameters={'Thickeningrate'}; end
990	end
991	if isfield(structmd,'cm_responses') & ismember(structmd.cm_responses(end,end),[165:170 383 388 389]),
992	pos=find(structmd.cm_responses==166); md.inversion.cost_functions(pos)=101;
993	pos=find(structmd.cm_responses==167); md.inversion.cost_functions(pos)=102;
994	pos=find(structmd.cm_responses==168); md.inversion.cost_functions(pos)=103;
995	pos=find(structmd.cm_responses==169); md.inversion.cost_functions(pos)=104;
996	pos=find(structmd.cm_responses==170); md.inversion.cost_functions(pos)=105;
997	pos=find(structmd.cm_responses==165); md.inversion.cost_functions(pos)=201;
998	pos=find(structmd.cm_responses==389); md.inversion.cost_functions(pos)=501;
999	pos=find(structmd.cm_responses==388); md.inversion.cost_functions(pos)=502;
1000	pos=find(structmd.cm_responses==382); md.inversion.cost_functions(pos)=503;
1001	end
1002
1003	if isfield(structmd,'artificial_diffusivity') & structmd.artificial_diffusivity==2,
1004	md.thermal.stabilization=2;
1005	md.masstransport.stabilization=1;
1006	md.balancethickness.stabilization=1;
1007	end
1008	if isnumeric(md.masstransport.hydrostatic_adjustment)
1009	if md.masstransport.hydrostatic_adjustment==269,
1010	md.masstransport.hydrostatic_adjustment='Incremental';
1011	else
1012	md.masstransport.hydrostatic_adjustment='Absolute';
1013	end
1014	end
1015
1016	%New fields
1017	if ~isfield(structmd,'upperelements') & isa(md.mesh,'mesh3dprisms')
1018	md.mesh.upperelements=transpose(1:md.mesh.numberofelements)+md.mesh.numberofelements2d;
1019	md.mesh.upperelements(end-md.mesh.numberofelements2d+1:end)=NaN;
1020	end
1021	if ~isfield(structmd,'lowerelements') & isa(md.mesh,'mesh3dprisms')
1022	md.mesh.lowerelements=transpose(1:md.mesh.numberofelements)-md.mesh.numberofelements2d;
1023	md.mesh.lowerelements(1:md.mesh.numberofelements2d)=NaN;
1024	end
1025	if ~isfield(structmd,'diagnostic_ref');
1026	md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
1027	end
1028	if ~isfield(structmd,'loadingforce');
1029	md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
1030	end
1031
1032	%2013 August 9
1033	if isfield(structmd,'prognostic') & isa(structmd.prognostic,'prognostic'),
1034	disp('Recovering old prognostic class');
1035	md.masstransport=masstransport(structmd.prognostic);
1036	end
1037	%2013 August 9
1038	if isfield(structmd,'diagnostic') & (isa(structmd.diagnostic,'diagnostic') \|\| isa(structmd.diagnostic,'stressbalance')),
1039	disp('Recovering old diagnostic class');
1040	md.stressbalance=stressbalance(structmd.diagnostic);
1041	end
1042	end% }}}
1043	function md = setdefaultparameters(md) % {{{
1044
1045	%initialize subclasses
1046	md.mesh = mesh2d();
1047	md.mask = mask();
1048	md.constants = constants();
1049	md.geometry = geometry();
1050	md.initialization = initialization();
1051	md.smb = SMBforcing();
1052	md.basalforcings = basalforcings();
1053	md.friction = friction();
1054	md.rifts = rifts();
1055	md.timestepping = timestepping();
1056	md.groundingline = groundingline();
1057	md.materials = matice();
1058	md.damage = damage();
1059	md.flowequation = flowequation();
1060	md.debug = debug();
1061	md.verbose = verbose();
1062	md.settings = settings();
1063	md.toolkits = toolkits();
1064	md.cluster = generic();
1065	md.balancethickness = balancethickness();
1066	md.stressbalance = stressbalance();
1067	md.hydrology = hydrologyshreve();
1068	md.masstransport = masstransport();
1069	md.thermal = thermal();
1070	md.steadystate = steadystate();
1071	md.transient = transient();
1072	md.calving = calving();
1073	md.gia = gia();
1074	md.autodiff = autodiff();
1075	md.flaim = flaim();
1076	md.inversion = inversion();
1077	md.qmu = qmu();
1078	md.radaroverlay = radaroverlay();
1079	md.results = struct();
1080	md.outputdefinition = outputdefinition();
1081	md.miscellaneous = miscellaneous();
1082	md.private = private();
1083	end
1084	%}}}
1085	function md = tetras(md,varargin) % {{{
1086	%TETRAS - split 3d prismatic mesh into 3 tetrahedrons
1087	%
1088	% Usage:
1089	% md=tetra(md)
1090
1091	if ~isa(md.mesh,'mesh3dprisms')
1092	error('mesh is not a 3d prismatic mesh');
1093	end
1094
1095	%Initialize tetra mesh
1096	md.mesh=mesh3dtetras(md.mesh);
1097
1098	%Subdivision from Philipp Furnstahl (http://studierstube.icg.tugraz.at/thesis/fuernstahl_thesis.pdf)
1099	steiner = 0;
1100	nbv = md.mesh.numberofvertices;
1101	nbt = 3*md.mesh.numberofelements;
1102	elements = zeros(nbt,4);
1103	for i=1:md.mesh.numberofelements
1104	v1=md.mesh.elements(i,1); v2=md.mesh.elements(i,2); v3=md.mesh.elements(i,3);
1105	v4=md.mesh.elements(i,4); v5=md.mesh.elements(i,5); v6=md.mesh.elements(i,6);
1106	if(min(v2,v4)<min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)<min(v2,v6)),
1107	steiner = steiner+1; nbv = nbv+1; nbt = nbt+5; v7 = nbv;
1108	md.mesh.x=[md.mesh.x; mean(md.mesh.x(md.mesh.elements(i,:)))];
1109	md.mesh.y=[md.mesh.y; mean(md.mesh.y(md.mesh.elements(i,:)))];
1110	md.mesh.z=[md.mesh.z; mean(md.mesh.z(md.mesh.elements(i,:)))];
1111	elements(3*(i-1)+1,:) = [v1 v2 v3 v7];
1112	elements(3*(i-1)+2,:) = [v1 v2 v4 v7];
1113	elements(3*(i-1)+3,:) = [v2 v4 v5 v7];
1114	elements(end+1,:) = [v2 v3 v5 v7];
1115	elements(end+1,:) = [v3 v5 v6 v7];
1116	elements(end+1,:) = [v1 v3 v6 v7];
1117	elements(end+1,:) = [v1 v4 v6 v7];
1118	elements(end+1,:) = [v4 v5 v6 v7];
1119	elseif(min(v2,v4)<min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)>min(v2,v6)),
1120	elements(3*(i-1)+1,:) = [v1 v2 v4 v6];
1121	elements(3*(i-1)+2,:) = [v2 v4 v5 v6];
1122	elements(3*(i-1)+3,:) = [v1 v2 v3 v6];
1123	elseif(min(v2,v4)<min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)<min(v2,v6)),
1124	elements(3*(i-1)+1,:) = [v1 v2 v3 v4];
1125	elements(3*(i-1)+2,:) = [v2 v3 v4 v5];
1126	elements(3*(i-1)+3,:) = [v3 v4 v5 v6];
1127	elseif(min(v2,v4)<min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)>min(v2,v6)),
1128	elements(3*(i-1)+1,:) = [v1 v2 v3 v4];
1129	elements(3*(i-1)+2,:) = [v2 v4 v5 v6];
1130	elements(3*(i-1)+3,:) = [v2 v3 v4 v6];
1131	elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)<min(v2,v6)),
1132	elements(3*(i-1)+1,:) = [v1 v4 v5 v6];
1133	elements(3*(i-1)+2,:) = [v1 v2 v3 v5];
1134	elements(3*(i-1)+3,:) = [v1 v3 v5 v6];
1135	elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)>min(v2,v6)),
1136	elements(3*(i-1)+1,:) = [v1 v4 v5 v6];
1137	elements(3*(i-1)+2,:) = [v1 v2 v5 v6];
1138	elements(3*(i-1)+3,:) = [v1 v2 v3 v6];
1139	elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)<min(v2,v6)),
1140	elements(3*(i-1)+1,:) = [v1 v3 v4 v5];
1141	elements(3*(i-1)+2,:) = [v1 v2 v3 v5];
1142	elements(3*(i-1)+3,:) = [v3 v4 v5 v6];
1143	elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)<min(v2,v6)),
1144	elements(3*(i-1)+1,:) = [v1 v5 v6 v4];
1145	elements(3*(i-1)+2,:) = [v1 v2 v3 v5];
1146	elements(3*(i-1)+3,:) = [v5 v6 v3 v1];
1147	elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)>min(v2,v6)),
1148	steiner = steiner+1; nbv = nbv+1; nbt = nbt+5; v7 = nbv;
1149	md.mesh.x=[md.mesh.x; mean(md.mesh.x(md.mesh.elements(i,:)))];
1150	md.mesh.y=[md.mesh.y; mean(md.mesh.y(md.mesh.elements(i,:)))];
1151	md.mesh.z=[md.mesh.z; mean(md.mesh.z(md.mesh.elements(i,:)))];
1152	elements(3*(i-1)+1,:) = [v1 v2 v3 v7];
1153	elements(3*(i-1)+2,:) = [v1 v4 v5 v7];
1154	elements(3*(i-1)+3,:) = [v1 v2 v5 v7];
1155	elements(end+1,:) = [v2 v5 v6 v7];
1156	elements(end+1,:) = [v2 v3 v6 v7];
1157	elements(end+1,:) = [v3 v4 v6 v7];
1158	elements(end+1,:) = [v1 v3 v4 v7];
1159	elements(end+1,:) = [v4 v5 v6 v7];
1160	else
1161	error('Case not supported'); %not supposed to happen!
1162	end
1163	%Reorder elements to make sure they are direct
1164	for j=1:3
1165	element = elements(3*(i-1)+j,:);
1166	matrix = [md.mesh.x(element), md.mesh.y(element), md.mesh.z(element), ones(4,1)];
1167	if det(matrix)>0,
1168	elements(3*(i-1)+j,1)=element(2);
1169	elements(3*(i-1)+j,2)=element(1);
1170	end
1171	end
1172	end
1173	%%Split in 3 tetras
1174	%subelement1 = [1 2 3 5];
1175	%subelement2 = [4 6 5 1];
1176	%subelement3 = [5 6 3 1];
1177	%elements=[md.mesh.elements(:,subelement1);md.mesh.elements(:,subelement2);md.mesh.elements(:,subelement3)];
1178	if steiner==0,
1179	disp('No Steiner point required to split prismatic mesh into tets');
1180	else
1181	disp([num2str(steiner) ' Steiner points had to be included'])
1182	error('Steiner point not supported yet');
1183	end
1184
1185	pos_elements = repmat([1:md.mesh.numberofelements]',3,1);
1186
1187	md.mesh.elements=elements;
1188	md.mesh.numberofelements=size(elements,1);
1189
1190	%p and q (same deal, except for element that are on the bedrock: )
1191	if ~isnan(md.friction.p),
1192	md.friction.p=md.friction.p(pos_elements);
1193	md.friction.q=md.friction.q(pos_elements);
1194	end
1195
1196	%elementstype
1197	if ~isnan(md.flowequation.element_equation)
1198	oldelements_type=md.flowequation.element_equation;
1199	md.flowequation.element_equation=md.flowequation.element_equation(pos_elements);
1200	end
1201
1202	%connectivity
1203	md.mesh.elementconnectivity=NaN;
1204
1205	%materials
1206	if ~isnan(md.materials.rheology_n),
1207	md.materials.rheology_n=md.materials.rheology_n(pos_elements);
1208	end
1209
1210	%increase connectivity if less than 25:
1211	if md.mesh.average_vertex_connectivity<=25,
1212	md.mesh.average_vertex_connectivity=100;
1213	end
1214	end % }}}
1215	function disp(self) % {{{
1216	disp(sprintf('%19s: %-22s -- %s','mesh' ,['[1x1 ' class(self.mesh) ']'],'mesh properties'));
1217	disp(sprintf('%19s: %-22s -- %s','mask' ,['[1x1 ' class(self.mask) ']'],'defines grounded and floating elements'));
1218	disp(sprintf('%19s: %-22s -- %s','geometry' ,['[1x1 ' class(self.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...'));
1219	disp(sprintf('%19s: %-22s -- %s','constants' ,['[1x1 ' class(self.constants) ']'],'physical constants'));
1220	disp(sprintf('%19s: %-22s -- %s','smb' ,['[1x1 ' class(self.smb) ']'],'surface mass balance'));
1221	disp(sprintf('%19s: %-22s -- %s','basalforcings' ,['[1x1 ' class(self.basalforcings) ']'],'bed forcings'));
1222	disp(sprintf('%19s: %-22s -- %s','materials' ,['[1x1 ' class(self.materials) ']'],'material properties'));
1223	disp(sprintf('%19s: %-22s -- %s','damage' ,['[1x1 ' class(self.damage) ']'],'parameters for damage evolution solution'));
1224	disp(sprintf('%19s: %-22s -- %s','friction' ,['[1x1 ' class(self.friction) ']'],'basal friction/drag properties'));
1225	disp(sprintf('%19s: %-22s -- %s','flowequation' ,['[1x1 ' class(self.flowequation) ']'],'flow equations'));
1226	disp(sprintf('%19s: %-22s -- %s','timestepping' ,['[1x1 ' class(self.timestepping) ']'],'time stepping for transient models'));
1227	disp(sprintf('%19s: %-22s -- %s','initialization' ,['[1x1 ' class(self.initialization) ']'],'initial guess/state'));
1228	disp(sprintf('%19s: %-22s -- %s','rifts' ,['[1x1 ' class(self.rifts) ']'],'rifts properties'));
1229	disp(sprintf('%19s: %-22s -- %s','debug' ,['[1x1 ' class(self.debug) ']'],'debugging tools (valgrind, gprof)'));
1230	disp(sprintf('%19s: %-22s -- %s','verbose' ,['[1x1 ' class(self.verbose) ']'],'verbosity level in solve'));
1231	disp(sprintf('%19s: %-22s -- %s','settings' ,['[1x1 ' class(self.settings) ']'],'settings properties'));
1232	disp(sprintf('%19s: %-22s -- %s','toolkits' ,['[1x1 ' class(self.toolkits) ']'],'PETSc options for each solution'));
1233	disp(sprintf('%19s: %-22s -- %s','cluster' ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of cpus...)'));
1234	disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(self.balancethickness) ']'],'parameters for balancethickness solution'));
1235	disp(sprintf('%19s: %-22s -- %s','stressbalance' ,['[1x1 ' class(self.stressbalance) ']'],'parameters for stressbalance solution'));
1236	disp(sprintf('%19s: %-22s -- %s','groundingline' ,['[1x1 ' class(self.groundingline) ']'],'parameters for groundingline solution'));
1237	disp(sprintf('%19s: %-22s -- %s','hydrology' ,['[1x1 ' class(self.hydrology) ']'],'parameters for hydrology solution'));
1238	disp(sprintf('%19s: %-22s -- %s','masstransport' ,['[1x1 ' class(self.masstransport) ']'],'parameters for masstransport solution'));
1239	disp(sprintf('%19s: %-22s -- %s','thermal' ,['[1x1 ' class(self.thermal) ']'],'parameters for thermal solution'));
1240	disp(sprintf('%19s: %-22s -- %s','steadystate' ,['[1x1 ' class(self.steadystate) ']'],'parameters for steadystate solution'));
1241	disp(sprintf('%19s: %-22s -- %s','transient' ,['[1x1 ' class(self.transient) ']'],'parameters for transient solution'));
1242	disp(sprintf('%19s: %-22s -- %s','calving' ,['[1x1 ' class(self.calving) ']'],'parameters for calving'));
1243	disp(sprintf('%19s: %-22s -- %s','gia' ,['[1x1 ' class(self.gia) ']'],'parameters for gia solution'));
1244	disp(sprintf('%19s: %-22s -- %s','autodiff' ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters'));
1245	disp(sprintf('%19s: %-22s -- %s','flaim' ,['[1x1 ' class(self.flaim) ']'],'flaim parameters'));
1246	disp(sprintf('%19s: %-22s -- %s','inversion' ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods'));
1247	disp(sprintf('%19s: %-22s -- %s','qmu' ,['[1x1 ' class(self.qmu) ']'],'dakota properties'));
1248	disp(sprintf('%19s: %-22s -- %s','outputdefinition',['[1x1 ' class(self.outputdefinition) ']'],'output definition'));
1249	disp(sprintf('%19s: %-22s -- %s','results' ,['[1x1 ' class(self.results) ']'],'model results'));
1250	disp(sprintf('%19s: %-22s -- %s','radaroverlay' ,['[1x1 ' class(self.radaroverlay) ']'],'radar image for plot overlay'));
1251	disp(sprintf('%19s: %-22s -- %s','miscellaneous' ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields'));
1252	end % }}}
1253	function memory(self) % {{{
1254
1255	disp(sprintf('\nMemory imprint:\n'));
1256
1257	fields=properties('model');
1258	mem=0;
1259
1260	for i=1:length(fields),
1261	field=self.(fields{i});
1262	s=whos('field');
1263	mem=mem+s.bytes/1e6;
1264	disp(sprintf('%19s: %6.2f Mb',fields{i},s.bytes/1e6));
1265	end
1266	disp(sprintf('%19s--%10s','--------------','--------------'));
1267	disp(sprintf('%19s: %g Mb','Total',mem));
1268	end % }}}
1269	function netcdf(self,filename) % {{{
1270	%NETCDF - save model as netcdf
1271	%
1272	% Usage:
1273	% netcdf(md,filename)
1274	%
1275	% Example:
1276	% netcdf(md,'model.nc');
1277
1278	disp('Saving model as NetCDF');
1279	%1. Create NetCDF file
1280	ncid=netcdf.create(filename,'CLOBBER');
1281	netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Conventions','CF-1.4');
1282	netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Title',['ISSM model (' self.miscellaneous.name ')']);
1283	netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Author',getenv('USER'));
1284	netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Date',datestr(now));
1285
1286	%Preallocate variable id, needed to write variables in netcdf file
1287	var_id=zeros(1000,1);%preallocate
1288
1289	for step=1:2,
1290	counter=0;
1291	[var_id,counter]=structtonc(ncid,'md',self,0,var_id,counter,step);
1292	if step==1, netcdf.endDef(ncid); end
1293	end
1294
1295	if counter>1000,
1296	warning(['preallocation of var_id need to be updated from ' num2str(1000) ' to ' num2str(counter)]);
1297	end
1298
1299	netcdf.close(ncid)
1300	end % }}}
1301	function xylim(self) % {{{
1302
1303	xlim([min(self.mesh.x) max(self.mesh.x)]);
1304	ylim([min(self.mesh.y) max(self.mesh.y)])
1305	end % }}}
1306	function md=upload(md) % {{{
1307	%the goal of this routine is to upload the model onto a server, and to empty it.
1308	%So first, save the model with a unique name and upload the file to the server:
1309	random_part=fix(rand(1)*10000);
1310	id=[md.miscellaneous.name '-' regexprep(datestr(now),'[^\w'']','') '-' num2str(random_part) '-' getenv('USER') '-' oshostname() '.upload'];
1311	eval(['save ' id ' md']);
1312
1313	%Now, upload the file:
1314	issmscpout(md.settings.upload_server,md.settings.upload_path,md.settings.upload_login,md.settings.upload_port,{id},1);
1315
1316	%Now, empty this model of everything except settings, and record name of file we just uploaded!
1317	settings_back=md.settings;
1318	md=model();
1319	md.settings=settings_back;
1320	md.settings.upload_filename=id;
1321
1322	%get locally rid of file that was uploaded
1323	eval(['delete ' id]);
1324
1325	end % }}}
1326	function md=download(md) % {{{
1327
1328	%the goal of this routine is to download the internals of the current model from a server, because
1329	%this model is empty, except for the settings which tell us where to go and find this model!
1330
1331	%Download the file:
1332	issmscpin(md.settings.upload_server, md.settings.upload_login, md.settings.upload_port, md.settings.upload_path, {md.settings.upload_filename});
1333
1334	name=md.settings.upload_filename;
1335
1336	%Now, load this model:
1337	md=loadmodel(md.settings.upload_filename);
1338
1339	%get locally rid of file that was downloaded
1340	eval(['delete ' name]);
1341
1342	end % }}}
1343	end
1344	end

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: issm/trunk-jpl/src/m/classes/model.m@ 19527

Download in other formats: