Context Navigation

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

Last change on this file since 24469 was 24469, checked in by Eric.Larour, 5 years ago
CHG: new dsl class to take care of the dynamic sea level in our sealevelrise_core. We made a new field to the model class to fit the dsl fields. The fields come from the slr steric rates class essentially. We also added GetStericRate and GetDynamicRate to the sea level core, which splits what used to be only the steric rate field.
File size: 69.9 KB

Rev	Line
[8926]	1	%MODEL class definition
	2	%
	3	% Usage:
	4	% md = model(varargin)
	5
	6	classdef 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;
[19527]	15	smb = 0;
[13692]	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;
[24469]	24	dsl = 0;
[19984]	25	slr = 0;
[9778]	26
[13692]	27	debug = 0;
	28	verbose = 0;
	29	settings = 0;
[14621]	30	toolkits = 0;
[13692]	31	cluster = 0;
[9778]	32
[13692]	33	balancethickness = 0;
[17757]	34	stressbalance = 0;
[13692]	35	groundingline = 0;
	36	hydrology = 0;
[17757]	37	masstransport = 0;
[13692]	38	thermal = 0;
	39	steadystate = 0;
	40	transient = 0;
[22958]	41	levelset = 0;
[18757]	42	calving = 0;
[23652]	43	frontalforcings = 0;
[22958]	44	love = 0;
	45	gia = 0;
[21260]	46	esa = 0;
[22958]	47
[13692]	48	autodiff = 0;
	49	inversion = 0;
	50	qmu = 0;
[22955]	51	amr = 0;
[13692]	52	results = 0;
[16388]	53	outputdefinition = 0;
[13692]	54	radaroverlay = 0;
	55	miscellaneous = 0;
	56	private = 0;
[9778]	57
[13692]	58	%}}}
	59	end
	60	methods (Static)
	61	function md = loadobj(md) % {{{
	62	% This function is directly called by matlab when a model object is
	63	% loaded. If the input is a struct it is an old version of model and
	64	% old fields must be recovered (make sure they are in the deprecated
	65	% model properties)
[8926]	66
[13692]	67	if verLessThan('matlab','7.9'),
	68	disp('Warning: your matlab version is old and there is a risk that load does not work correctly');
	69	disp(' if the model is not loaded correctly, rename temporarily loadobj so that matlab does not use it');
[8952]	70
[13692]	71	% This is a Matlab bug: all the fields of md have their default value
	72	% Example of error message:
	73	% Warning: Error loading an object of class 'model':
	74	% Undefined function or method 'exist' for input arguments of type 'cell'
	75	%
	76	% This has been fixed in MATLAB 7.9 (R2009b) and later versions
	77	end
[8952]	78
[13692]	79	if isstruct(md)
	80	disp('Recovering model object from a previous version');
	81	md = structtomodel(model,md);
	82	end
[13239]	83
[13692]	84	%2012 August 4th
	85	if isa(md.materials,'materials'),
	86	disp('Recovering old materials');
	87	if numel(md.materials.rheology_Z)==1 & isnan(md.materials.rheology_Z),
[13718]	88	md.materials=matice(md.materials);
	89	else
[13692]	90	md.materials=matdamageice(md.materials);
	91	end
	92	end
[14559]	93	%2013 April 12
[15771]	94	if numel(md.stressbalance.loadingforce==1)
	95	md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
[14559]	96	end
[14618]	97	%2013 April 17
	98	if isa(md.hydrology,'hydrology'),
	99	disp('Recovering old hydrology class');
	100	md.hydrology=hydrologyshreve(md.materials);
	101	end
[16356]	102	%2013 October 9
	103	if ~isa(md.damage,'damage'),
	104	md.damage=damage();
	105	md.damage.D=zeros(md.mesh.numberofvertices,1);
	106	md.damage.spcdamage=NaN*ones(md.mesh.numberofvertices,1);
	107	end
[16822]	108	%2013 November 18
	109	if ~isa(md.outputdefinition,'outputdefinition'),
	110	md.outputdefinition=outputdefinition();
	111	end
[17558]	112	%2014 March 26th
	113	if isa(md.mesh,'mesh'),
[17566]	114	disp('Recovering old mesh class');
	115	if isprop(md.mesh,'dimension'),
	116	if md.mesh.dimension==2,
	117	md.mesh=mesh2d(md.mesh);
	118	else
	119	md.mesh=mesh3dprisms(md.mesh);
	120	end
[17558]	121	else
[17566]	122	md.mesh=mesh2dvertical(md.mesh);
[17558]	123	end
	124	end
[18775]	125	%2014 November 12
[21702]	126	if isa(md.calving,'double'); md.calving=calving(); end
[20059]	127	%2016 February 3
[21702]	128	if isa(md.slr,'double'); md.slr=slr(); end
[21260]	129	%2016 October 11
[21702]	130	if isa(md.esa,'double'); md.esa=esa(); end
[21545]	131	%2017 February 10th
[23434]	132	if isa(md.settings,'settings'), %this 'isa' verification: 2018 October 24th
	133	if md.settings.solver_residue_threshold==0,
	134	md.settings.solver_residue_threshold = 1e-6;
	135	end
[21545]	136	end
[21667]	137	%2017 April 10th
[21702]	138	if isa(md.gia,'gia'), md.gia=giaivins(); end
	139	%2017 May 4th
	140	if isa(md.amr,'double'); md.amr=amr(); end
[22019]	141	%2017 Aug 29th
[22040]	142	if isa(md.love,'double'); md.love=fourierlove(); end
[22194]	143	%2017 Oct 26th
	144	if isa(md.calving,'calvingdev')
	145	disp('Warning: calvingdev is now calvingvonmises');
	146	md.calving=calvingvonmises(md.calving);
	147	end
[22955]	148	%2017 Dec 21st (needs to be here)
	149	if isempty(md.settings)
	150	disp('Warning: md.settings had to be reset, make sure to adjust md.settings.output_frequency and other fields');
	151	md.settings = issmsettings();
	152	end
[23503]	153	%2018 Dec 1st
	154	if md.settings.sb_coupling_frequency==0
	155	md.settings.sb_coupling_frequency=1;
	156	end
[23652]	157	%2019 Jan..
	158	if isa(md.frontalforcings,'double');
	159	if(~isnan(md.calving.meltingrate))
	160	disp('Warning: md.calving.meltingrate is now in md.frontalforcings');
	161	end
	162	md.frontalforcings=frontalforcings(md.calving);
	163	end
[23758]	164	%2019 Feb 26
	165	if isa(md.settings.results_on_nodes,'double')
	166	if md.settings.results_on_nodes == 0
	167	md.settings.results_on_nodes = {};
	168	else
	169	md.settings.results_on_nodes = {'all'};
	170	end
	171	end
[23814]	172	%2019 Mar 28
[23984]	173	if isa(md.smb,'SMBcomponents') \| isa(md.smb,'SMBmeltcomponents') \| isa(md.smb,'SMBforcing') \| isa(md.smb,'SMBgemb')
	174	if isa(md.smb.isclimatology,'double')
	175	if prod(size(md.smb.isclimatology)) ~= 1
	176	md.smb.isclimatology = 0;
	177	end
[23814]	178	end
	179	end
[13692]	180	end% }}}
	181	end
	182	methods
	183	function md = model(varargin) % {{{
[8926]	184
[13692]	185	switch nargin
	186	case 0
	187	md=setdefaultparameters(md);
[18492]	188	case 1
[19087]	189	error('model constructor not supported yet');
[18503]	190
[13692]	191	otherwise
	192	error('model constructor error message: 0 of 1 argument only in input.');
	193	end
	194	end
	195	%}}}
	196	function md = checkmessage(md,string) % {{{
	197	if(nargout~=1) error('wrong usage, model must be an output'); end
	198	disp(['model not consistent: ' string]);
	199	md.private.isconsistent=false;
	200	end
	201	%}}}
	202	function md = collapse(md)% {{{
	203	%COLLAPSE - collapses a 3d mesh into a 2d mesh
	204	%
	205	% This routine collapses a 3d model into a 2d model
	206	% and collapses all the fileds of the 3d model by
	207	% taking their depth-averaged values
	208	%
	209	% Usage:
	210	% md=collapse(md)
	211	%
	212	% See also: EXTRUDE, MODELEXTRACT
[13005]	213
[13692]	214	%Check that the model is really a 3d model
[17687]	215	if ~strcmp(md.mesh.elementtype(),'Penta'),
[13692]	216	error('collapse error message: only 3d mesh can be collapsed')
	217	end
[13005]	218
[17724]	219	%Start with changing all the fields from the 3d mesh
[13005]	220
[21827]	221	%dealing with the friction law
[13692]	222	%drag is limited to nodes that are on the bedrock.
[18775]	223	if isa(md.friction,'friction'),
	224	md.friction.coefficient=project2d(md,md.friction.coefficient,1);
	225	md.friction.p=project2d(md,md.friction.p,1);
	226	md.friction.q=project2d(md,md.friction.q,1);
[21819]	227	elseif isa(md.friction,'frictioncoulomb'),
	228	md.friction.coefficient=project2d(md,md.friction.coefficient,1);
	229	md.friction.coefficientcoulomb=project2d(md,md.friction.coefficientcoulomb,1);
	230	md.friction.p=project2d(md,md.friction.p,1);
	231	md.friction.q=project2d(md,md.friction.q,1);
[18775]	232	elseif isa(md.friction,'frictionhydro'),
	233	md.friction.q=project2d(md,md.friction.q,1);
	234	md.friction.C=project2d(md,md.friction.C,1);
	235	md.friction.As=project2d(md,md.friction.As,1);
[18798]	236	md.friction.effective_pressure=project2d(md,md.friction.effective_pressure,1);
[18775]	237	elseif isa(md.friction,'frictionwaterlayer'),
	238	md.friction.coefficient=project2d(md,md.friction.coefficient,1);
	239	md.friction.p=project2d(md,md.friction.p,1);
	240	md.friction.q=project2d(md,md.friction.q,1);
	241	md.friction.water_layer=project2d(md,md.friction.water_layer,1);
	242	elseif isa(md.friction,'frictionweertman'),
	243	md.friction.C=project2d(md,md.friction.C,1);
	244	md.friction.m=project2d(md,md.friction.m,1);
[19720]	245	elseif isa(md.friction,'frictionweertmantemp'),
	246	md.friction.C=project2d(md,md.friction.C,1);
	247	md.friction.m=project2d(md,md.friction.m,1);
	248	else
	249	disp('friction type not supported');
[21827]	250	end
[13005]	251
[13692]	252	%observations
	253	if ~isnan(md.inversion.vx_obs), md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); end;
	254	if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); end;
	255	if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end;
[24465]	256	if ~isnan(md.inversion.thickness_obs), md.inversion.thickness_obs=project2d(md,md.inversion.thickness_obs,md.mesh.numberoflayers); end;
[13692]	257	if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end;
	258	if numel(md.inversion.min_parameters)>1, md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end;
	259	if numel(md.inversion.max_parameters)>1, md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end;
[19527]	260	if isa(md.smb,'SMBforcing') & ~isnan(md.smb.mass_balance),
	261	md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers);
	262	elseif isa(md.smb,'SMBhenning') & ~isnan(md.smb.smbref),
	263	md.smb.smbref=project2d(md,md.smb.smbref,md.mesh.numberoflayers);
