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