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