[13692]	264	end;
[13005]	265
[13692]	266	%results
	267	if ~isnan(md.initialization.vx),md.initialization.vx=DepthAverage(md,md.initialization.vx);end;
	268	if ~isnan(md.initialization.vy),md.initialization.vy=DepthAverage(md,md.initialization.vy);end;
	269	if ~isnan(md.initialization.vz),md.initialization.vz=DepthAverage(md,md.initialization.vz);end;
	270	if ~isnan(md.initialization.vel),md.initialization.vel=DepthAverage(md,md.initialization.vel);end;
[18578]	271	if ~isnan(md.initialization.temperature),md.initialization.temperature=DepthAverage(md,md.initialization.temperature);end;
[18506]	272	if ~isnan(md.initialization.pressure),md.initialization.pressure=project2d(md,md.initialization.pressure,1);end;
[18578]	273	if ~isnan(md.initialization.sediment_head),md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1);end;
	274	if ~isnan(md.initialization.epl_head),md.initialization.epl_head=project2d(md,md.initialization.epl_head,1);end;
	275	if ~isnan(md.initialization.epl_thickness),md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1);end;
[21417]	276	if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project2d(md,md.initialization.waterfraction,1);end;
	277	if ~isnan(md.initialization.watercolumn),md.initialization.watercolumn=project2d(md,md.initialization.watercolumn,1);end;
[21530]	278	%giaivins
[21584]	279	if ~isnan(md.gia.mantle_viscosity), md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1); end
	280	if ~isnan(md.gia.lithosphere_thickness), md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1); end
[15021]	281
[13692]	282	%elementstype
	283	if ~isnan(md.flowequation.element_equation)
	284	md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1);
	285	md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1);
[15564]	286	md.flowequation.borderSSA=project2d(md,md.flowequation.borderSSA,1);
	287	md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1);
	288	md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1);
[13692]	289	end
[13005]	290
[13692]	291	%boundary conditions
[15771]	292	md.stressbalance.spcvx=project2d(md,md.stressbalance.spcvx,md.mesh.numberoflayers);
	293	md.stressbalance.spcvy=project2d(md,md.stressbalance.spcvy,md.mesh.numberoflayers);
	294	md.stressbalance.spcvz=project2d(md,md.stressbalance.spcvz,md.mesh.numberoflayers);
	295	md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers);
	296	md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers);
[15767]	297	md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers);
[21544]	298	if numel(md.damage.spcdamage)>1, md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers); end
[21417]	299	if numel(md.levelset.spclevelset)>1, md.levelset.spclevelset=project2d(md,md.levelset.spclevelset,md.mesh.numberoflayers); end
[13692]	300	md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);
[13005]	301
[18506]	302	% Hydrologydc variables
	303	if isa(md.hydrology,'hydrologydc');
	304	md.hydrology.spcsediment_head=project2d(md,md.hydrology.spcsediment_head,1);
	305	md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1);
	306	md.hydrology.sediment_transmitivity=project2d(md,md.hydrology.sediment_transmitivity,1);
	307	md.hydrology.basal_moulin_input=project2d(md,md.hydrology.basal_moulin_input,1);
	308	if(md.hydrology.isefficientlayer==1)
	309	md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1);
[21808]	310	end
	311	end
[18506]	312
[13692]	313	%materials
	314	md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B);
	315	md.materials.rheology_n=project2d(md,md.materials.rheology_n,1);
[16160]	316
	317	%damage:
[17940]	318	if md.damage.isdamage,
	319	md.damage.D=DepthAverage(md,md.damage.D);
	320	end
[13005]	321
[13692]	322	%special for thermal modeling:
[18378]	323	if ~isnan(md.basalforcings.groundedice_melting_rate),
	324	md.basalforcings.groundedice_melting_rate=project2d(md,md.basalforcings.groundedice_melting_rate,1);
	325	end
[21417]	326	if isprop(md.basalforcings,'floatingice_melting_rate') & ~isnan(md.basalforcings.floatingice_melting_rate),
[18378]	327	md.basalforcings.floatingice_melting_rate=project2d(md,md.basalforcings.floatingice_melting_rate,1);
	328	end
[13692]	329	md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); %bedrock only gets geothermal flux
[13005]	330
[21417]	331	if isprop(md.calving,'coeff') & ~isnan(md.calving.coeff),
	332	md.calving.coeff=project2d(md,md.calving.coeff,1);
	333	end
[23652]	334	if isprop(md.frontalforcings,'meltingrate') & ~isnan(md.frontalforcings.meltingrate),
	335	md.frontalforcings.meltingrate=project2d(md,md.frontalforcings.meltingrate,1);
[21417]	336	end
	337
[13692]	338	%update of connectivity matrix
	339	md.mesh.average_vertex_connectivity=25;
[13005]	340
[13692]	341	%Collapse the mesh
	342	nodes2d=md.mesh.numberofvertices2d;
	343	elements2d=md.mesh.numberofelements2d;
[13005]	344
[13692]	345	%parameters
	346	md.geometry.surface=project2d(md,md.geometry.surface,1);
	347	md.geometry.thickness=project2d(md,md.geometry.thickness,1);
[17590]	348	md.geometry.base=project2d(md,md.geometry.base,1);
[18480]	349	if ~isnan(md.geometry.bed),
	350	md.geometry.bed=project2d(md,md.geometry.bed,1);
	351	end
[17991]	352
[18378]	353	if ~isnan(md.mask.groundedice_levelset),
	354	md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1);
	355	end
	356	if ~isnan(md.mask.ice_levelset),
	357	md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1);
	358	end
[13005]	359
[22955]	360	%lat long
	361	if numel(md.mesh.lat) ==md.mesh.numberofvertices, md.mesh.lat=project2d(md,md.mesh.lat,1); end
	362	if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end
	363
[21808]	364	%outputdefinitions
	365	for i=1:length(md.outputdefinition.definitions)
	366	if isobject(md.outputdefinition.definitions{i})
	367	%get subfields
	368	solutionsubfields=fields(md.outputdefinition.definitions{i});
	369	for j=1:length(solutionsubfields),
	370	field=md.outputdefinition.definitions{i}.(solutionsubfields{j});
	371	if length(field)==md.mesh.numberofvertices \| length(field)==md.mesh.numberofelements,
	372	md.outputdefinition.definitions{i}.(solutionsubfields{j})=project2d(md,md.outputdefinition.definitions{i}.(solutionsubfields{j}),1);
	373	end
	374	end
	375	end
	376	end
	377
[13692]	378	%Initialize with the 2d mesh
[17724]	379	mesh=mesh2d();
	380	mesh.x=md.mesh.x2d;
	381	mesh.y=md.mesh.y2d;
	382	mesh.numberofvertices=md.mesh.numberofvertices2d;
	383	mesh.numberofelements=md.mesh.numberofelements2d;
	384	mesh.elements=md.mesh.elements2d;
[21499]	385	if numel(md.mesh.lat) ==md.mesh.numberofvertices, mesh.lat=project2d(md,md.mesh.lat,1); end
	386	if numel(md.mesh.long)==md.mesh.numberofvertices, mesh.long=project2d(md,md.mesh.long,1); end
	387	mesh.epsg=md.mesh.epsg;
[22324]	388	if numel(md.mesh.scale_factor)==md.mesh.numberofvertices, mesh.scale_factor=project2d(md,md.mesh.scale_factor,1); end
[17991]	389	if ~isnan(md.mesh.vertexonboundary), mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1); end
	390	if ~isnan(md.mesh.elementconnectivity), mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1); end
