Changeset 13692


Ignore:
Timestamp:
10/16/12 09:56:05 (12 years ago)
Author:
jschierm
Message:

CHG: Get rid of tab+space indent.

Location:
issm/trunk-jpl/src/m/classes/model
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/classes/model/model.m

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

  • issm/trunk-jpl/src/m/classes/model/model.py

    r13470 r13692  
    170170        # }}}
    171171
     172        """
     173        function md2 = extract(md,area) % {{{
     174                %extract - extract a model according to an Argus contour or flag list
     175                %
     176                %   This routine extracts a submodel from a bigger model with respect to a given contour
     177                %   md must be followed by the corresponding exp file or flags list
     178                %   It can either be a domain file (argus type, .exp extension), or an array of element flags.
     179                %   If user wants every element outside the domain to be
     180                %   extract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp');
     181                %   an empty string '' will be considered as an empty domain
     182                %   a string 'all' will be considered as the entire domain
     183                %   add an argument 0 if you do not want the elements to be checked (faster)
     184                %
     185                %   Usage:
     186                %      md2=extract(md,area);
     187                %
     188                %   Examples:
     189                %      md2=extract(md,'Domain.exp');
     190                %      md2=extract(md,md.mask.elementonfloatingice);
     191                %
     192                %   See also: EXTRUDE, COLLAPSE
     193
     194                %copy model
     195                md1=md;
     196
     197                %some checks
     198                if ((nargin~=2) | (nargout~=1)),
     199                        help extract
     200                        error('extract error message: bad usage');
     201                end
     202
     203                %get check option
     204                if (nargin==3 & varargin{1}==0),
     205                        checkoutline=0;
     206                else
     207                        checkoutline=1;
     208                end
     209
     210                %get elements that are inside area
     211                flag_elem=FlagElements(md1,area);
     212                if ~any(flag_elem),
     213                        error('extracted model is empty');
     214                end
     215
     216                %kick out all elements with 3 dirichlets
     217                spc_elem=find(~flag_elem);
     218                spc_node=sort(unique(md1.mesh.elements(spc_elem,:)));
     219                flag=ones(md1.mesh.numberofvertices,1);
     220                flag(spc_node)=0;
     221                pos=find(sum(flag(md1.mesh.elements),2)==0);
     222                flag_elem(pos)=0;
     223
     224                %extracted elements and nodes lists
     225                pos_elem=find(flag_elem);
     226                pos_node=sort(unique(md1.mesh.elements(pos_elem,:)));
     227
     228                %keep track of some fields
     229                numberofvertices1=md1.mesh.numberofvertices;
     230                numberofelements1=md1.mesh.numberofelements;
     231                numberofvertices2=length(pos_node);
     232                numberofelements2=length(pos_elem);
     233                flag_node=zeros(numberofvertices1,1);
     234                flag_node(pos_node)=1;
     235
     236                %Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
     237                Pelem=zeros(numberofelements1,1);
     238                Pelem(pos_elem)=[1:numberofelements2]';
     239                Pnode=zeros(numberofvertices1,1);
     240                Pnode(pos_node)=[1:numberofvertices2]';
     241
     242                %renumber the elements (some node won't exist anymore)
     243                elements_1=md1.mesh.elements;
     244                elements_2=elements_1(pos_elem,:);
     245                elements_2(:,1)=Pnode(elements_2(:,1));
     246                elements_2(:,2)=Pnode(elements_2(:,2));
     247                elements_2(:,3)=Pnode(elements_2(:,3));
     248                if md1.mesh.dimension==3,
     249                        elements_2(:,4)=Pnode(elements_2(:,4));
     250                        elements_2(:,5)=Pnode(elements_2(:,5));
     251                        elements_2(:,6)=Pnode(elements_2(:,6));
     252                end
     253
     254                %OK, now create the new model !
     255
     256                %take every fields from model
     257                md2=md1;
     258
     259                %automatically modify fields
     260
     261                %loop over model fields
     262                model_fields=fields(md1);
     263                for i=1:length(model_fields),
     264                        %get field
     265                        field=md1.(model_fields{i});
     266                        fieldsize=size(field);
     267                        if isobject(field), %recursive call
     268                                object_fields=fields(md1.(model_fields{i}));
     269                                for j=1:length(object_fields),
     270                                        %get field
     271                                        field=md1.(model_fields{i}).(object_fields{j});
     272                                        fieldsize=size(field);
     273                                        %size = number of nodes * n
     274                                        if fieldsize(1)==numberofvertices1
     275                                                md2.(model_fields{i}).(object_fields{j})=field(pos_node,:);
     276                                        elseif (fieldsize(1)==numberofvertices1+1)
     277                                                md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)];
     278                                                %size = number of elements * n
     279                                        elseif fieldsize(1)==numberofelements1
     280                                                md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:);
     281                                        end
     282                                end
     283                        else
     284                                %size = number of nodes * n
     285                                if fieldsize(1)==numberofvertices1
     286                                        md2.(model_fields{i})=field(pos_node,:);
     287                                elseif (fieldsize(1)==numberofvertices1+1)
     288                                        md2.(model_fields{i})=[field(pos_node,:); field(end,:)];
     289                                        %size = number of elements * n
     290                                elseif fieldsize(1)==numberofelements1
     291                                        md2.(model_fields{i})=field(pos_elem,:);
     292                                end
     293                        end
     294                end
     295
     296                %modify some specific fields
     297
     298                %Mesh
     299                md2.mesh.numberofelements=numberofelements2;
     300                md2.mesh.numberofvertices=numberofvertices2;
     301                md2.mesh.elements=elements_2;
     302
     303                %mesh.uppervertex mesh.lowervertex
     304                if md1.mesh.dimension==3
     305                        md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node);
     306                        pos=find(~isnan(md2.mesh.uppervertex));
     307                        md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos));
     308
     309                        md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node);
     310                        pos=find(~isnan(md2.mesh.lowervertex));
     311                        md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos));
     312
     313                        md2.mesh.upperelements=md1.mesh.upperelements(pos_elem);
     314                        pos=find(~isnan(md2.mesh.upperelements));
     315                        md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos));
     316
     317                        md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem);
     318                        pos=find(~isnan(md2.mesh.lowerelements));
     319                        md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos));
     320                end
     321
     322                %Initial 2d mesh
     323                if md1.mesh.dimension==3
     324                        flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d);
     325                        pos_elem_2d=find(flag_elem_2d);
     326                        flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d);
     327                        pos_node_2d=find(flag_node_2d);
     328
     329                        md2.mesh.numberofelements2d=length(pos_elem_2d);
     330                        md2.mesh.numberofvertices2d=length(pos_node_2d);
     331                        md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:);
     332                        md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1));
     333                        md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2));
     334                        md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3));
     335
     336                        md2.mesh.x2d=md1.mesh.x(pos_node_2d);
     337                        md2.mesh.y2d=md1.mesh.y(pos_node_2d);
     338                end
     339
     340                %Edges
     341                if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs...
     342                        %renumber first two columns
     343                        pos=find(md2.mesh.edges(:,4)~=-1);
     344                        md2.mesh.edges(:  ,1)=Pnode(md2.mesh.edges(:,1));
     345                        md2.mesh.edges(:  ,2)=Pnode(md2.mesh.edges(:,2));
     346                        md2.mesh.edges(:  ,3)=Pelem(md2.mesh.edges(:,3));
     347                        md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4));
     348                        %remove edges when the 2 vertices are not in the domain.
     349                        md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:);
     350                        %Replace all zeros by -1 in the last two columns;
     351                        pos=find(md2.mesh.edges(:,3)==0);
     352                        md2.mesh.edges(pos,3)=-1;
     353                        pos=find(md2.mesh.edges(:,4)==0);
     354                        md2.mesh.edges(pos,4)=-1;
     355                        %Invert -1 on the third column with last column (Also invert first two columns!!)
     356                        pos=find(md2.mesh.edges(:,3)==-1);
     357                        md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4);
     358                        md2.mesh.edges(pos,4)=-1;
     359                        values=md2.mesh.edges(pos,2);
     360                        md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1);
     361                        md2.mesh.edges(pos,1)=values;
     362                        %Finally remove edges that do not belong to any element
     363                        pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1);
     364                        md2.mesh.edges(pos,:)=[];
     365                end
     366
     367                %Penalties
     368                if ~isnan(md2.diagnostic.vertex_pairing),
     369                        for i=1:size(md1.diagnostic.vertex_pairing,1);
     370                                md2.diagnostic.vertex_pairing(i,:)=Pnode(md1.diagnostic.vertex_pairing(i,:));
     371                        end
     372                        md2.diagnostic.vertex_pairing=md2.diagnostic.vertex_pairing(find(md2.diagnostic.vertex_pairing(:,1)),:);
     373                end
     374                if ~isnan(md2.prognostic.vertex_pairing),
     375                        for i=1:size(md1.prognostic.vertex_pairing,1);
     376                                md2.prognostic.vertex_pairing(i,:)=Pnode(md1.prognostic.vertex_pairing(i,:));
     377                        end
     378                        md2.prognostic.vertex_pairing=md2.prognostic.vertex_pairing(find(md2.prognostic.vertex_pairing(:,1)),:);
     379                end
     380
     381                %recreate segments
     382                if md1.mesh.dimension==2
     383                        md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices);
     384                        md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
     385                        md2.mesh.segments=contourenvelope(md2);
     386                        md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
     387                else
     388                        %First do the connectivity for the contourenvelope in 2d
     389                        md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d);
     390                        md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity);
     391                        md2.mesh.segments=contourenvelope(md2);
     392                        md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
     393                        md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1);
     394                        %Then do it for 3d as usual
     395                        md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices);
     396                        md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
     397                end
     398
     399                %Boundary conditions: Dirichlets on new boundary
     400                %Catch the elements that have not been extracted
     401                orphans_elem=find(~flag_elem);
     402                orphans_node=unique(md1.mesh.elements(orphans_elem,:))';
     403                %Figure out which node are on the boundary between md2 and md1
     404                nodestoflag1=intersect(orphans_node,pos_node);
     405                nodestoflag2=Pnode(nodestoflag1);
     406                if numel(md1.diagnostic.spcvx)>1 & numel(md1.diagnostic.spcvy)>2 & numel(md1.diagnostic.spcvz)>2,
     407                        if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1
     408                                md2.diagnostic.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2);
     409                                md2.diagnostic.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2);
     410                        else
     411                                md2.diagnostic.spcvx(nodestoflag2)=NaN;
     412                                md2.diagnostic.spcvy(nodestoflag2)=NaN;
     413                                disp(' ')
     414                                disp('!! extract warning: spc values should be checked !!')
     415                                disp(' ')
     416                        end
     417                        %put 0 for vz
     418                        md2.diagnostic.spcvz(nodestoflag2)=0;
     419                end
     420                if ~isnan(md1.thermal.spctemperature),
     421                        md2.thermal.spctemperature(nodestoflag2,1)=1;
     422                end
     423
     424                %Diagnostic
     425                if ~isnan(md2.diagnostic.icefront)
     426                        md2.diagnostic.icefront(:,1)=Pnode(md1.diagnostic.icefront(:,1));
     427                        md2.diagnostic.icefront(:,2)=Pnode(md1.diagnostic.icefront(:,2));
     428                        md2.diagnostic.icefront(:,end-1)=Pelem(md1.diagnostic.icefront(:,end-1));
     429                        if md1.mesh.dimension==3
     430                                md2.diagnostic.icefront(:,3)=Pnode(md1.diagnostic.icefront(:,3));
     431                                md2.diagnostic.icefront(:,4)=Pnode(md1.diagnostic.icefront(:,4));
     432                        end
     433                        md2.diagnostic.icefront=md2.diagnostic.icefront(find(md2.diagnostic.icefront(:,1) & md2.diagnostic.icefront(:,2) & md2.diagnostic.icefront(:,end)),:);
     434                end
     435
     436                %Results fields
     437                if isstruct(md1.results),
     438                        md2.results=struct();
     439                        solutionfields=fields(md1.results);
     440                        for i=1:length(solutionfields),
     441                                %get subfields
     442                                solutionsubfields=fields(md1.results.(solutionfields{i}));
     443                                for j=1:length(solutionsubfields),
     444                                        field=md1.results.(solutionfields{i}).(solutionsubfields{j});
     445                                        if length(field)==numberofvertices1,
     446                                                md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_node);
     447                                        elseif length(field)==numberofelements1,
     448                                                md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_elem);
     449                                        else
     450                                                md2.results.(solutionfields{i}).(solutionsubfields{j})=field;
     451                                        end
     452                                end
     453                        end
     454                end
     455
     456                %Keep track of pos_node and pos_elem
     457                md2.mesh.extractedvertices=pos_node;
     458                md2.mesh.extractedelements=pos_elem;
     459        end % }}}
     460        """
     461
    172462        def extrude(md,*args):    # {{{
    173463                """
Note: See TracChangeset for help on using the changeset viewer.