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