[17724]	391	md.mesh=mesh;
[18738]	392	md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
	393	md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
[19955]	394	md.mesh.segments=contourenvelope(md.mesh);
[13005]	395
[13692]	396	end % }}}
[22955]	397	function md2 = extract(md,area,varargin) % {{{
[13692]	398	%extract - extract a model according to an Argus contour or flag list
	399	%
	400	% This routine extracts a submodel from a bigger model with respect to a given contour
	401	% md must be followed by the corresponding exp file or flags list
	402	% It can either be a domain file (argus type, .exp extension), or an array of element flags.
	403	% If user wants every element outside the domain to be
[15564]	404	% extract2d, add '~' to the name of the domain file (ex: '~HO.exp');
[13692]	405	% an empty string '' will be considered as an empty domain
	406	% a string 'all' will be considered as the entire domain
	407	%
	408	% Usage:
	409	% md2=extract(md,area);
	410	%
	411	% Examples:
	412	% md2=extract(md,'Domain.exp');
	413	%
	414	% See also: EXTRUDE, COLLAPSE
[13005]	415
[13692]	416	%copy model
	417	md1=md;
[13005]	418
[22955]	419	%recover optoins:
	420	options=pairoptions(varargin{:});
	421
[13692]	422	%some checks
[22955]	423	if ((nargin<2) \| (nargout~=1)),
[13692]	424	help extract
	425	error('extract error message: bad usage');
	426	end
[13005]	427
[13692]	428	%get elements that are inside area
	429	flag_elem=FlagElements(md1,area);
	430	if ~any(flag_elem),
	431	error('extracted model is empty');
	432	end
[13005]	433
[13692]	434	%kick out all elements with 3 dirichlets
[22955]	435	if getfieldvalue(options,'spccheck',1)
	436	spc_elem=find(~flag_elem);
	437	spc_node=sort(unique(md1.mesh.elements(spc_elem,:)));
	438	flag=ones(md1.mesh.numberofvertices,1);
	439	flag(spc_node)=0;
	440	pos=find(sum(flag(md1.mesh.elements),2)==0);
	441	flag_elem(pos)=0;
	442	end
[13005]	443
[13692]	444	%extracted elements and nodes lists
	445	pos_elem=find(flag_elem);
	446	pos_node=sort(unique(md1.mesh.elements(pos_elem,:)));
[13005]	447
[13692]	448	%keep track of some fields
	449	numberofvertices1=md1.mesh.numberofvertices;
	450	numberofelements1=md1.mesh.numberofelements;
	451	numberofvertices2=length(pos_node);
	452	numberofelements2=length(pos_elem);
	453	flag_node=zeros(numberofvertices1,1);
	454	flag_node(pos_node)=1;
[13005]	455
[13692]	456	%Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
	457	Pelem=zeros(numberofelements1,1);
	458	Pelem(pos_elem)=[1:numberofelements2]';
	459	Pnode=zeros(numberofvertices1,1);
	460	Pnode(pos_node)=[1:numberofvertices2]';
[13005]	461
[13857]	462	%renumber the elements (some nodes won't exist anymore)
[13692]	463	elements_1=md1.mesh.elements;
	464	elements_2=elements_1(pos_elem,:);
	465	elements_2(:,1)=Pnode(elements_2(:,1));
	466	elements_2(:,2)=Pnode(elements_2(:,2));
	467	elements_2(:,3)=Pnode(elements_2(:,3));
[17558]	468	if isa(md1.mesh,'mesh3dprisms'),
[13692]	469	elements_2(:,4)=Pnode(elements_2(:,4));
	470	elements_2(:,5)=Pnode(elements_2(:,5));
	471	elements_2(:,6)=Pnode(elements_2(:,6));
	472	end
[13005]	473
[13857]	474	%OK, now create the new model!
[13005]	475
[13857]	476	%take every field from model
[13692]	477	md2=md1;
[13005]	478
[13692]	479	%automatically modify fields
[13005]	480
[13692]	481	%loop over model fields
	482	model_fields=fields(md1);
	483	for i=1:length(model_fields),
	484	%get field
	485	field=md1.(model_fields{i});
	486	fieldsize=size(field);
	487	if isobject(field), %recursive call
	488	object_fields=fields(md1.(model_fields{i}));
	489	for j=1:length(object_fields),
	490	%get field
	491	field=md1.(model_fields{i}).(object_fields{j});
	492	fieldsize=size(field);
	493	%size = number of nodes * n
	494	if fieldsize(1)==numberofvertices1
	495	md2.(model_fields{i}).(object_fields{j})=field(pos_node,:);
	496	elseif (fieldsize(1)==numberofvertices1+1)
	497	md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)];
[13857]	498	%size = number of elements * n
[13692]	499	elseif fieldsize(1)==numberofelements1
	500	md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:);
[21428]	501	elseif (fieldsize(1)==numberofelements1+1)
	502	md2.(model_fields{i}).(object_fields{j})=[field(pos_elem,:); field(end,:)];
[13692]	503	end
	504	end
	505	else
	506	%size = number of nodes * n
	507	if fieldsize(1)==numberofvertices1
	508	md2.(model_fields{i})=field(pos_node,:);
	509	elseif (fieldsize(1)==numberofvertices1+1)
	510	md2.(model_fields{i})=[field(pos_node,:); field(end,:)];
[13857]	511	%size = number of elements * n
[13692]	512	elseif fieldsize(1)==numberofelements1
	513	md2.(model_fields{i})=field(pos_elem,:);
[21428]	514	elseif (fieldsize(1)==numberofelements1+1)
	515	md2.(model_fields{i})=[field(pos_elem,:); field(end,:)];
[13692]	516	end
	517	end
	518	end
[13005]	519
[13692]	520	%modify some specific fields
[13005]	521
[13692]	522	%Mesh
	523	md2.mesh.numberofelements=numberofelements2;
	524	md2.mesh.numberofvertices=numberofvertices2;
	525	md2.mesh.elements=elements_2;
[13005]	526
[13692]	527	%mesh.uppervertex mesh.lowervertex
[17558]	528	if isa(md1.mesh,'mesh3dprisms'),
[13692]	529	md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node);
	530	pos=find(~isnan(md2.mesh.uppervertex));
	531	md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos));
[13005]	532
[13692]	533	md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node);
	534	pos=find(~isnan(md2.mesh.lowervertex));
	535	md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos));
[13005]	536
[13692]	537	md2.mesh.upperelements=md1.mesh.upperelements(pos_elem);
	538	pos=find(~isnan(md2.mesh.upperelements));
	539	md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos));
[13005]	540
[13692]	541	md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem);
	542	pos=find(~isnan(md2.mesh.lowerelements));
	543	md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos));
	544	end
[13005]	545
[13692]	546	%Initial 2d mesh
[17558]	547	if isa(md1.mesh,'mesh3dprisms'),
[13692]	548	flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d);
	549	pos_elem_2d=find(flag_elem_2d);
	550	flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d);
	551	pos_node_2d=find(flag_node_2d);
[13005]	552
[13692]	553	md2.mesh.numberofelements2d=length(pos_elem_2d);
	554	md2.mesh.numberofvertices2d=length(pos_node_2d);
	555	md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:);
	556	md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1));
	557	md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2));
	558	md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3));
[13005]	559
[13692]	560	md2.mesh.x2d=md1.mesh.x(pos_node_2d);
	561	md2.mesh.y2d=md1.mesh.y(pos_node_2d);
	562	end
[13005]	563
[13692]	564	%Edges
[17686]	565	if(dimension(md.mesh)==2),
[17563]	566	if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs...
	567	%renumber first two columns
	568	pos=find(md2.mesh.edges(:,4)~=-1);
	569	md2.mesh.edges(: ,1)=Pnode(md2.mesh.edges(:,1));
	570	md2.mesh.edges(: ,2)=Pnode(md2.mesh.edges(:,2));
	571	md2.mesh.edges(: ,3)=Pelem(md2.mesh.edges(:,3));
	572	md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4));
	573	%remove edges when the 2 vertices are not in the domain.
	574	md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:);
	575	%Replace all zeros by -1 in the last two columns
	576	pos=find(md2.mesh.edges(:,3)==0);
	577	md2.mesh.edges(pos,3)=-1;
	578	pos=find(md2.mesh.edges(:,4)==0);
	579	md2.mesh.edges(pos,4)=-1;
	580	%Invert -1 on the third column with last column (Also invert first two columns!!)
	581	pos=find(md2.mesh.edges(:,3)==-1);
	582	md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4);
	583	md2.mesh.edges(pos,4)=-1;
	584	values=md2.mesh.edges(pos,2);
	585	md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1);
	586	md2.mesh.edges(pos,1)=values;
	587	%Finally remove edges that do not belong to any element
	588	pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1);
	589	md2.mesh.edges(pos,:)=[];
	590	end
[13692]	591	end
[13005]	592
[13692]	593	%Penalties
[15771]	594	if ~isnan(md2.stressbalance.vertex_pairing),
	595	for i=1:size(md1.stressbalance.vertex_pairing,1);
	596	md2.stressbalance.vertex_pairing(i,:)=Pnode(md1.stressbalance.vertex_pairing(i,:));
