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

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

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