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

Last change on this file since 18775 was 18775, checked in by bdef, 10 years ago

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