[13692]	597	end
[15771]	598	md2.stressbalance.vertex_pairing=md2.stressbalance.vertex_pairing(find(md2.stressbalance.vertex_pairing(:,1)),:);
[13692]	599	end
[15767]	600	if ~isnan(md2.masstransport.vertex_pairing),
	601	for i=1:size(md1.masstransport.vertex_pairing,1);
	602	md2.masstransport.vertex_pairing(i,:)=Pnode(md1.masstransport.vertex_pairing(i,:));
[13692]	603	end
[15767]	604	md2.masstransport.vertex_pairing=md2.masstransport.vertex_pairing(find(md2.masstransport.vertex_pairing(:,1)),:);
[13692]	605	end
[13005]	606
[13692]	607	%recreate segments
[20322]	608	if isa(md1.mesh,'mesh2d') \| isa(md1.mesh','mesh3dsurface'),
[13692]	609	md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices);
	610	md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
[19957]	611	md2.mesh.segments=contourenvelope(md2.mesh);
[13692]	612	md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
	613	else
	614	%First do the connectivity for the contourenvelope in 2d
	615	md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d);
	616	md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity);
[19957]	617	segments=contourenvelope(md2.mesh);
[17565]	618	md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(segments(:,1:2))=1;
[13692]	619	md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1);
	620	%Then do it for 3d as usual
	621	md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices);
	622	md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
	623	end
[13005]	624
[13692]	625	%Boundary conditions: Dirichlets on new boundary
	626	%Catch the elements that have not been extracted
	627	orphans_elem=find(~flag_elem);
	628	orphans_node=unique(md1.mesh.elements(orphans_elem,:))';
	629	%Figure out which node are on the boundary between md2 and md1
	630	nodestoflag1=intersect(orphans_node,pos_node);
	631	nodestoflag2=Pnode(nodestoflag1);
[23771]	632	if numel(md1.stressbalance.spcvx)>1 & numel(md1.stressbalance.spcvy)>1 & numel(md1.stressbalance.spcvz)>1,
[13692]	633	if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1
[15771]	634	md2.stressbalance.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2);
	635	md2.stressbalance.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2);
[13692]	636	else
[15771]	637	md2.stressbalance.spcvx(nodestoflag2)=NaN;
	638	md2.stressbalance.spcvy(nodestoflag2)=NaN;
[13692]	639	disp(' ')
	640	disp('!! extract warning: spc values should be checked !!')
	641	disp(' ')
	642	end
	643	%put 0 for vz
[15771]	644	md2.stressbalance.spcvz(nodestoflag2)=0;
[13692]	645	end
	646	if ~isnan(md1.thermal.spctemperature),
	647	md2.thermal.spctemperature(nodestoflag2,1)=1;
	648	end
[13005]	649
[13692]	650	%Results fields
	651	if isstruct(md1.results),
	652	md2.results=struct();
	653	solutionfields=fields(md1.results);
	654	for i=1:length(solutionfields),
[14230]	655	if isstruct(md1.results.(solutionfields{i}))
	656	%get subfields
[24134]	657	% loop over time steps
	658	for p=1:length(md1.results.(solutionfields{i}))
	659	current = md1.results.(solutionfields{i})(p);
	660	solutionsubfields=fields(current);
	661	for j=1:length(solutionsubfields),
	662	field=md1.results.(solutionfields{i})(p).(solutionsubfields{j});
[14230]	663	if length(field)==numberofvertices1,
[24134]	664	md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_node);
[14230]	665	elseif length(field)==numberofelements1,
[24134]	666	md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_elem);
[14230]	667	else
[24134]	668	md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field;
[14230]	669	end
[24134]	670	end
[14230]	671	end
	672	else
	673	field=md1.results.(solutionfields{i});
[13692]	674	if length(field)==numberofvertices1,
[14230]	675	md2.results.(solutionfields{i})=field(pos_node);
[13692]	676	elseif length(field)==numberofelements1,
[14230]	677	md2.results.(solutionfields{i})=field(pos_elem);
[13692]	678	else
[14230]	679	md2.results.(solutionfields{i})=field;
[13692]	680	end
	681	end
	682	end
	683	end
[13005]	684
[21808]	685	%OutputDefinitions fields
	686	for i=1:length(md1.outputdefinition.definitions),
	687	if isobject(md1.outputdefinition.definitions{i})
	688	%get subfields
	689	solutionsubfields=fields(md1.outputdefinition.definitions{i});
	690	for j=1:length(solutionsubfields),
	691	field=md1.outputdefinition.definitions{i}.(solutionsubfields{j});
	692	if length(field)==numberofvertices1,
	693	md2.outputdefinition.definitions{i}.(solutionsubfields{j})=field(pos_node);
	694	elseif length(field)==numberofelements1,
	695	md2.outputdefinition.definitions{i}.(solutionsubfields{j})=field(pos_elem);
	696	end
	697	end
	698	end
	699	end
	700
[13692]	701	%Keep track of pos_node and pos_elem
	702	md2.mesh.extractedvertices=pos_node;
	703	md2.mesh.extractedelements=pos_elem;
	704	end % }}}
	705	function md = extrude(md,varargin) % {{{
	706	%EXTRUDE - vertically extrude a 2d mesh
	707	%
	708	% vertically extrude a 2d mesh and create corresponding 3d mesh.
	709	% The vertical distribution can:
	710	% - follow a polynomial law
	711	% - follow two polynomial laws, one for the lower part and one for the upper part of the mesh
	712	% - be discribed by a list of coefficients (between 0 and 1)
	713	%
	714	%
	715	% Usage:
	716	% md=extrude(md,numlayers,extrusionexponent);
	717	% md=extrude(md,numlayers,lowerexponent,upperexponent);
	718	% md=extrude(md,listofcoefficients);
	719	%
	720	% Example:
[18216]	721	% md=extrude(md,15,1.3);
	722	% md=extrude(md,15,1.3,1.2);
[13692]	723	% md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]);
	724	%
	725	% See also: MODELEXTRACT, COLLAPSE
[13005]	726
[13692]	727	%some checks on list of arguments
	728	if ((nargin>4) \| (nargin<2) \| (nargout~=1)),
	729	help extrude;
	730	error('extrude error message');
	731	end
[13005]	732
[13692]	733	%Extrude the mesh
	734	if nargin==2, %list of coefficients
	735	clist=varargin{1};
	736	if any(clist<0) \| any(clist>1),
	737	error('extrusioncoefficients must be between 0 and 1');
	738	end
	739	extrusionlist=sort(unique([clist(:);0;1]));
	740	numlayers=length(extrusionlist);
	741	elseif nargin==3, %one polynomial law
	742	if varargin{2}<=0,
	743	help extrude;
	744	error('extrusionexponent must be >=0');
	745	end
	746	numlayers=varargin{1};
	747	extrusionlist=((0:1:numlayers-1)/(numlayers-1)).^varargin{2};
	748	elseif nargin==4, %two polynomial laws
	749	numlayers=varargin{1};
	750	lowerexp=varargin{2};
	751	upperexp=varargin{3};
[13005]	752
[13692]	753	if varargin{2}<=0 \| varargin{3}<=0,
	754	help extrude;
	755	error('lower and upper extrusionexponents must be >=0');
	756	end
[13005]	757
[13692]	758	lowerextrusionlist=[(0:2/(numlayers-1):1).^lowerexp]/2;
	759	upperextrusionlist=[(0:2/(numlayers-1):1).^upperexp]/2;
	760	extrusionlist=sort(unique([lowerextrusionlist 1-upperextrusionlist]));
[13005]	761
[13692]	762	end
[13005]	763
[13692]	764	if numlayers<2,
	765	error('number of layers should be at least 2');
	766	end
[17686]	767	if strcmp(md.mesh.domaintype(),'3D')
[13692]	768	error('Cannot extrude a 3d mesh (extrude cannot be called more than once)');
	769	end
[13005]	770
[13692]	771	%Initialize with the 2d mesh
[17558]	772	mesh2d = md.mesh;
	773	md.mesh=mesh3dprisms();
	774	md.mesh.x = mesh2d.x;
	775	md.mesh.y = mesh2d.y;
	776	md.mesh.elements = mesh2d.elements;
	777	md.mesh.numberofelements = mesh2d.numberofelements;
	778	md.mesh.numberofvertices = mesh2d.numberofvertices;
	779
	780	md.mesh.lat = mesh2d.lat;
	781	md.mesh.long = mesh2d.long;
[18558]	782	md.mesh.epsg = mesh2d.epsg;
[22324]	783	md.mesh.scale_factor = mesh2d.scale_factor;
[17558]	784
	785	md.mesh.vertexonboundary = mesh2d.vertexonboundary;
	786	md.mesh.vertexconnectivity = mesh2d.vertexconnectivity;
	787	md.mesh.elementconnectivity = mesh2d.elementconnectivity;
	788	md.mesh.average_vertex_connectivity = mesh2d.average_vertex_connectivity;
	789
	790	md.mesh.extractedvertices = mesh2d.extractedvertices;
	791	md.mesh.extractedelements = mesh2d.extractedelements;
	792
[13692]	793	x3d=[];
	794	y3d=[];
	795	z3d=[]; %the lower node is on the bed
	796	thickness3d=md.geometry.thickness; %thickness and bed for these nodes
[17590]	797	bed3d=md.geometry.base;
[13005]	798
[13692]	799	%Create the new layers
	800	for i=1:numlayers,
	801	x3d=[x3d; md.mesh.x];
	802	y3d=[y3d; md.mesh.y];
	803	%nodes are distributed between bed and surface accordingly to the given exponent
	804	z3d=[z3d; bed3d+thickness3d*extrusionlist(i)];
	805	end
	806	number_nodes3d=size(x3d,1); %number of 3d nodes for the non extruded part of the mesh
[13005]	807
[13692]	808	%Extrude elements
	809	elements3d=[];
	810	for i=1:numlayers-1,
	811	elements3d=[elements3d;[md.mesh.elements+(i-1)md.mesh.numberofvertices md.mesh.elements+imd.mesh.numberofvertices]]; %Create the elements of the 3d mesh for the non extruded part
	812	end
	813	number_el3d=size(elements3d,1); %number of 3d nodes for the non extruded part of the mesh
[13005]	814
[13692]	815	%Keep a trace of lower and upper nodes
[17590]	816	lowervertex=NaN*ones(number_nodes3d,1);
	817	uppervertex=NaN*ones(number_nodes3d,1);
	818	lowervertex(md.mesh.numberofvertices+1:end)=1:(numlayers-1)*md.mesh.numberofvertices;
	819	uppervertex(1:(numlayers-1)*md.mesh.numberofvertices)=md.mesh.numberofvertices+1:number_nodes3d;
	820	md.mesh.lowervertex=lowervertex;
	821	md.mesh.uppervertex=uppervertex;
[13005]	822
[13692]	823	%same for lower and upper elements
[17590]	824	lowerelements=NaN*ones(number_el3d,1);
	825	upperelements=NaN*ones(number_el3d,1);
	826	lowerelements(md.mesh.numberofelements+1:end)=1:(numlayers-2)*md.mesh.numberofelements;
	827	upperelements(1:(numlayers-2)md.mesh.numberofelements)=md.mesh.numberofelements+1:(numlayers-1)md.mesh.numberofelements;
	828	md.mesh.lowerelements=lowerelements;
	829	md.mesh.upperelements=upperelements;
[13005]	830
[13692]	831	%Save old mesh
	832	md.mesh.x2d=md.mesh.x;
	833	md.mesh.y2d=md.mesh.y;
	834	md.mesh.elements2d=md.mesh.elements;
	835	md.mesh.numberofelements2d=md.mesh.numberofelements;
	836	md.mesh.numberofvertices2d=md.mesh.numberofvertices;
[13005]	837
[13692]	838	%Build global 3d mesh
	839	md.mesh.elements=elements3d;
	840	md.mesh.x=x3d;
	841	md.mesh.y=y3d;
	842	md.mesh.z=z3d;
	843	md.mesh.numberofelements=number_el3d;
	844	md.mesh.numberofvertices=number_nodes3d;
	845	md.mesh.numberoflayers=numlayers;
[13005]	846
[13692]	847	%Ok, now deal with the other fields from the 2d mesh:
[13005]	848
[19048]	849	%bedinfo and surface info
	850	md.mesh.vertexonbase=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1);
	851	md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers);
	852	md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node');
	853
[13692]	854	%lat long
	855	md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node');
	856	md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node');
[22324]	857	md.mesh.scale_factor=project3d(md,'vector',md.mesh.scale_factor,'type','node');
[13005]	858
[19048]	859	md.geometry=extrude(md.geometry,md);
	860	md.friction = extrude(md.friction,md);
	861	md.inversion = extrude(md.inversion,md);
[19527]	862	md.smb = extrude(md.smb,md);
[19048]	863	md.initialization = extrude(md.initialization,md);
[13005]	864
[19050]	865	md.flowequation=md.flowequation.extrude(md);
[19048]	866	md.stressbalance=extrude(md.stressbalance,md);
[19050]	867	md.thermal=md.thermal.extrude(md);
	868	md.masstransport=md.masstransport.extrude(md);
[20460]	869	md.levelset=extrude(md.levelset,md);
[19048]	870	md.calving=extrude(md.calving,md);
[23652]	871	md.frontalforcings=extrude(md.frontalforcings,md);
[19048]	872	md.hydrology = extrude(md.hydrology,md);
[22955]	873	md.slr = extrude(md.slr,md);
[24469]	874	md.dsl = extrude(md.dsl,md);
[13005]	875
[13692]	876	%connectivity
[17991]	877	if ~isnan(md.mesh.elementconnectivity)
	878	md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1);
	879	md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN;
	880	for i=2:numlayers-1,
	881	md.mesh.elementconnectivity((i-1)md.mesh.numberofelements2d+1:(i)md.mesh.numberofelements2d,:)...
	882	=md.mesh.elementconnectivity((i-1)md.mesh.numberofelements2d+1:(i)md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d;
	883	end
	884	md.mesh.elementconnectivity(find(isnan(md.mesh.elementconnectivity)))=0;
[13692]	885	end
[13005]	886
[19048]	887	md.materials=extrude(md.materials,md);
	888	md.damage=extrude(md.damage,md);
	889	md.mask=extrude(md.mask,md);
	890	md.qmu=extrude(md.qmu,md);
	891	md.basalforcings=extrude(md.basalforcings,md);
[21808]	892	md.outputdefinition=extrude(md.outputdefinition,md);
[13005]	893
[13692]	894	%increase connectivity if less than 25:
	895	if md.mesh.average_vertex_connectivity<=25,
	896	md.mesh.average_vertex_connectivity=100;
	897	end
[13005]	898	end % }}}
[13692]	899	function md = structtomodel(md,structmd) % {{{
[8952]	900
[13692]	901	if ~isstruct(structmd) error('input model is not a structure'); end
[8952]	902
[13692]	903	%loaded model is a struct, initialize output and recover all fields
	904	md = structtoobj(model,structmd);
[8952]	905
[13692]	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
[10452]	909
[13692]	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
[18378]	914	if isfield(structmd,'melting'), md.basalforcings.floatingice_melting_rate=structmd.melting; end
[18068]	915	if isfield(structmd,'melting_rate'), md.basalforcings.floatingice_melting_rate=structmd.melting_rate; end
[18378]	916	if isfield(structmd,'melting_rate'), md.basalforcings.groundedice_melting_rate=structmd.melting_rate; end
[19527]	917	if isfield(structmd,'accumulation'), md.smb.mass_balance=structmd.accumulation; end
[13692]	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
[17610]	922	if isfield(structmd,'gridonbase'), md.mesh.vertexonbase=structmd.gridonbase; end
[13692]	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
[14621]	926	if isfield(structmd,'petscoptions') & ~isempty(structmd.petscoptions), md.toolkits=structmd.petscoptions; end
[13692]	927	if isfield(structmd,'g'), md.constants.g=structmd.g; end
	928	if isfield(structmd,'yts'), md.constants.yts=structmd.yts; end
[19527]	929	if isfield(structmd,'surface_mass_balance'), md.smb.mass_balance=structmd.surface_mass_balance; end
[18068]	930	if isfield(structmd,'basal_melting_rate'), md.basalforcings.floatingice_melting_rate=structmd.basal_melting_rate; end
[13692]	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
[15767]	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
[13692]	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
[16160]	973	if isfield(structmd,'rheology_Z'), md.damage.D=(1-structmd.rheology_Z); end
[13692]	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
[15564]	977	if isfield(structmd,'isSIA'), md.flowequation.isSIA=structmd.isSIA; end
	978	if isfield(structmd,'isFS'), md.flowequation.isFS=structmd.isFS; end
[13692]	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
[15771]	983	if isfield(structmd,'isdiagnostic'), md.transient.isstressbalance=structmd.isdiagnostic; end
[15768]	984	if isfield(structmd,'isprognostic'), md.transient.ismasstransport=structmd.isprognostic; end
[13692]	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
[17590]	1005	if isfield(structmd,'bed'), md.geometry.base=structmd.bed; end
[13692]	1006	if isfield(structmd,'thickness'), md.geometry.thickness=structmd.thickness; end
[17590]	1007	if isfield(structmd,'bathymetry'), md.geometry.bed=structmd.bathymetry; end
[13692]	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
[22324]	1015	if isfield(structmd,'scale_factor'), md.mesh.scale_factor=structmd.scale_factor; end
[13692]	1016	if isfield(structmd,'segments'), md.mesh.segments=structmd.segments; end
	1017	if isfield(structmd,'segmentmarkers'), md.mesh.segmentmarkers=structmd.segmentmarkers; end
	1018	if isfield(structmd,'numlayers'), md.mesh.numberoflayers=structmd.numlayers; end
	1019	if isfield(structmd,'numberofelements'), md.mesh.numberofelements=structmd.numberofelements; end
	1020	if isfield(structmd,'numberofvertices'), md.mesh.numberofvertices=structmd.numberofvertices; end
	1021	if isfield(structmd,'numberofnodes'), md.mesh.numberofvertices=structmd.numberofnodes; end
	1022	if isfield(structmd,'numberofedges'), md.mesh.numberofedges=structmd.numberofedges; end
	1023	if isfield(structmd,'numberofelements2d'), md.mesh.numberofelements2d=structmd.numberofelements2d; end
	1024	if isfield(structmd,'numberofnodes2d'), md.mesh.numberofvertices2d=structmd.numberofnodes2d; end
	1025	if isfield(structmd,'nodeconnectivity'), md.mesh.vertexconnectivity=structmd.nodeconnectivity; end
	1026	if isfield(structmd,'elementconnectivity'), md.mesh.elementconnectivity=structmd.elementconnectivity; end
	1027	if isfield(structmd,'uppernodes'), md.mesh.uppervertex=structmd.uppernodes; end
	1028	if isfield(structmd,'lowernodes'), md.mesh.lowervertex=structmd.lowernodes; end
	1029	if isfield(structmd,'upperelements'), md.mesh.upperelements=structmd.upperelements; end
	1030	if isfield(structmd,'lowerelements'), md.mesh.lowerelements=structmd.lowerelements; end
	1031	if isfield(structmd,'nodeonsurface'), md.mesh.vertexonsurface=structmd.nodeonsurface; end
[17610]	1032	if isfield(structmd,'nodeonbase'), md.mesh.vertexonbase=structmd.nodeonbase; end
[13692]	1033	if isfield(structmd,'elements2d'), md.mesh.elements2d=structmd.elements2d; end
	1034	if isfield(structmd,'y2d'), md.mesh.y2d=structmd.y2d; end
	1035	if isfield(structmd,'x2d'), md.mesh.x2d=structmd.x2d; end
	1036	if isfield(structmd,'elements'), md.mesh.elements=structmd.elements; end
[13717]	1037	if isfield(structmd,'edges'),
	1038	md.mesh.edges=structmd.edges;
	1039	md.mesh.edges(isnan(md.mesh.edges))=-1;
	1040	end
[13692]	1041	if isfield(structmd,'y'), md.mesh.y=structmd.y; end
	1042	if isfield(structmd,'x'), md.mesh.x=structmd.x; end
	1043	if isfield(structmd,'z'), md.mesh.z=structmd.z; end
[15771]	1044	if isfield(structmd,'diagnostic_ref'), md.stressbalance.referential=structmd.diagnostic_ref; end
[13692]	1045	if isfield(structmd,'npart'); md.qmu.numberofpartitions=structmd.npart; end
	1046	if isfield(structmd,'part'); md.qmu.partition=structmd.part; end
[13646]	1047
[13692]	1048	if isnumeric(md.verbose),
	1049	md.verbose=verbose;
	1050	end
[15768]	1051
[13692]	1052	if isfield(structmd,'spcvelocity'),
[15771]	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);
[13692]	1059	end
	1060	if isfield(structmd,'spcvx'),
[15771]	1061	md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
	1062	pos=find(~isnan(structmd.spcvx)); md.stressbalance.spcvx(pos)=structmd.spcvx(pos);
[13692]	1063	end
	1064	if isfield(structmd,'spcvy'),
[15771]	1065	md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
	1066	pos=find(~isnan(structmd.spcvy)); md.stressbalance.spcvy(pos)=structmd.spcvy(pos);
[13692]	1067	end
	1068	if isfield(structmd,'spcvz'),
[15771]	1069	md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
	1070	pos=find(~isnan(structmd.spcvz)); md.stressbalance.spcvz(pos)=structmd.spcvz(pos);
[13692]	1071	end
[14620]	1072	if isfield(structmd,'pressureload'),
	1073	if ~isempty(structmd.pressureload) & ismember(structmd.pressureload(end,end),[118 119 120]),
[15771]	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;
[14620]	1077	end
[13692]	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
[17941]	1106	if (structmd.groundingline_migration==273), md.groundingline.migration='AggressiveMigration'; end
[13692]	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
[11659]	1125
[13692]	1126	if isfield(structmd,'artificial_diffusivity') & structmd.artificial_diffusivity==2,
	1127	md.thermal.stabilization=2;
[15767]	1128	md.masstransport.stabilization=1;
[13692]	1129	md.balancethickness.stabilization=1;
	1130	end
[15767]	1131	if isnumeric(md.masstransport.hydrostatic_adjustment)
	1132	if md.masstransport.hydrostatic_adjustment==269,
	1133	md.masstransport.hydrostatic_adjustment='Incremental';
[13692]	1134	else
[15767]	1135	md.masstransport.hydrostatic_adjustment='Absolute';
[13692]	1136	end
	1137	end
[8952]	1138
[13692]	1139	%New fields
[19124]	1140	if ~isfield(structmd,'upperelements') & isa(md.mesh,'mesh3dprisms')
[13692]	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
[19124]	1144	if ~isfield(structmd,'lowerelements') & isa(md.mesh,'mesh3dprisms')
[13692]	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');
[15771]	1149	md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
[13692]	1150	end
[14529]	1151	if ~isfield(structmd,'loadingforce');
[15771]	1152	md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
[14529]	1153	end
[15768]	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
[15771]	1160	%2013 August 9
[15775]	1161	if isfield(structmd,'diagnostic') & (isa(structmd.diagnostic,'diagnostic') \|\| isa(structmd.diagnostic,'stressbalance')),
[15771]	1162	disp('Recovering old diagnostic class');
[15775]	1163	md.stressbalance=stressbalance(structmd.diagnostic);
[15771]	1164	end
[19642]	1165	%2014 January 9th
[21148]	1166	if isfield(structmd,'surfaceforcings') & isa(md.smb,'surfaceforcings'),
[19642]	1167	disp('Recovering old surfaceforcings class');
	1168	mass_balance=structmd.surfaceforcings.mass_balance;
	1169	md.smb=SMB();
	1170	md.smb.mass_balance=mass_balance;
	1171	end
	1172	%2015 September 10
[21148]	1173	if isfield(structmd,'surfaceforcings') & isa(structmd.surfaceforcings,'SMB'),
[19642]	1174	disp('Recovering old SMB class');
	1175	md.smb=SMBforcing(structmd.surfaceforcings);
	1176	end
[21148]	1177	if isfield(structmd,'surfaceforcings') & isa(structmd.surfaceforcings,'SMBhenning'),
	1178	disp('Recovering old SMBhenning class');
	1179	md.smb=SMBhenning(structmd.surfaceforcings);
	1180	end
[13692]	1181	end% }}}
	1182	function md = setdefaultparameters(md) % {{{
[8926]	1183
[13692]	1184	%initialize subclasses
[17558]	1185	md.mesh = mesh2d();
[13692]	1186	md.mask = mask();
	1187	md.constants = constants();
	1188	md.geometry = geometry();
	1189	md.initialization = initialization();
[19527]	1190	md.smb = SMBforcing();
[13692]	1191	md.basalforcings = basalforcings();
	1192	md.friction = friction();
	1193	md.rifts = rifts();
[19984]	1194	md.slr = slr();
[24469]	1195	md.dsl = dsl();
[13692]	1196	md.timestepping = timestepping();
	1197	md.groundingline = groundingline();
	1198	md.materials = matice();
[16160]	1199	md.damage = damage();
[13692]	1200	md.flowequation = flowequation();
	1201	md.debug = debug();
[14558]	1202	md.verbose = verbose();
[22296]	1203	md.settings = issmsettings();
[17932]	1204	md.toolkits = toolkits();
[13692]	1205	md.cluster = generic();
	1206	md.balancethickness = balancethickness();
[17079]	1207	md.stressbalance = stressbalance();
[14555]	1208	md.hydrology = hydrologyshreve();
[17079]	1209	md.masstransport = masstransport();
[13692]	1210	md.thermal = thermal();
	1211	md.steadystate = steadystate();
	1212	md.transient = transient();
[22040]	1213	md.levelset = levelset();
[18757]	1214	md.calving = calving();
[23652]	1215	md.frontalforcings = frontalforcings();
[22955]	1216	md.gia = giaivins();
	1217	md.esa = esa();
[22040]	1218	md.love = fourierlove();
[21260]	1219	md.esa = esa();
[13692]	1220	md.autodiff = autodiff();
	1221	md.inversion = inversion();
	1222	md.qmu = qmu();
[22955]	1223	md.amr = amr();
[13692]	1224	md.radaroverlay = radaroverlay();
	1225	md.results = struct();
[16388]	1226	md.outputdefinition = outputdefinition();
[13692]	1227	md.miscellaneous = miscellaneous();
	1228	md.private = private();
	1229	end
	1230	%}}}
[17483]	1231	function md = tetras(md,varargin) % {{{
	1232	%TETRAS - split 3d prismatic mesh into 3 tetrahedrons
	1233	%
	1234	% Usage:
	1235	% md=tetra(md)
	1236
	1237	if ~isa(md.mesh,'mesh3dprisms')
	1238	error('mesh is not a 3d prismatic mesh');
	1239	end
	1240
	1241	%Initialize tetra mesh
	1242	md.mesh=mesh3dtetras(md.mesh);
	1243
[17754]	1244	%Subdivision from Philipp Furnstahl (http://studierstube.icg.tugraz.at/thesis/fuernstahl_thesis.pdf)
	1245	steiner = 0;
	1246	nbv = md.mesh.numberofvertices;
	1247	nbt = 3*md.mesh.numberofelements;
	1248	elements = zeros(nbt,4);
	1249	for i=1:md.mesh.numberofelements
	1250	v1=md.mesh.elements(i,1); v2=md.mesh.elements(i,2); v3=md.mesh.elements(i,3);
	1251	v4=md.mesh.elements(i,4); v5=md.mesh.elements(i,5); v6=md.mesh.elements(i,6);
	1252	if(min(v2,v4)<min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)<min(v2,v6)),
	1253	steiner = steiner+1; nbv = nbv+1; nbt = nbt+5; v7 = nbv;
	1254	md.mesh.x=[md.mesh.x; mean(md.mesh.x(md.mesh.elements(i,:)))];
	1255	md.mesh.y=[md.mesh.y; mean(md.mesh.y(md.mesh.elements(i,:)))];
	1256	md.mesh.z=[md.mesh.z; mean(md.mesh.z(md.mesh.elements(i,:)))];
	1257	elements(3*(i-1)+1,:) = [v1 v2 v3 v7];
	1258	elements(3*(i-1)+2,:) = [v1 v2 v4 v7];
	1259	elements(3*(i-1)+3,:) = [v2 v4 v5 v7];
	1260	elements(end+1,:) = [v2 v3 v5 v7];
	1261	elements(end+1,:) = [v3 v5 v6 v7];
	1262	elements(end+1,:) = [v1 v3 v6 v7];
	1263	elements(end+1,:) = [v1 v4 v6 v7];
	1264	elements(end+1,:) = [v4 v5 v6 v7];
	1265	elseif(min(v2,v4)<min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)>min(v2,v6)),
	1266	elements(3*(i-1)+1,:) = [v1 v2 v4 v6];
	1267	elements(3*(i-1)+2,:) = [v2 v4 v5 v6];
	1268	elements(3*(i-1)+3,:) = [v1 v2 v3 v6];
	1269	elseif(min(v2,v4)<min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)<min(v2,v6)),
	1270	elements(3*(i-1)+1,:) = [v1 v2 v3 v4];
	1271	elements(3*(i-1)+2,:) = [v2 v3 v4 v5];
	1272	elements(3*(i-1)+3,:) = [v3 v4 v5 v6];
	1273	elseif(min(v2,v4)<min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)>min(v2,v6)),
	1274	elements(3*(i-1)+1,:) = [v1 v2 v3 v4];
	1275	elements(3*(i-1)+2,:) = [v2 v4 v5 v6];
	1276	elements(3*(i-1)+3,:) = [v2 v3 v4 v6];
[18142]	1277	elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)<min(v2,v6)),
[17754]	1278	elements(3*(i-1)+1,:) = [v1 v4 v5 v6];
	1279	elements(3*(i-1)+2,:) = [v1 v2 v3 v5];
	1280	elements(3*(i-1)+3,:) = [v1 v3 v5 v6];
	1281	elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)>min(v2,v6)),
	1282	elements(3*(i-1)+1,:) = [v1 v4 v5 v6];
	1283	elements(3*(i-1)+2,:) = [v1 v2 v5 v6];
	1284	elements(3*(i-1)+3,:) = [v1 v2 v3 v6];
	1285	elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)<min(v2,v6)),
	1286	elements(3*(i-1)+1,:) = [v1 v3 v4 v5];
	1287	elements(3*(i-1)+2,:) = [v1 v2 v3 v5];
	1288	elements(3*(i-1)+3,:) = [v3 v4 v5 v6];
	1289	elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)<min(v2,v6)),
	1290	elements(3*(i-1)+1,:) = [v1 v5 v6 v4];
	1291	elements(3*(i-1)+2,:) = [v1 v2 v3 v5];
	1292	elements(3*(i-1)+3,:) = [v5 v6 v3 v1];
	1293	elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)>min(v2,v6)),
	1294	steiner = steiner+1; nbv = nbv+1; nbt = nbt+5; v7 = nbv;
	1295	md.mesh.x=[md.mesh.x; mean(md.mesh.x(md.mesh.elements(i,:)))];
	1296	md.mesh.y=[md.mesh.y; mean(md.mesh.y(md.mesh.elements(i,:)))];
	1297	md.mesh.z=[md.mesh.z; mean(md.mesh.z(md.mesh.elements(i,:)))];
	1298	elements(3*(i-1)+1,:) = [v1 v2 v3 v7];
	1299	elements(3*(i-1)+2,:) = [v1 v4 v5 v7];
	1300	elements(3*(i-1)+3,:) = [v1 v2 v5 v7];
	1301	elements(end+1,:) = [v2 v5 v6 v7];
	1302	elements(end+1,:) = [v2 v3 v6 v7];
	1303	elements(end+1,:) = [v3 v4 v6 v7];
	1304	elements(end+1,:) = [v1 v3 v4 v7];
	1305	elements(end+1,:) = [v4 v5 v6 v7];
	1306	else
	1307	error('Case not supported'); %not supposed to happen!
	1308	end
[17774]	1309	%Reorder elements to make sure they are direct
	1310	for j=1:3
	1311	element = elements(3*(i-1)+j,:);
	1312	matrix = [md.mesh.x(element), md.mesh.y(element), md.mesh.z(element), ones(4,1)];
	1313	if det(matrix)>0,
	1314	elements(3*(i-1)+j,1)=element(2);
	1315	elements(3*(i-1)+j,2)=element(1);
	1316	end
	1317	end
[17754]	1318	end
	1319	%%Split in 3 tetras
	1320	%subelement1 = [1 2 3 5];
	1321	%subelement2 = [4 6 5 1];
	1322	%subelement3 = [5 6 3 1];
	1323	%elements=[md.mesh.elements(:,subelement1);md.mesh.elements(:,subelement2);md.mesh.elements(:,subelement3)];
[17774]	1324	if steiner==0,
	1325	disp('No Steiner point required to split prismatic mesh into tets');
	1326	else
	1327	disp([num2str(steiner) ' Steiner points had to be included'])
	1328	error('Steiner point not supported yet');
	1329	end
[17754]	1330
[17483]	1331	pos_elements = repmat([1:md.mesh.numberofelements]',3,1);
	1332
	1333	md.mesh.elements=elements;
	1334	md.mesh.numberofelements=size(elements,1);
	1335
	1336	%p and q (same deal, except for element that are on the bedrock: )
[17774]	1337	if ~isnan(md.friction.p),
	1338	md.friction.p=md.friction.p(pos_elements);
	1339	md.friction.q=md.friction.q(pos_elements);
	1340	end
[17483]	1341
	1342	%elementstype
	1343	if ~isnan(md.flowequation.element_equation)
	1344	oldelements_type=md.flowequation.element_equation;
	1345	md.flowequation.element_equation=md.flowequation.element_equation(pos_elements);
	1346	end
	1347
	1348	%connectivity
	1349	md.mesh.elementconnectivity=NaN;
	1350
	1351	%materials
[17774]	1352	if ~isnan(md.materials.rheology_n),
	1353	md.materials.rheology_n=md.materials.rheology_n(pos_elements);
	1354	end
[17483]	1355
	1356	%increase connectivity if less than 25:
	1357	if md.mesh.average_vertex_connectivity<=25,
	1358	md.mesh.average_vertex_connectivity=100;
	1359	end
	1360	end % }}}
[19040]	1361	function disp(self) % {{{
	1362	disp(sprintf('%19s: %-22s -- %s','mesh' ,['[1x1 ' class(self.mesh) ']'],'mesh properties'));
	1363	disp(sprintf('%19s: %-22s -- %s','mask' ,['[1x1 ' class(self.mask) ']'],'defines grounded and floating elements'));
	1364	disp(sprintf('%19s: %-22s -- %s','geometry' ,['[1x1 ' class(self.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...'));
	1365	disp(sprintf('%19s: %-22s -- %s','constants' ,['[1x1 ' class(self.constants) ']'],'physical constants'));
[19527]	1366	disp(sprintf('%19s: %-22s -- %s','smb' ,['[1x1 ' class(self.smb) ']'],'surface mass balance'));
[19040]	1367	disp(sprintf('%19s: %-22s -- %s','basalforcings' ,['[1x1 ' class(self.basalforcings) ']'],'bed forcings'));
	1368	disp(sprintf('%19s: %-22s -- %s','materials' ,['[1x1 ' class(self.materials) ']'],'material properties'));
	1369	disp(sprintf('%19s: %-22s -- %s','damage' ,['[1x1 ' class(self.damage) ']'],'parameters for damage evolution solution'));
	1370	disp(sprintf('%19s: %-22s -- %s','friction' ,['[1x1 ' class(self.friction) ']'],'basal friction/drag properties'));
	1371	disp(sprintf('%19s: %-22s -- %s','flowequation' ,['[1x1 ' class(self.flowequation) ']'],'flow equations'));
	1372	disp(sprintf('%19s: %-22s -- %s','timestepping' ,['[1x1 ' class(self.timestepping) ']'],'time stepping for transient models'));
	1373	disp(sprintf('%19s: %-22s -- %s','initialization' ,['[1x1 ' class(self.initialization) ']'],'initial guess/state'));
	1374	disp(sprintf('%19s: %-22s -- %s','rifts' ,['[1x1 ' class(self.rifts) ']'],'rifts properties'));
[19984]	1375	disp(sprintf('%19s: %-22s -- %s','slr' ,['[1x1 ' class(self.slr) ']'],'slr forcings'));
[24469]	1376	disp(sprintf('%19s: %-22s -- %s','dsl' ,['[1x1 ' class(self.dsl) ']'],'dynamic sea-level '));
[19040]	1377	disp(sprintf('%19s: %-22s -- %s','debug' ,['[1x1 ' class(self.debug) ']'],'debugging tools (valgrind, gprof)'));
	1378	disp(sprintf('%19s: %-22s -- %s','verbose' ,['[1x1 ' class(self.verbose) ']'],'verbosity level in solve'));
	1379	disp(sprintf('%19s: %-22s -- %s','settings' ,['[1x1 ' class(self.settings) ']'],'settings properties'));
	1380	disp(sprintf('%19s: %-22s -- %s','toolkits' ,['[1x1 ' class(self.toolkits) ']'],'PETSc options for each solution'));
	1381	disp(sprintf('%19s: %-22s -- %s','cluster' ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of cpus...)'));
	1382	disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(self.balancethickness) ']'],'parameters for balancethickness solution'));
	1383	disp(sprintf('%19s: %-22s -- %s','stressbalance' ,['[1x1 ' class(self.stressbalance) ']'],'parameters for stressbalance solution'));
	1384	disp(sprintf('%19s: %-22s -- %s','groundingline' ,['[1x1 ' class(self.groundingline) ']'],'parameters for groundingline solution'));
	1385	disp(sprintf('%19s: %-22s -- %s','hydrology' ,['[1x1 ' class(self.hydrology) ']'],'parameters for hydrology solution'));
	1386	disp(sprintf('%19s: %-22s -- %s','masstransport' ,['[1x1 ' class(self.masstransport) ']'],'parameters for masstransport solution'));
	1387	disp(sprintf('%19s: %-22s -- %s','thermal' ,['[1x1 ' class(self.thermal) ']'],'parameters for thermal solution'));
	1388	disp(sprintf('%19s: %-22s -- %s','steadystate' ,['[1x1 ' class(self.steadystate) ']'],'parameters for steadystate solution'));
	1389	disp(sprintf('%19s: %-22s -- %s','transient' ,['[1x1 ' class(self.transient) ']'],'parameters for transient solution'));
[21544]	1390	disp(sprintf('%19s: %-22s -- %s','levelset' ,['[1x1 ' class(self.levelset) ']'],'parameters for moving boundaries (level-set method)'));
[19040]	1391	disp(sprintf('%19s: %-22s -- %s','calving' ,['[1x1 ' class(self.calving) ']'],'parameters for calving'));
[23652]	1392	disp(sprintf('%19s: %-22s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings'));
[22955]	1393	disp(sprintf('%19s: %-22s -- %s','gia' ,['[1x1 ' class(self.gia) ']'],'parameters for gia solution'));
[21260]	1394	disp(sprintf('%19s: %-22s -- %s','esa' ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution'));
[22955]	1395	disp(sprintf('%19s: %-22s -- %s','love' ,['[1x1 ' class(self.love) ']'],'parameters for love solution'));
[19040]	1396	disp(sprintf('%19s: %-22s -- %s','autodiff' ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters'));
	1397	disp(sprintf('%19s: %-22s -- %s','inversion' ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods'));
	1398	disp(sprintf('%19s: %-22s -- %s','qmu' ,['[1x1 ' class(self.qmu) ']'],'dakota properties'));
[21674]	1399	disp(sprintf('%19s: %-22s -- %s','amr' ,['[1x1 ' class(self.amr) ']'],'adaptive mesh refinement properties'));
[19040]	1400	disp(sprintf('%19s: %-22s -- %s','outputdefinition',['[1x1 ' class(self.outputdefinition) ']'],'output definition'));
	1401	disp(sprintf('%19s: %-22s -- %s','results' ,['[1x1 ' class(self.results) ']'],'model results'));
	1402	disp(sprintf('%19s: %-22s -- %s','radaroverlay' ,['[1x1 ' class(self.radaroverlay) ']'],'radar image for plot overlay'));
	1403	disp(sprintf('%19s: %-22s -- %s','miscellaneous' ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields'));
[13692]	1404	end % }}}
[19040]	1405	function memory(self) % {{{
[15106]	1406
[14611]	1407	disp(sprintf('\nMemory imprint:\n'));
[14307]	1408
[14603]	1409	fields=properties('model');
	1410	mem=0;
[15106]	1411
[14603]	1412	for i=1:length(fields),
[19040]	1413	field=self.(fields{i});
[14603]	1414	s=whos('field');
	1415	mem=mem+s.bytes/1e6;
[14611]	1416	disp(sprintf('%19s: %6.2f Mb',fields{i},s.bytes/1e6));
[14307]	1417	end
[14603]	1418	disp(sprintf('%19s--%10s','--------------','--------------'));
	1419	disp(sprintf('%19s: %g Mb','Total',mem));
[14307]	1420	end % }}}
[19040]	1421	function netcdf(self,filename) % {{{
[14603]	1422	%NETCDF - save model as netcdf
	1423	%
	1424	% Usage:
	1425	% netcdf(md,filename)
	1426	%
	1427	% Example:
	1428	% netcdf(md,'model.nc');
	1429
	1430	disp('Saving model as NetCDF');
	1431	%1. Create NetCDF file
	1432	ncid=netcdf.create(filename,'CLOBBER');
	1433	netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Conventions','CF-1.4');
[19040]	1434	netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Title',['ISSM model (' self.miscellaneous.name ')']);
[14603]	1435	netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Author',getenv('USER'));
	1436	netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Date',datestr(now));
	1437
[14611]	1438	%Preallocate variable id, needed to write variables in netcdf file
[14631]	1439	var_id=zeros(1000,1);%preallocate
[14609]	1440
	1441	for step=1:2,
	1442	counter=0;
[19040]	1443	[var_id,counter]=structtonc(ncid,'md',self,0,var_id,counter,step);
[14609]	1444	if step==1, netcdf.endDef(ncid); end
[14603]	1445	end
[14611]	1446
[14631]	1447	if counter>1000,
	1448	warning(['preallocation of var_id need to be updated from ' num2str(1000) ' to ' num2str(counter)]);
[14609]	1449	end
	1450
	1451	netcdf.close(ncid)
[14603]	1452	end % }}}
[19040]	1453	function xylim(self) % {{{
[14405]	1454
[19040]	1455	xlim([min(self.mesh.x) max(self.mesh.x)]);
	1456	ylim([min(self.mesh.y) max(self.mesh.y)])
[14405]	1457	end % }}}
[15316]	1458	function md=upload(md) % {{{
	1459	%the goal of this routine is to upload the model onto a server, and to empty it.
	1460	%So first, save the model with a unique name and upload the file to the server:
	1461	random_part=fix(rand(1)*10000);
	1462	id=[md.miscellaneous.name '-' regexprep(datestr(now),'[^\w'']','') '-' num2str(random_part) '-' getenv('USER') '-' oshostname() '.upload'];
	1463	eval(['save ' id ' md']);
	1464
	1465	%Now, upload the file:
	1466	issmscpout(md.settings.upload_server,md.settings.upload_path,md.settings.upload_login,md.settings.upload_port,{id},1);
	1467
	1468	%Now, empty this model of everything except settings, and record name of file we just uploaded!
	1469	settings_back=md.settings;
	1470	md=model();
	1471	md.settings=settings_back;
	1472	md.settings.upload_filename=id;
	1473
	1474	%get locally rid of file that was uploaded
	1475	eval(['delete ' id]);
	1476
	1477	end % }}}
	1478	function md=download(md) % {{{
[15643]	1479
[15316]	1480	%the goal of this routine is to download the internals of the current model from a server, because
	1481	%this model is empty, except for the settings which tell us where to go and find this model!
[15643]	1482
[15316]	1483	%Download the file:
	1484	issmscpin(md.settings.upload_server, md.settings.upload_login, md.settings.upload_port, md.settings.upload_path, {md.settings.upload_filename});
	1485
	1486	name=md.settings.upload_filename;
	1487
	1488	%Now, load this model:
	1489	md=loadmodel(md.settings.upload_filename);
	1490
	1491	%get locally rid of file that was downloaded
	1492	eval(['delete ' name]);
	1493
	1494	end % }}}
[22955]	1495	function savemodeljs(md,modelname,websiteroot,varargin) % {{{
	1496
[19879]	1497	%the goal of this routine is to save the model as a javascript array that can be included in any html
	1498	%file:
	1499
[22955]	1500	options=pairoptions(varargin{:});
	1501	optimization=getfieldvalue(options,'optimize',0);
	1502
	1503
[19879]	1504	%disp:
[19894]	1505	disp(['saving model ''' modelname ''' in file ' websiteroot '/js/' modelname '.js']);
[19879]	1506
	1507	%open file for writing and declare the model:
[19894]	1508	fid=fopen([websiteroot '/js/' modelname '.js'],'w');
[19879]	1509	fprintf(fid,'var %s=new model();\n',modelname);
	1510
	1511	%now go through all the classes and fwrite all the corresponding fields:
	1512
	1513	fields=properties('model');
	1514	for i=1:length(fields),
	1515	field=fields{i};
	1516
	1517	%Some properties do not need to be saved
[19894]	1518	if ismember(field,{'results','cluster' }),
[19879]	1519	continue;
	1520	end
	1521
[22955]	1522	%some optimization:
	1523	if optimization==1,
	1524	%optimize for plotting only:
	1525	if ~ismember(field,{'geometry','mesh','mask'}),
	1526	continue;
	1527	end
	1528	end
	1529
[19879]	1530	%Check that current field is an object
	1531	if ~isobject(md.(field))
	1532	error(['field ''' char(field) ''' is not an object']);
	1533	end
	1534
	1535	%savemodeljs for current object
	1536	%disp(['javascript saving ' field '...']);
	1537	savemodeljs(md.(field),fid,modelname);
	1538	end
	1539
	1540	%done, close file:
	1541	fclose(fid);
	1542	end
[13692]	1543	end
[8926]	1544	end

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

Download in other formats: