Index: /issm/trunk-jpl/src/m/classes/model/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model/model.m (revision 13691) +++ /issm/trunk-jpl/src/m/classes/model/model.m (revision 13692) @@ -5,1119 +5,1119 @@ classdef model - properties (SetAccess=public) %Model fields - % {{{ - %Careful here: no other class should be used as default value this is a bug of matlab - mesh = 0; - mask = 0; - - geometry = 0; - constants = 0; - surfaceforcings = 0; - basalforcings = 0; - materials = 0; - friction = 0; - flowequation = 0; - timestepping = 0; - initialization = 0; - rifts = 0; - - debug = 0; - verbose = 0; - settings = 0; - solver = 0; - cluster = 0; - - balancethickness = 0; - diagnostic = 0; - groundingline = 0; - hydrology = 0; - prognostic = 0; - thermal = 0; - steadystate = 0; - transient = 0; - - autodiff = 0; - flaim = 0; - inversion = 0; - qmu = 0; - - results = 0; - radaroverlay = 0; - miscellaneous = 0; - private = 0; - - %}}} - end - methods (Static) - function md = loadobj(md) % {{{ - % This function is directly called by matlab when a model object is - % loaded. If the input is a struct it is an old version of model and - % old fields must be recovered (make sure they are in the deprecated - % model properties) - - if verLessThan('matlab','7.9'), - disp('Warning: your matlab version is old and there is a risk that load does not work correctly'); - disp(' if the model is not loaded correctly, rename temporarily loadobj so that matlab does not use it'); - - % This is a Matlab bug: all the fields of md have their default value - % Example of error message: - % Warning: Error loading an object of class 'model': - % Undefined function or method 'exist' for input arguments of type 'cell' - % - % This has been fixed in MATLAB 7.9 (R2009b) and later versions - end - - if isstruct(md) - disp('Recovering model object from a previous version'); - md = structtomodel(model,md); - end - - %2012 August 4th - if isa(md.materials,'materials'), - disp('Recovering old materials'); - if numel(md.materials.rheology_Z)==1 & isnan(md.materials.rheology_Z), - md.materials=matice(md.materials); - else - md.materials=matdamageice(md.materials); - end - end - - end% }}} - end - methods - function md = model(varargin) % {{{ - - switch nargin - case 0 - md=setdefaultparameters(md); - otherwise - error('model constructor error message: 0 of 1 argument only in input.'); - end - end - %}}} - function md = checkmessage(md,string) % {{{ - if(nargout~=1) error('wrong usage, model must be an output'); end - disp(['model not consistent: ' string]); - md.private.isconsistent=false; - end - %}}} - function md = collapse(md)% {{{ - %COLLAPSE - collapses a 3d mesh into a 2d mesh - % - % This routine collapses a 3d model into a 2d model - % and collapses all the fileds of the 3d model by - % taking their depth-averaged values - % - % Usage: - % md=collapse(md) - % - % See also: EXTRUDE, MODELEXTRACT - - %Check that the model is really a 3d model - if ~md.mesh.dimension==3, - error('collapse error message: only 3d mesh can be collapsed') - end - - %Start with changing alle the fields from the 3d mesh - - %drag is limited to nodes that are on the bedrock. - md.friction.coefficient=project2d(md,md.friction.coefficient,1); - - %p and q (same deal, except for element that are on the bedrock: ) - md.friction.p=project2d(md,md.friction.p,1); - md.friction.q=project2d(md,md.friction.q,1); - - %observations - if ~isnan(md.inversion.vx_obs), md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); end; - if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); end; - if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end; - if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end; - if numel(md.inversion.min_parameters)>1, md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end; - if numel(md.inversion.max_parameters)>1, md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end; - if ~isnan(md.surfaceforcings.mass_balance), - md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers); - end; - if ~isnan(md.balancethickness.thickening_rate), md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.mesh.numberoflayers); end; - - %results - if ~isnan(md.initialization.vx),md.initialization.vx=DepthAverage(md,md.initialization.vx);end; - if ~isnan(md.initialization.vy),md.initialization.vy=DepthAverage(md,md.initialization.vy);end; - if ~isnan(md.initialization.vz),md.initialization.vz=DepthAverage(md,md.initialization.vz);end; - if ~isnan(md.initialization.vel),md.initialization.vel=DepthAverage(md,md.initialization.vel);end; - if ~isnan(md.initialization.temperature),md.initialization.temperature=DepthAverage(md,md.initialization.temperature);end; - - %bedinfo and surface info - md.mesh.elementonbed=ones(md.mesh.numberofelements2d,1); - md.mesh.elementonsurface=ones(md.mesh.numberofelements2d,1); - md.mesh.vertexonbed=ones(md.mesh.numberofvertices2d,1); - md.mesh.vertexonsurface=ones(md.mesh.numberofvertices2d,1); - - %elementstype - if ~isnan(md.flowequation.element_equation) - md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1); - md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1); - md.flowequation.bordermacayeal=project2d(md,md.flowequation.bordermacayeal,1); - md.flowequation.borderpattyn=project2d(md,md.flowequation.borderpattyn,1); - md.flowequation.borderstokes=project2d(md,md.flowequation.borderstokes,1); - end - - %boundary conditions - md.diagnostic.spcvx=project2d(md,md.diagnostic.spcvx,md.mesh.numberoflayers); - md.diagnostic.spcvy=project2d(md,md.diagnostic.spcvy,md.mesh.numberoflayers); - md.diagnostic.spcvz=project2d(md,md.diagnostic.spcvz,md.mesh.numberoflayers); - md.diagnostic.referential=project2d(md,md.diagnostic.referential,md.mesh.numberoflayers); - md.prognostic.spcthickness=project2d(md,md.prognostic.spcthickness,md.mesh.numberoflayers); - md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers); - - %Extrusion of Neumann BC - if ~isnan(md.diagnostic.icefront), - numberofneumann2d=size(md.diagnostic.icefront,1)/(md.mesh.numberoflayers-1); - md.diagnostic.icefront=[md.diagnostic.icefront(1:numberofneumann2d,1:2) md.diagnostic.icefront(1:numberofneumann2d,5:6)]; %Add two columns on the first layer - end - - %materials - md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B); - md.materials.rheology_n=project2d(md,md.materials.rheology_n,1); - if isa(md.materials,'matdamageice') - md.materials.rheology_Z=DepthAverage(md,md.materials.rheology_Z); - end - - %special for thermal modeling: - md.basalforcings.melting_rate=project2d(md,md.basalforcings.melting_rate,1); - md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); %bedrock only gets geothermal flux - - %update of connectivity matrix - md.mesh.average_vertex_connectivity=25; - - %Collapse the mesh - nodes2d=md.mesh.numberofvertices2d; - elements2d=md.mesh.numberofelements2d; - - %parameters - md.geometry.surface=project2d(md,md.geometry.surface,1); - md.geometry.thickness=project2d(md,md.geometry.thickness,1); - md.geometry.bed=project2d(md,md.geometry.bed,1); - md.geometry.bathymetry=project2d(md,md.geometry.bathymetry,1); - md.mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1); - md.mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1); - md.mask.elementonfloatingice=project2d(md,md.mask.elementonfloatingice,1); - md.mask.vertexonfloatingice=project2d(md,md.mask.vertexonfloatingice,1); - md.mask.elementongroundedice=project2d(md,md.mask.elementongroundedice,1); - md.mask.vertexongroundedice=project2d(md,md.mask.vertexongroundedice,1); - md.mask.elementonwater=project2d(md,md.mask.elementonwater,1); - md.mask.vertexonwater=project2d(md,md.mask.vertexonwater,1); - - %lat long - if numel(md.mesh.lat) ==md.mesh.numberofvertices, md.mesh.lat=project2d(md,md.mesh.lat,1); end - if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end - - %Initialize with the 2d mesh - md.mesh.x=md.mesh.x2d; - md.mesh.y=md.mesh.y2d; - md.mesh.z=zeros(size(md.mesh.x2d)); - md.mesh.numberofvertices=md.mesh.numberofvertices2d; - md.mesh.numberofelements=md.mesh.numberofelements2d; - md.mesh.elements=md.mesh.elements2d; - - %Keep a trace of lower and upper nodes - md.mesh.lowervertex=NaN; - md.mesh.uppervertex=NaN; - md.mesh.lowerelements=NaN; - md.mesh.upperelements=NaN; - - %Remove old mesh - md.mesh.x2d=NaN; - md.mesh.y2d=NaN; - md.mesh.elements2d=NaN; - md.mesh.numberofelements2d=md.mesh.numberofelements; - md.mesh.numberofvertices2d=md.mesh.numberofvertices; - md.mesh.numberoflayers=0; - - %Update mesh type - md.mesh.dimension=2; - end % }}} - function md2 = extract(md,area) % {{{ - %extract - extract a model according to an Argus contour or flag list - % - % This routine extracts a submodel from a bigger model with respect to a given contour - % md must be followed by the corresponding exp file or flags list - % It can either be a domain file (argus type, .exp extension), or an array of element flags. - % If user wants every element outside the domain to be - % extract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp'); - % an empty string '' will be considered as an empty domain - % a string 'all' will be considered as the entire domain - % add an argument 0 if you do not want the elements to be checked (faster) - % - % Usage: - % md2=extract(md,area); - % - % Examples: - % md2=extract(md,'Domain.exp'); - % md2=extract(md,md.mask.elementonfloatingice); - % - % See also: EXTRUDE, COLLAPSE - - %copy model - md1=md; - - %some checks - if ((nargin~=2) | (nargout~=1)), - help extract - error('extract error message: bad usage'); - end - - %get check option - if (nargin==3 & varargin{1}==0), - checkoutline=0; - else - checkoutline=1; - end - - %get elements that are inside area - flag_elem=FlagElements(md1,area); - if ~any(flag_elem), - error('extracted model is empty'); - end - - %kick out all elements with 3 dirichlets - spc_elem=find(~flag_elem); - spc_node=sort(unique(md1.mesh.elements(spc_elem,:))); - flag=ones(md1.mesh.numberofvertices,1); - flag(spc_node)=0; - pos=find(sum(flag(md1.mesh.elements),2)==0); - flag_elem(pos)=0; - - %extracted elements and nodes lists - pos_elem=find(flag_elem); - pos_node=sort(unique(md1.mesh.elements(pos_elem,:))); - - %keep track of some fields - numberofvertices1=md1.mesh.numberofvertices; - numberofelements1=md1.mesh.numberofelements; - numberofvertices2=length(pos_node); - numberofelements2=length(pos_elem); - flag_node=zeros(numberofvertices1,1); - flag_node(pos_node)=1; - - %Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) - Pelem=zeros(numberofelements1,1); - Pelem(pos_elem)=[1:numberofelements2]'; - Pnode=zeros(numberofvertices1,1); - Pnode(pos_node)=[1:numberofvertices2]'; - - %renumber the elements (some node won't exist anymore) - elements_1=md1.mesh.elements; - elements_2=elements_1(pos_elem,:); - elements_2(:,1)=Pnode(elements_2(:,1)); - elements_2(:,2)=Pnode(elements_2(:,2)); - elements_2(:,3)=Pnode(elements_2(:,3)); - if md1.mesh.dimension==3, - elements_2(:,4)=Pnode(elements_2(:,4)); - elements_2(:,5)=Pnode(elements_2(:,5)); - elements_2(:,6)=Pnode(elements_2(:,6)); - end - - %OK, now create the new model ! - - %take every fields from model - md2=md1; - - %automatically modify fields - - %loop over model fields - model_fields=fields(md1); - for i=1:length(model_fields), - %get field - field=md1.(model_fields{i}); - fieldsize=size(field); - if isobject(field), %recursive call - object_fields=fields(md1.(model_fields{i})); - for j=1:length(object_fields), - %get field - field=md1.(model_fields{i}).(object_fields{j}); - fieldsize=size(field); - %size = number of nodes * n - if fieldsize(1)==numberofvertices1 - md2.(model_fields{i}).(object_fields{j})=field(pos_node,:); - elseif (fieldsize(1)==numberofvertices1+1) - md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)]; - %size = number of elements * n - elseif fieldsize(1)==numberofelements1 - md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:); - end - end - else - %size = number of nodes * n - if fieldsize(1)==numberofvertices1 - md2.(model_fields{i})=field(pos_node,:); - elseif (fieldsize(1)==numberofvertices1+1) - md2.(model_fields{i})=[field(pos_node,:); field(end,:)]; - %size = number of elements * n - elseif fieldsize(1)==numberofelements1 - md2.(model_fields{i})=field(pos_elem,:); - end - end - end - - %modify some specific fields - - %Mesh - md2.mesh.numberofelements=numberofelements2; - md2.mesh.numberofvertices=numberofvertices2; - md2.mesh.elements=elements_2; - - %mesh.uppervertex mesh.lowervertex - if md1.mesh.dimension==3 - md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node); - pos=find(~isnan(md2.mesh.uppervertex)); - md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos)); - - md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node); - pos=find(~isnan(md2.mesh.lowervertex)); - md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos)); - - md2.mesh.upperelements=md1.mesh.upperelements(pos_elem); - pos=find(~isnan(md2.mesh.upperelements)); - md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos)); - - md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem); - pos=find(~isnan(md2.mesh.lowerelements)); - md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos)); - end - - %Initial 2d mesh - if md1.mesh.dimension==3 - flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d); - pos_elem_2d=find(flag_elem_2d); - flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d); - pos_node_2d=find(flag_node_2d); - - md2.mesh.numberofelements2d=length(pos_elem_2d); - md2.mesh.numberofvertices2d=length(pos_node_2d); - md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:); - md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1)); - md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2)); - md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3)); - - md2.mesh.x2d=md1.mesh.x(pos_node_2d); - md2.mesh.y2d=md1.mesh.y(pos_node_2d); - end - - %Edges - if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs... - %renumber first two columns - pos=find(md2.mesh.edges(:,4)~=-1); - md2.mesh.edges(: ,1)=Pnode(md2.mesh.edges(:,1)); - md2.mesh.edges(: ,2)=Pnode(md2.mesh.edges(:,2)); - md2.mesh.edges(: ,3)=Pelem(md2.mesh.edges(:,3)); - md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4)); - %remove edges when the 2 vertices are not in the domain. - md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:); - %Replace all zeros by -1 in the last two columns; - pos=find(md2.mesh.edges(:,3)==0); - md2.mesh.edges(pos,3)=-1; - pos=find(md2.mesh.edges(:,4)==0); - md2.mesh.edges(pos,4)=-1; - %Invert -1 on the third column with last column (Also invert first two columns!!) - pos=find(md2.mesh.edges(:,3)==-1); - md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4); - md2.mesh.edges(pos,4)=-1; - values=md2.mesh.edges(pos,2); - md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1); - md2.mesh.edges(pos,1)=values; - %Finally remove edges that do not belong to any element - pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1); - md2.mesh.edges(pos,:)=[]; - end - - %Penalties - if ~isnan(md2.diagnostic.vertex_pairing), - for i=1:size(md1.diagnostic.vertex_pairing,1); - md2.diagnostic.vertex_pairing(i,:)=Pnode(md1.diagnostic.vertex_pairing(i,:)); - end - md2.diagnostic.vertex_pairing=md2.diagnostic.vertex_pairing(find(md2.diagnostic.vertex_pairing(:,1)),:); - end - if ~isnan(md2.prognostic.vertex_pairing), - for i=1:size(md1.prognostic.vertex_pairing,1); - md2.prognostic.vertex_pairing(i,:)=Pnode(md1.prognostic.vertex_pairing(i,:)); - end - md2.prognostic.vertex_pairing=md2.prognostic.vertex_pairing(find(md2.prognostic.vertex_pairing(:,1)),:); - end - - %recreate segments - if md1.mesh.dimension==2 - md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); - md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); - md2.mesh.segments=contourenvelope(md2); - md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; - else - %First do the connectivity for the contourenvelope in 2d - md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d); - md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity); - md2.mesh.segments=contourenvelope(md2); - md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; - md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1); - %Then do it for 3d as usual - md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); - md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); - end - - %Boundary conditions: Dirichlets on new boundary - %Catch the elements that have not been extracted - orphans_elem=find(~flag_elem); - orphans_node=unique(md1.mesh.elements(orphans_elem,:))'; - %Figure out which node are on the boundary between md2 and md1 - nodestoflag1=intersect(orphans_node,pos_node); - nodestoflag2=Pnode(nodestoflag1); - if numel(md1.diagnostic.spcvx)>1 & numel(md1.diagnostic.spcvy)>2 & numel(md1.diagnostic.spcvz)>2, - if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1 - md2.diagnostic.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2); - md2.diagnostic.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2); - else - md2.diagnostic.spcvx(nodestoflag2)=NaN; - md2.diagnostic.spcvy(nodestoflag2)=NaN; - disp(' ') - disp('!! extract warning: spc values should be checked !!') - disp(' ') - end - %put 0 for vz - md2.diagnostic.spcvz(nodestoflag2)=0; - end - if ~isnan(md1.thermal.spctemperature), - md2.thermal.spctemperature(nodestoflag2,1)=1; - end - - %Diagnostic - if ~isnan(md2.diagnostic.icefront) - md2.diagnostic.icefront(:,1)=Pnode(md1.diagnostic.icefront(:,1)); - md2.diagnostic.icefront(:,2)=Pnode(md1.diagnostic.icefront(:,2)); - md2.diagnostic.icefront(:,end-1)=Pelem(md1.diagnostic.icefront(:,end-1)); - if md1.mesh.dimension==3 - md2.diagnostic.icefront(:,3)=Pnode(md1.diagnostic.icefront(:,3)); - md2.diagnostic.icefront(:,4)=Pnode(md1.diagnostic.icefront(:,4)); - end - md2.diagnostic.icefront=md2.diagnostic.icefront(find(md2.diagnostic.icefront(:,1) & md2.diagnostic.icefront(:,2) & md2.diagnostic.icefront(:,end)),:); - end - - %Results fields - if isstruct(md1.results), - md2.results=struct(); - solutionfields=fields(md1.results); - for i=1:length(solutionfields), - %get subfields - solutionsubfields=fields(md1.results.(solutionfields{i})); - for j=1:length(solutionsubfields), - field=md1.results.(solutionfields{i}).(solutionsubfields{j}); - if length(field)==numberofvertices1, - md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_node); - elseif length(field)==numberofelements1, - md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_elem); - else - md2.results.(solutionfields{i}).(solutionsubfields{j})=field; - end - end - end - end - - %Keep track of pos_node and pos_elem - md2.mesh.extractedvertices=pos_node; - md2.mesh.extractedelements=pos_elem; - end % }}} - function md = extrude(md,varargin) % {{{ - %EXTRUDE - vertically extrude a 2d mesh - % - % vertically extrude a 2d mesh and create corresponding 3d mesh. - % The vertical distribution can: - % - follow a polynomial law - % - follow two polynomial laws, one for the lower part and one for the upper part of the mesh - % - be discribed by a list of coefficients (between 0 and 1) - % - % - % Usage: - % md=extrude(md,numlayers,extrusionexponent); - % md=extrude(md,numlayers,lowerexponent,upperexponent); - % md=extrude(md,listofcoefficients); - % - % Example: - % md=extrude(md,8,3); - % md=extrude(md,8,3,2); - % md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]); - % - % See also: MODELEXTRACT, COLLAPSE - - %some checks on list of arguments - if ((nargin>4) | (nargin<2) | (nargout~=1)), - help extrude; - error('extrude error message'); - end - - %Extrude the mesh - if nargin==2, %list of coefficients - clist=varargin{1}; - if any(clist<0) | any(clist>1), - error('extrusioncoefficients must be between 0 and 1'); - end - extrusionlist=sort(unique([clist(:);0;1])); - numlayers=length(extrusionlist); - elseif nargin==3, %one polynomial law - if varargin{2}<=0, - help extrude; - error('extrusionexponent must be >=0'); - end - numlayers=varargin{1}; - extrusionlist=((0:1:numlayers-1)/(numlayers-1)).^varargin{2}; - elseif nargin==4, %two polynomial laws - numlayers=varargin{1}; - lowerexp=varargin{2}; - upperexp=varargin{3}; - - if varargin{2}<=0 | varargin{3}<=0, - help extrude; - error('lower and upper extrusionexponents must be >=0'); - end - - lowerextrusionlist=[(0:2/(numlayers-1):1).^lowerexp]/2; - upperextrusionlist=[(0:2/(numlayers-1):1).^upperexp]/2; - extrusionlist=sort(unique([lowerextrusionlist 1-upperextrusionlist])); - - end - - if numlayers<2, - error('number of layers should be at least 2'); - end - if md.mesh.dimension==3, - error('Cannot extrude a 3d mesh (extrude cannot be called more than once)'); - end - - %Initialize with the 2d mesh - x3d=[]; - y3d=[]; - z3d=[]; %the lower node is on the bed - thickness3d=md.geometry.thickness; %thickness and bed for these nodes - bed3d=md.geometry.bed; - - %Create the new layers - for i=1:numlayers, - x3d=[x3d; md.mesh.x]; - y3d=[y3d; md.mesh.y]; - %nodes are distributed between bed and surface accordingly to the given exponent - z3d=[z3d; bed3d+thickness3d*extrusionlist(i)]; - end - number_nodes3d=size(x3d,1); %number of 3d nodes for the non extruded part of the mesh - - %Extrude elements - elements3d=[]; - for i=1:numlayers-1, - 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 - end - number_el3d=size(elements3d,1); %number of 3d nodes for the non extruded part of the mesh - - %Keep a trace of lower and upper nodes - mesh.lowervertex=NaN*ones(number_nodes3d,1); - mesh.uppervertex=NaN*ones(number_nodes3d,1); - mesh.lowervertex(md.mesh.numberofvertices+1:end)=1:(numlayers-1)*md.mesh.numberofvertices; - mesh.uppervertex(1:(numlayers-1)*md.mesh.numberofvertices)=md.mesh.numberofvertices+1:number_nodes3d; - md.mesh.lowervertex=mesh.lowervertex; - md.mesh.uppervertex=mesh.uppervertex; - - %same for lower and upper elements - mesh.lowerelements=NaN*ones(number_el3d,1); - mesh.upperelements=NaN*ones(number_el3d,1); - mesh.lowerelements(md.mesh.numberofelements+1:end)=1:(numlayers-2)*md.mesh.numberofelements; - mesh.upperelements(1:(numlayers-2)*md.mesh.numberofelements)=md.mesh.numberofelements+1:(numlayers-1)*md.mesh.numberofelements; - md.mesh.lowerelements=mesh.lowerelements; - md.mesh.upperelements=mesh.upperelements; - - %Save old mesh - md.mesh.x2d=md.mesh.x; - md.mesh.y2d=md.mesh.y; - md.mesh.elements2d=md.mesh.elements; - md.mesh.numberofelements2d=md.mesh.numberofelements; - md.mesh.numberofvertices2d=md.mesh.numberofvertices; - - %Update mesh type - md.mesh.dimension=3; - - %Build global 3d mesh - md.mesh.elements=elements3d; - md.mesh.x=x3d; - md.mesh.y=y3d; - md.mesh.z=z3d; - md.mesh.numberofelements=number_el3d; - md.mesh.numberofvertices=number_nodes3d; - md.mesh.numberoflayers=numlayers; - - %Ok, now deal with the other fields from the 2d mesh: - - %lat long - md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node'); - md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node'); - - %drag coefficient is limited to nodes that are on the bedrock. - md.friction.coefficient=project3d(md,'vector',md.friction.coefficient,'type','node','layer',1); - - %p and q (same deal, except for element that are on the bedrock: ) - md.friction.p=project3d(md,'vector',md.friction.p,'type','element'); - md.friction.q=project3d(md,'vector',md.friction.q,'type','element'); - - %observations - md.inversion.vx_obs=project3d(md,'vector',md.inversion.vx_obs,'type','node'); - md.inversion.vy_obs=project3d(md,'vector',md.inversion.vy_obs,'type','node'); - md.inversion.vel_obs=project3d(md,'vector',md.inversion.vel_obs,'type','node'); - md.surfaceforcings.mass_balance=project3d(md,'vector',md.surfaceforcings.mass_balance,'type','node'); - md.surfaceforcings.precipitation=project3d(md,'vector',md.surfaceforcings.precipitation,'type','node'); - md.balancethickness.thickening_rate=project3d(md,'vector',md.balancethickness.thickening_rate,'type','node'); - md.surfaceforcings.monthlytemperatures=project3d(md,'vector',md.surfaceforcings.monthlytemperatures,'type','node'); - - %results - if ~isnan(md.initialization.vx),md.initialization.vx=project3d(md,'vector',md.initialization.vx,'type','node');end; - if ~isnan(md.initialization.vy),md.initialization.vy=project3d(md,'vector',md.initialization.vy,'type','node');end; - if ~isnan(md.initialization.vz),md.initialization.vz=project3d(md,'vector',md.initialization.vz,'type','node');end; - if ~isnan(md.initialization.vel),md.initialization.vel=project3d(md,'vector',md.initialization.vel,'type','node');end; - if ~isnan(md.initialization.temperature),md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node');end; - if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node');end; - - %bedinfo and surface info - md.mesh.elementonbed=project3d(md,'vector',ones(md.mesh.numberofelements2d,1),'type','element','layer',1); - md.mesh.elementonsurface=project3d(md,'vector',ones(md.mesh.numberofelements2d,1),'type','element','layer',md.mesh.numberoflayers-1); - md.mesh.vertexonbed=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1); - md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers); - - %elementstype - if ~isnan(md.flowequation.element_equation) - oldelements_type=md.flowequation.element_equation; - md.flowequation.element_equation=zeros(number_el3d,1); - md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element'); - end - - %verticestype - if ~isnan(md.flowequation.vertex_equation) - oldvertices_type=md.flowequation.vertex_equation; - md.flowequation.vertex_equation=zeros(number_nodes3d,1); - md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node'); - end - md.flowequation.bordermacayeal=project3d(md,'vector',md.flowequation.bordermacayeal,'type','node'); - md.flowequation.borderpattyn=project3d(md,'vector',md.flowequation.borderpattyn,'type','node'); - md.flowequation.borderstokes=project3d(md,'vector',md.flowequation.borderstokes,'type','node'); - - %boundary conditions - md.diagnostic.spcvx=project3d(md,'vector',md.diagnostic.spcvx,'type','node'); - md.diagnostic.spcvy=project3d(md,'vector',md.diagnostic.spcvy,'type','node'); - md.diagnostic.spcvz=project3d(md,'vector',md.diagnostic.spcvz,'type','node'); - md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN); - md.prognostic.spcthickness=project3d(md,'vector',md.prognostic.spcthickness,'type','node'); - md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node'); - md.diagnostic.referential=project3d(md,'vector',md.diagnostic.referential,'type','node'); - - %in 3d, pressureload: [node1 node2 node3 node4 element] - 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 - pressureload=[]; - for i=1:numlayers-1, - pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d pressureload_layer1(:,6)]; - end - md.diagnostic.icefront=pressureload; - - %connectivity - md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1); - md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN; - for i=2:numlayers-1, - md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)... - =md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d; - end - md.mesh.elementconnectivity(find(isnan(md.mesh.elementconnectivity)))=0; - - %materials - md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node'); - md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element'); - if isa(md.materials,'matdamageice') - md.materials.rheology_Z=project3d(md,'vector',md.materials.rheology_Z,'type','node'); - end - - %parameters - md.geometry.surface=project3d(md,'vector',md.geometry.surface,'type','node'); - md.geometry.thickness=project3d(md,'vector',md.geometry.thickness,'type','node'); - md.geometry.hydrostatic_ratio=project3d(md,'vector',md.geometry.hydrostatic_ratio,'type','node'); - md.geometry.bed=project3d(md,'vector',md.geometry.bed,'type','node'); - md.geometry.bathymetry=project3d(md,'vector',md.geometry.bathymetry,'type','node'); - md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node'); - md.mask.elementonfloatingice=project3d(md,'vector',md.mask.elementonfloatingice,'type','element'); - md.mask.vertexonfloatingice=project3d(md,'vector',md.mask.vertexonfloatingice,'type','node'); - md.mask.elementongroundedice=project3d(md,'vector',md.mask.elementongroundedice,'type','element'); - md.mask.vertexongroundedice=project3d(md,'vector',md.mask.vertexongroundedice,'type','node'); - md.mask.elementonwater=project3d(md,'vector',md.mask.elementonwater,'type','element'); - md.mask.vertexonwater=project3d(md,'vector',md.mask.vertexonwater,'type','node'); - if ~isnan(md.inversion.cost_functions_coefficients),md.inversion.cost_functions_coefficients=project3d(md,'vector',md.inversion.cost_functions_coefficients,'type','node');end; - if ~isnan(md.inversion.min_parameters),md.inversion.min_parameters=project3d(md,'vector',md.inversion.min_parameters,'type','node');end; - if ~isnan(md.inversion.max_parameters),md.inversion.max_parameters=project3d(md,'vector',md.inversion.max_parameters,'type','node');end; - if ~isnan(md.qmu.partition),md.qmu.partition=project3d(md,'vector',md.qmu.partition','type','node');end - if(md.surfaceforcings.isdelta18o),md.surfaceforcings.temperatures_lgm=project3d(md,'vector',md.surfaceforcings.temperatures_lgm,'type','node');end - if(md.surfaceforcings.isdelta18o),md.surfaceforcings.temperatures_presentday=project3d(md,'vector',md.surfaceforcings.temperatures_presentday,'type','node');end - if(md.surfaceforcings.isdelta18o),md.surfaceforcings.precipitations_presentday=project3d(md,'vector',md.surfaceforcings.precipitations_presentday,'type','node');end - - %Put lithostatic pressure if there is an existing pressure - if ~isnan(md.initialization.pressure), - md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z); - end - - %special for thermal modeling: - md.basalforcings.melting_rate=project3d(md,'vector',md.basalforcings.melting_rate,'type','node','layer',1); - if ~isnan(md.basalforcings.geothermalflux) - md.basalforcings.geothermalflux=project3d(md,'vector',md.basalforcings.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux - end - - %increase connectivity if less than 25: - if md.mesh.average_vertex_connectivity<=25, - md.mesh.average_vertex_connectivity=100; - end + properties (SetAccess=public) %Model fields + % {{{ + %Careful here: no other class should be used as default value this is a bug of matlab + mesh = 0; + mask = 0; + + geometry = 0; + constants = 0; + surfaceforcings = 0; + basalforcings = 0; + materials = 0; + friction = 0; + flowequation = 0; + timestepping = 0; + initialization = 0; + rifts = 0; + + debug = 0; + verbose = 0; + settings = 0; + solver = 0; + cluster = 0; + + balancethickness = 0; + diagnostic = 0; + groundingline = 0; + hydrology = 0; + prognostic = 0; + thermal = 0; + steadystate = 0; + transient = 0; + + autodiff = 0; + flaim = 0; + inversion = 0; + qmu = 0; + + results = 0; + radaroverlay = 0; + miscellaneous = 0; + private = 0; + + %}}} + end + methods (Static) + function md = loadobj(md) % {{{ + % This function is directly called by matlab when a model object is + % loaded. If the input is a struct it is an old version of model and + % old fields must be recovered (make sure they are in the deprecated + % model properties) + + if verLessThan('matlab','7.9'), + disp('Warning: your matlab version is old and there is a risk that load does not work correctly'); + disp(' if the model is not loaded correctly, rename temporarily loadobj so that matlab does not use it'); + + % This is a Matlab bug: all the fields of md have their default value + % Example of error message: + % Warning: Error loading an object of class 'model': + % Undefined function or method 'exist' for input arguments of type 'cell' + % + % This has been fixed in MATLAB 7.9 (R2009b) and later versions + end + + if isstruct(md) + disp('Recovering model object from a previous version'); + md = structtomodel(model,md); + end + + %2012 August 4th + if isa(md.materials,'materials'), + disp('Recovering old materials'); + if numel(md.materials.rheology_Z)==1 & isnan(md.materials.rheology_Z), + md.materials=matice(md.materials); + else + md.materials=matdamageice(md.materials); + end + end + + end% }}} + end + methods + function md = model(varargin) % {{{ + + switch nargin + case 0 + md=setdefaultparameters(md); + otherwise + error('model constructor error message: 0 of 1 argument only in input.'); + end + end + %}}} + function md = checkmessage(md,string) % {{{ + if(nargout~=1) error('wrong usage, model must be an output'); end + disp(['model not consistent: ' string]); + md.private.isconsistent=false; + end + %}}} + function md = collapse(md)% {{{ + %COLLAPSE - collapses a 3d mesh into a 2d mesh + % + % This routine collapses a 3d model into a 2d model + % and collapses all the fileds of the 3d model by + % taking their depth-averaged values + % + % Usage: + % md=collapse(md) + % + % See also: EXTRUDE, MODELEXTRACT + + %Check that the model is really a 3d model + if ~md.mesh.dimension==3, + error('collapse error message: only 3d mesh can be collapsed') + end + + %Start with changing alle the fields from the 3d mesh + + %drag is limited to nodes that are on the bedrock. + md.friction.coefficient=project2d(md,md.friction.coefficient,1); + + %p and q (same deal, except for element that are on the bedrock: ) + md.friction.p=project2d(md,md.friction.p,1); + md.friction.q=project2d(md,md.friction.q,1); + + %observations + if ~isnan(md.inversion.vx_obs), md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); end; + if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); end; + if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end; + if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end; + if numel(md.inversion.min_parameters)>1, md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end; + if numel(md.inversion.max_parameters)>1, md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end; + if ~isnan(md.surfaceforcings.mass_balance), + md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers); + end; + if ~isnan(md.balancethickness.thickening_rate), md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.mesh.numberoflayers); end; + + %results + if ~isnan(md.initialization.vx),md.initialization.vx=DepthAverage(md,md.initialization.vx);end; + if ~isnan(md.initialization.vy),md.initialization.vy=DepthAverage(md,md.initialization.vy);end; + if ~isnan(md.initialization.vz),md.initialization.vz=DepthAverage(md,md.initialization.vz);end; + if ~isnan(md.initialization.vel),md.initialization.vel=DepthAverage(md,md.initialization.vel);end; + if ~isnan(md.initialization.temperature),md.initialization.temperature=DepthAverage(md,md.initialization.temperature);end; + + %bedinfo and surface info + md.mesh.elementonbed=ones(md.mesh.numberofelements2d,1); + md.mesh.elementonsurface=ones(md.mesh.numberofelements2d,1); + md.mesh.vertexonbed=ones(md.mesh.numberofvertices2d,1); + md.mesh.vertexonsurface=ones(md.mesh.numberofvertices2d,1); + + %elementstype + if ~isnan(md.flowequation.element_equation) + md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1); + md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1); + md.flowequation.bordermacayeal=project2d(md,md.flowequation.bordermacayeal,1); + md.flowequation.borderpattyn=project2d(md,md.flowequation.borderpattyn,1); + md.flowequation.borderstokes=project2d(md,md.flowequation.borderstokes,1); + end + + %boundary conditions + md.diagnostic.spcvx=project2d(md,md.diagnostic.spcvx,md.mesh.numberoflayers); + md.diagnostic.spcvy=project2d(md,md.diagnostic.spcvy,md.mesh.numberoflayers); + md.diagnostic.spcvz=project2d(md,md.diagnostic.spcvz,md.mesh.numberoflayers); + md.diagnostic.referential=project2d(md,md.diagnostic.referential,md.mesh.numberoflayers); + md.prognostic.spcthickness=project2d(md,md.prognostic.spcthickness,md.mesh.numberoflayers); + md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers); + + %Extrusion of Neumann BC + if ~isnan(md.diagnostic.icefront), + numberofneumann2d=size(md.diagnostic.icefront,1)/(md.mesh.numberoflayers-1); + md.diagnostic.icefront=[md.diagnostic.icefront(1:numberofneumann2d,1:2) md.diagnostic.icefront(1:numberofneumann2d,5:6)]; %Add two columns on the first layer + end + + %materials + md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B); + md.materials.rheology_n=project2d(md,md.materials.rheology_n,1); + if isa(md.materials,'matdamageice') + md.materials.rheology_Z=DepthAverage(md,md.materials.rheology_Z); + end + + %special for thermal modeling: + md.basalforcings.melting_rate=project2d(md,md.basalforcings.melting_rate,1); + md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); %bedrock only gets geothermal flux + + %update of connectivity matrix + md.mesh.average_vertex_connectivity=25; + + %Collapse the mesh + nodes2d=md.mesh.numberofvertices2d; + elements2d=md.mesh.numberofelements2d; + + %parameters + md.geometry.surface=project2d(md,md.geometry.surface,1); + md.geometry.thickness=project2d(md,md.geometry.thickness,1); + md.geometry.bed=project2d(md,md.geometry.bed,1); + md.geometry.bathymetry=project2d(md,md.geometry.bathymetry,1); + md.mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1); + md.mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1); + md.mask.elementonfloatingice=project2d(md,md.mask.elementonfloatingice,1); + md.mask.vertexonfloatingice=project2d(md,md.mask.vertexonfloatingice,1); + md.mask.elementongroundedice=project2d(md,md.mask.elementongroundedice,1); + md.mask.vertexongroundedice=project2d(md,md.mask.vertexongroundedice,1); + md.mask.elementonwater=project2d(md,md.mask.elementonwater,1); + md.mask.vertexonwater=project2d(md,md.mask.vertexonwater,1); + + %lat long + if numel(md.mesh.lat) ==md.mesh.numberofvertices, md.mesh.lat=project2d(md,md.mesh.lat,1); end + if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end + + %Initialize with the 2d mesh + md.mesh.x=md.mesh.x2d; + md.mesh.y=md.mesh.y2d; + md.mesh.z=zeros(size(md.mesh.x2d)); + md.mesh.numberofvertices=md.mesh.numberofvertices2d; + md.mesh.numberofelements=md.mesh.numberofelements2d; + md.mesh.elements=md.mesh.elements2d; + + %Keep a trace of lower and upper nodes + md.mesh.lowervertex=NaN; + md.mesh.uppervertex=NaN; + md.mesh.lowerelements=NaN; + md.mesh.upperelements=NaN; + + %Remove old mesh + md.mesh.x2d=NaN; + md.mesh.y2d=NaN; + md.mesh.elements2d=NaN; + md.mesh.numberofelements2d=md.mesh.numberofelements; + md.mesh.numberofvertices2d=md.mesh.numberofvertices; + md.mesh.numberoflayers=0; + + %Update mesh type + md.mesh.dimension=2; + end % }}} + function md2 = extract(md,area) % {{{ + %extract - extract a model according to an Argus contour or flag list + % + % This routine extracts a submodel from a bigger model with respect to a given contour + % md must be followed by the corresponding exp file or flags list + % It can either be a domain file (argus type, .exp extension), or an array of element flags. + % If user wants every element outside the domain to be + % extract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp'); + % an empty string '' will be considered as an empty domain + % a string 'all' will be considered as the entire domain + % add an argument 0 if you do not want the elements to be checked (faster) + % + % Usage: + % md2=extract(md,area); + % + % Examples: + % md2=extract(md,'Domain.exp'); + % md2=extract(md,md.mask.elementonfloatingice); + % + % See also: EXTRUDE, COLLAPSE + + %copy model + md1=md; + + %some checks + if ((nargin~=2) | (nargout~=1)), + help extract + error('extract error message: bad usage'); + end + + %get check option + if (nargin==3 & varargin{1}==0), + checkoutline=0; + else + checkoutline=1; + end + + %get elements that are inside area + flag_elem=FlagElements(md1,area); + if ~any(flag_elem), + error('extracted model is empty'); + end + + %kick out all elements with 3 dirichlets + spc_elem=find(~flag_elem); + spc_node=sort(unique(md1.mesh.elements(spc_elem,:))); + flag=ones(md1.mesh.numberofvertices,1); + flag(spc_node)=0; + pos=find(sum(flag(md1.mesh.elements),2)==0); + flag_elem(pos)=0; + + %extracted elements and nodes lists + pos_elem=find(flag_elem); + pos_node=sort(unique(md1.mesh.elements(pos_elem,:))); + + %keep track of some fields + numberofvertices1=md1.mesh.numberofvertices; + numberofelements1=md1.mesh.numberofelements; + numberofvertices2=length(pos_node); + numberofelements2=length(pos_elem); + flag_node=zeros(numberofvertices1,1); + flag_node(pos_node)=1; + + %Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) + Pelem=zeros(numberofelements1,1); + Pelem(pos_elem)=[1:numberofelements2]'; + Pnode=zeros(numberofvertices1,1); + Pnode(pos_node)=[1:numberofvertices2]'; + + %renumber the elements (some node won't exist anymore) + elements_1=md1.mesh.elements; + elements_2=elements_1(pos_elem,:); + elements_2(:,1)=Pnode(elements_2(:,1)); + elements_2(:,2)=Pnode(elements_2(:,2)); + elements_2(:,3)=Pnode(elements_2(:,3)); + if md1.mesh.dimension==3, + elements_2(:,4)=Pnode(elements_2(:,4)); + elements_2(:,5)=Pnode(elements_2(:,5)); + elements_2(:,6)=Pnode(elements_2(:,6)); + end + + %OK, now create the new model ! + + %take every fields from model + md2=md1; + + %automatically modify fields + + %loop over model fields + model_fields=fields(md1); + for i=1:length(model_fields), + %get field + field=md1.(model_fields{i}); + fieldsize=size(field); + if isobject(field), %recursive call + object_fields=fields(md1.(model_fields{i})); + for j=1:length(object_fields), + %get field + field=md1.(model_fields{i}).(object_fields{j}); + fieldsize=size(field); + %size = number of nodes * n + if fieldsize(1)==numberofvertices1 + md2.(model_fields{i}).(object_fields{j})=field(pos_node,:); + elseif (fieldsize(1)==numberofvertices1+1) + md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)]; + %size = number of elements * n + elseif fieldsize(1)==numberofelements1 + md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:); + end + end + else + %size = number of nodes * n + if fieldsize(1)==numberofvertices1 + md2.(model_fields{i})=field(pos_node,:); + elseif (fieldsize(1)==numberofvertices1+1) + md2.(model_fields{i})=[field(pos_node,:); field(end,:)]; + %size = number of elements * n + elseif fieldsize(1)==numberofelements1 + md2.(model_fields{i})=field(pos_elem,:); + end + end + end + + %modify some specific fields + + %Mesh + md2.mesh.numberofelements=numberofelements2; + md2.mesh.numberofvertices=numberofvertices2; + md2.mesh.elements=elements_2; + + %mesh.uppervertex mesh.lowervertex + if md1.mesh.dimension==3 + md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node); + pos=find(~isnan(md2.mesh.uppervertex)); + md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos)); + + md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node); + pos=find(~isnan(md2.mesh.lowervertex)); + md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos)); + + md2.mesh.upperelements=md1.mesh.upperelements(pos_elem); + pos=find(~isnan(md2.mesh.upperelements)); + md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos)); + + md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem); + pos=find(~isnan(md2.mesh.lowerelements)); + md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos)); + end + + %Initial 2d mesh + if md1.mesh.dimension==3 + flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d); + pos_elem_2d=find(flag_elem_2d); + flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d); + pos_node_2d=find(flag_node_2d); + + md2.mesh.numberofelements2d=length(pos_elem_2d); + md2.mesh.numberofvertices2d=length(pos_node_2d); + md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:); + md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1)); + md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2)); + md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3)); + + md2.mesh.x2d=md1.mesh.x(pos_node_2d); + md2.mesh.y2d=md1.mesh.y(pos_node_2d); + end + + %Edges + if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs... + %renumber first two columns + pos=find(md2.mesh.edges(:,4)~=-1); + md2.mesh.edges(: ,1)=Pnode(md2.mesh.edges(:,1)); + md2.mesh.edges(: ,2)=Pnode(md2.mesh.edges(:,2)); + md2.mesh.edges(: ,3)=Pelem(md2.mesh.edges(:,3)); + md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4)); + %remove edges when the 2 vertices are not in the domain. + md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:); + %Replace all zeros by -1 in the last two columns; + pos=find(md2.mesh.edges(:,3)==0); + md2.mesh.edges(pos,3)=-1; + pos=find(md2.mesh.edges(:,4)==0); + md2.mesh.edges(pos,4)=-1; + %Invert -1 on the third column with last column (Also invert first two columns!!) + pos=find(md2.mesh.edges(:,3)==-1); + md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4); + md2.mesh.edges(pos,4)=-1; + values=md2.mesh.edges(pos,2); + md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1); + md2.mesh.edges(pos,1)=values; + %Finally remove edges that do not belong to any element + pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1); + md2.mesh.edges(pos,:)=[]; + end + + %Penalties + if ~isnan(md2.diagnostic.vertex_pairing), + for i=1:size(md1.diagnostic.vertex_pairing,1); + md2.diagnostic.vertex_pairing(i,:)=Pnode(md1.diagnostic.vertex_pairing(i,:)); + end + md2.diagnostic.vertex_pairing=md2.diagnostic.vertex_pairing(find(md2.diagnostic.vertex_pairing(:,1)),:); + end + if ~isnan(md2.prognostic.vertex_pairing), + for i=1:size(md1.prognostic.vertex_pairing,1); + md2.prognostic.vertex_pairing(i,:)=Pnode(md1.prognostic.vertex_pairing(i,:)); + end + md2.prognostic.vertex_pairing=md2.prognostic.vertex_pairing(find(md2.prognostic.vertex_pairing(:,1)),:); + end + + %recreate segments + if md1.mesh.dimension==2 + md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); + md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); + md2.mesh.segments=contourenvelope(md2); + md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; + else + %First do the connectivity for the contourenvelope in 2d + md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d); + md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity); + md2.mesh.segments=contourenvelope(md2); + md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; + md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1); + %Then do it for 3d as usual + md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); + md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); + end + + %Boundary conditions: Dirichlets on new boundary + %Catch the elements that have not been extracted + orphans_elem=find(~flag_elem); + orphans_node=unique(md1.mesh.elements(orphans_elem,:))'; + %Figure out which node are on the boundary between md2 and md1 + nodestoflag1=intersect(orphans_node,pos_node); + nodestoflag2=Pnode(nodestoflag1); + if numel(md1.diagnostic.spcvx)>1 & numel(md1.diagnostic.spcvy)>2 & numel(md1.diagnostic.spcvz)>2, + if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1 + md2.diagnostic.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2); + md2.diagnostic.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2); + else + md2.diagnostic.spcvx(nodestoflag2)=NaN; + md2.diagnostic.spcvy(nodestoflag2)=NaN; + disp(' ') + disp('!! extract warning: spc values should be checked !!') + disp(' ') + end + %put 0 for vz + md2.diagnostic.spcvz(nodestoflag2)=0; + end + if ~isnan(md1.thermal.spctemperature), + md2.thermal.spctemperature(nodestoflag2,1)=1; + end + + %Diagnostic + if ~isnan(md2.diagnostic.icefront) + md2.diagnostic.icefront(:,1)=Pnode(md1.diagnostic.icefront(:,1)); + md2.diagnostic.icefront(:,2)=Pnode(md1.diagnostic.icefront(:,2)); + md2.diagnostic.icefront(:,end-1)=Pelem(md1.diagnostic.icefront(:,end-1)); + if md1.mesh.dimension==3 + md2.diagnostic.icefront(:,3)=Pnode(md1.diagnostic.icefront(:,3)); + md2.diagnostic.icefront(:,4)=Pnode(md1.diagnostic.icefront(:,4)); + end + md2.diagnostic.icefront=md2.diagnostic.icefront(find(md2.diagnostic.icefront(:,1) & md2.diagnostic.icefront(:,2) & md2.diagnostic.icefront(:,end)),:); + end + + %Results fields + if isstruct(md1.results), + md2.results=struct(); + solutionfields=fields(md1.results); + for i=1:length(solutionfields), + %get subfields + solutionsubfields=fields(md1.results.(solutionfields{i})); + for j=1:length(solutionsubfields), + field=md1.results.(solutionfields{i}).(solutionsubfields{j}); + if length(field)==numberofvertices1, + md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_node); + elseif length(field)==numberofelements1, + md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_elem); + else + md2.results.(solutionfields{i}).(solutionsubfields{j})=field; + end + end + end + end + + %Keep track of pos_node and pos_elem + md2.mesh.extractedvertices=pos_node; + md2.mesh.extractedelements=pos_elem; + end % }}} + function md = extrude(md,varargin) % {{{ + %EXTRUDE - vertically extrude a 2d mesh + % + % vertically extrude a 2d mesh and create corresponding 3d mesh. + % The vertical distribution can: + % - follow a polynomial law + % - follow two polynomial laws, one for the lower part and one for the upper part of the mesh + % - be discribed by a list of coefficients (between 0 and 1) + % + % + % Usage: + % md=extrude(md,numlayers,extrusionexponent); + % md=extrude(md,numlayers,lowerexponent,upperexponent); + % md=extrude(md,listofcoefficients); + % + % Example: + % md=extrude(md,8,3); + % md=extrude(md,8,3,2); + % md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]); + % + % See also: MODELEXTRACT, COLLAPSE + + %some checks on list of arguments + if ((nargin>4) | (nargin<2) | (nargout~=1)), + help extrude; + error('extrude error message'); + end + + %Extrude the mesh + if nargin==2, %list of coefficients + clist=varargin{1}; + if any(clist<0) | any(clist>1), + error('extrusioncoefficients must be between 0 and 1'); + end + extrusionlist=sort(unique([clist(:);0;1])); + numlayers=length(extrusionlist); + elseif nargin==3, %one polynomial law + if varargin{2}<=0, + help extrude; + error('extrusionexponent must be >=0'); + end + numlayers=varargin{1}; + extrusionlist=((0:1:numlayers-1)/(numlayers-1)).^varargin{2}; + elseif nargin==4, %two polynomial laws + numlayers=varargin{1}; + lowerexp=varargin{2}; + upperexp=varargin{3}; + + if varargin{2}<=0 | varargin{3}<=0, + help extrude; + error('lower and upper extrusionexponents must be >=0'); + end + + lowerextrusionlist=[(0:2/(numlayers-1):1).^lowerexp]/2; + upperextrusionlist=[(0:2/(numlayers-1):1).^upperexp]/2; + extrusionlist=sort(unique([lowerextrusionlist 1-upperextrusionlist])); + + end + + if numlayers<2, + error('number of layers should be at least 2'); + end + if md.mesh.dimension==3, + error('Cannot extrude a 3d mesh (extrude cannot be called more than once)'); + end + + %Initialize with the 2d mesh + x3d=[]; + y3d=[]; + z3d=[]; %the lower node is on the bed + thickness3d=md.geometry.thickness; %thickness and bed for these nodes + bed3d=md.geometry.bed; + + %Create the new layers + for i=1:numlayers, + x3d=[x3d; md.mesh.x]; + y3d=[y3d; md.mesh.y]; + %nodes are distributed between bed and surface accordingly to the given exponent + z3d=[z3d; bed3d+thickness3d*extrusionlist(i)]; + end + number_nodes3d=size(x3d,1); %number of 3d nodes for the non extruded part of the mesh + + %Extrude elements + elements3d=[]; + for i=1:numlayers-1, + 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 + end + number_el3d=size(elements3d,1); %number of 3d nodes for the non extruded part of the mesh + + %Keep a trace of lower and upper nodes + mesh.lowervertex=NaN*ones(number_nodes3d,1); + mesh.uppervertex=NaN*ones(number_nodes3d,1); + mesh.lowervertex(md.mesh.numberofvertices+1:end)=1:(numlayers-1)*md.mesh.numberofvertices; + mesh.uppervertex(1:(numlayers-1)*md.mesh.numberofvertices)=md.mesh.numberofvertices+1:number_nodes3d; + md.mesh.lowervertex=mesh.lowervertex; + md.mesh.uppervertex=mesh.uppervertex; + + %same for lower and upper elements + mesh.lowerelements=NaN*ones(number_el3d,1); + mesh.upperelements=NaN*ones(number_el3d,1); + mesh.lowerelements(md.mesh.numberofelements+1:end)=1:(numlayers-2)*md.mesh.numberofelements; + mesh.upperelements(1:(numlayers-2)*md.mesh.numberofelements)=md.mesh.numberofelements+1:(numlayers-1)*md.mesh.numberofelements; + md.mesh.lowerelements=mesh.lowerelements; + md.mesh.upperelements=mesh.upperelements; + + %Save old mesh + md.mesh.x2d=md.mesh.x; + md.mesh.y2d=md.mesh.y; + md.mesh.elements2d=md.mesh.elements; + md.mesh.numberofelements2d=md.mesh.numberofelements; + md.mesh.numberofvertices2d=md.mesh.numberofvertices; + + %Update mesh type + md.mesh.dimension=3; + + %Build global 3d mesh + md.mesh.elements=elements3d; + md.mesh.x=x3d; + md.mesh.y=y3d; + md.mesh.z=z3d; + md.mesh.numberofelements=number_el3d; + md.mesh.numberofvertices=number_nodes3d; + md.mesh.numberoflayers=numlayers; + + %Ok, now deal with the other fields from the 2d mesh: + + %lat long + md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node'); + md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node'); + + %drag coefficient is limited to nodes that are on the bedrock. + md.friction.coefficient=project3d(md,'vector',md.friction.coefficient,'type','node','layer',1); + + %p and q (same deal, except for element that are on the bedrock: ) + md.friction.p=project3d(md,'vector',md.friction.p,'type','element'); + md.friction.q=project3d(md,'vector',md.friction.q,'type','element'); + + %observations + md.inversion.vx_obs=project3d(md,'vector',md.inversion.vx_obs,'type','node'); + md.inversion.vy_obs=project3d(md,'vector',md.inversion.vy_obs,'type','node'); + md.inversion.vel_obs=project3d(md,'vector',md.inversion.vel_obs,'type','node'); + md.surfaceforcings.mass_balance=project3d(md,'vector',md.surfaceforcings.mass_balance,'type','node'); + md.surfaceforcings.precipitation=project3d(md,'vector',md.surfaceforcings.precipitation,'type','node'); + md.balancethickness.thickening_rate=project3d(md,'vector',md.balancethickness.thickening_rate,'type','node'); + md.surfaceforcings.monthlytemperatures=project3d(md,'vector',md.surfaceforcings.monthlytemperatures,'type','node'); + + %results + if ~isnan(md.initialization.vx),md.initialization.vx=project3d(md,'vector',md.initialization.vx,'type','node');end; + if ~isnan(md.initialization.vy),md.initialization.vy=project3d(md,'vector',md.initialization.vy,'type','node');end; + if ~isnan(md.initialization.vz),md.initialization.vz=project3d(md,'vector',md.initialization.vz,'type','node');end; + if ~isnan(md.initialization.vel),md.initialization.vel=project3d(md,'vector',md.initialization.vel,'type','node');end; + if ~isnan(md.initialization.temperature),md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node');end; + if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node');end; + + %bedinfo and surface info + md.mesh.elementonbed=project3d(md,'vector',ones(md.mesh.numberofelements2d,1),'type','element','layer',1); + md.mesh.elementonsurface=project3d(md,'vector',ones(md.mesh.numberofelements2d,1),'type','element','layer',md.mesh.numberoflayers-1); + md.mesh.vertexonbed=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1); + md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers); + + %elementstype + if ~isnan(md.flowequation.element_equation) + oldelements_type=md.flowequation.element_equation; + md.flowequation.element_equation=zeros(number_el3d,1); + md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element'); + end + + %verticestype + if ~isnan(md.flowequation.vertex_equation) + oldvertices_type=md.flowequation.vertex_equation; + md.flowequation.vertex_equation=zeros(number_nodes3d,1); + md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node'); + end + md.flowequation.bordermacayeal=project3d(md,'vector',md.flowequation.bordermacayeal,'type','node'); + md.flowequation.borderpattyn=project3d(md,'vector',md.flowequation.borderpattyn,'type','node'); + md.flowequation.borderstokes=project3d(md,'vector',md.flowequation.borderstokes,'type','node'); + + %boundary conditions + md.diagnostic.spcvx=project3d(md,'vector',md.diagnostic.spcvx,'type','node'); + md.diagnostic.spcvy=project3d(md,'vector',md.diagnostic.spcvy,'type','node'); + md.diagnostic.spcvz=project3d(md,'vector',md.diagnostic.spcvz,'type','node'); + md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN); + md.prognostic.spcthickness=project3d(md,'vector',md.prognostic.spcthickness,'type','node'); + md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node'); + md.diagnostic.referential=project3d(md,'vector',md.diagnostic.referential,'type','node'); + + %in 3d, pressureload: [node1 node2 node3 node4 element] + 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 + pressureload=[]; + for i=1:numlayers-1, + pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d pressureload_layer1(:,6)]; + end + md.diagnostic.icefront=pressureload; + + %connectivity + md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1); + md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN; + for i=2:numlayers-1, + md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)... + =md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d; + end + md.mesh.elementconnectivity(find(isnan(md.mesh.elementconnectivity)))=0; + + %materials + md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node'); + md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element'); + if isa(md.materials,'matdamageice') + md.materials.rheology_Z=project3d(md,'vector',md.materials.rheology_Z,'type','node'); + end + + %parameters + md.geometry.surface=project3d(md,'vector',md.geometry.surface,'type','node'); + md.geometry.thickness=project3d(md,'vector',md.geometry.thickness,'type','node'); + md.geometry.hydrostatic_ratio=project3d(md,'vector',md.geometry.hydrostatic_ratio,'type','node'); + md.geometry.bed=project3d(md,'vector',md.geometry.bed,'type','node'); + md.geometry.bathymetry=project3d(md,'vector',md.geometry.bathymetry,'type','node'); + md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node'); + md.mask.elementonfloatingice=project3d(md,'vector',md.mask.elementonfloatingice,'type','element'); + md.mask.vertexonfloatingice=project3d(md,'vector',md.mask.vertexonfloatingice,'type','node'); + md.mask.elementongroundedice=project3d(md,'vector',md.mask.elementongroundedice,'type','element'); + md.mask.vertexongroundedice=project3d(md,'vector',md.mask.vertexongroundedice,'type','node'); + md.mask.elementonwater=project3d(md,'vector',md.mask.elementonwater,'type','element'); + md.mask.vertexonwater=project3d(md,'vector',md.mask.vertexonwater,'type','node'); + if ~isnan(md.inversion.cost_functions_coefficients),md.inversion.cost_functions_coefficients=project3d(md,'vector',md.inversion.cost_functions_coefficients,'type','node');end; + if ~isnan(md.inversion.min_parameters),md.inversion.min_parameters=project3d(md,'vector',md.inversion.min_parameters,'type','node');end; + if ~isnan(md.inversion.max_parameters),md.inversion.max_parameters=project3d(md,'vector',md.inversion.max_parameters,'type','node');end; + if ~isnan(md.qmu.partition),md.qmu.partition=project3d(md,'vector',md.qmu.partition','type','node');end + if(md.surfaceforcings.isdelta18o),md.surfaceforcings.temperatures_lgm=project3d(md,'vector',md.surfaceforcings.temperatures_lgm,'type','node');end + if(md.surfaceforcings.isdelta18o),md.surfaceforcings.temperatures_presentday=project3d(md,'vector',md.surfaceforcings.temperatures_presentday,'type','node');end + if(md.surfaceforcings.isdelta18o),md.surfaceforcings.precipitations_presentday=project3d(md,'vector',md.surfaceforcings.precipitations_presentday,'type','node');end + + %Put lithostatic pressure if there is an existing pressure + if ~isnan(md.initialization.pressure), + md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z); + end + + %special for thermal modeling: + md.basalforcings.melting_rate=project3d(md,'vector',md.basalforcings.melting_rate,'type','node','layer',1); + if ~isnan(md.basalforcings.geothermalflux) + md.basalforcings.geothermalflux=project3d(md,'vector',md.basalforcings.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux + end + + %increase connectivity if less than 25: + if md.mesh.average_vertex_connectivity<=25, + md.mesh.average_vertex_connectivity=100; + end end % }}} - function md = structtomodel(md,structmd) % {{{ - - if ~isstruct(structmd) error('input model is not a structure'); end - - %loaded model is a struct, initialize output and recover all fields - md = structtoobj(model,structmd); - - %Old field now classes - if (isfield(structmd,'timestepping') & isnumeric(md.timestepping)), md.timestepping=timestepping(); end - if (isfield(structmd,'mask') & isnumeric(md.mask)),md.mask=mask(); end - - %Field name change - if isfield(structmd,'drag'), md.friction.coefficient=structmd.drag; end - if isfield(structmd,'p'), md.friction.p=structmd.p; end - if isfield(structmd,'q'), md.friction.q=structmd.p; end - if isfield(structmd,'melting'), md.basalforcings.melting_rate=structmd.melting; end - if isfield(structmd,'melting_rate'), md.basalforcings.melting_rate=structmd.melting_rate; end - if isfield(structmd,'accumulation'), md.surfaceforcings.mass_balance=structmd.accumulation; end - if isfield(structmd,'numberofgrids'), md.mesh.numberofvertices=structmd.numberofgrids; end - if isfield(structmd,'numberofgrids2d'), md.mesh.numberofvertices2d=structmd.numberofgrids2d; end - if isfield(structmd,'uppergrids'), md.mesh.uppervertex=structmd.uppergrids; end - if isfield(structmd,'lowergrids'), md.mesh.lowervertex=structmd.lowergrids; end - if isfield(structmd,'gridonbed'), md.mesh.vertexonbed=structmd.gridonbed; end - if isfield(structmd,'gridonsurface'), md.mesh.vertexonsurface=structmd.gridonsurface; end - if isfield(structmd,'extractedgrids'), md.mesh.extractedvertices=structmd.extractedgrids; end - if isfield(structmd,'gridoniceshelf'), md.mask.vertexonfloatingice=structmd.gridoniceshelf; end - if isfield(structmd,'gridonicesheet'), md.mask.vertexongroundedice=structmd.gridonicesheet; end - if isfield(structmd,'gridonwater'), md.mask.vertexonwater=structmd.gridonwater; end - if isfield(structmd,'gridonboundary'), md.mesh.vertexonboundary=structmd.gridonboundary; end - if isfield(structmd,'petscoptions') & ~isempty(structmd.petscoptions), md.solver=structmd.petscoptions; end - if isfield(structmd,'g'), md.constants.g=structmd.g; end - if isfield(structmd,'yts'), md.constants.yts=structmd.yts; end - if isfield(structmd,'surface_mass_balance'), md.surfaceforcings.mass_balance=structmd.surface_mass_balance; end - if isfield(structmd,'basal_melting_rate'), md.basalforcings.melting_rate=structmd.basal_melting_rate; end - if isfield(structmd,'basal_melting_rate_correction'), md.basalforcings.melting_rate_correction=structmd.basal_melting_rate_correction; end - if isfield(structmd,'geothermalflux'), md.basalforcings.geothermalflux=structmd.geothermalflux; end - if isfield(structmd,'drag'), md.friction.coefficient=structmd.drag; end - if isfield(structmd,'drag_coefficient'), md.friction.coefficient=structmd.drag_coefficient; end - if isfield(structmd,'drag_p'), md.friction.p=structmd.drag_p; end - if isfield(structmd,'drag_q'), md.friction.q=structmd.drag_q; end - if isfield(structmd,'riftproperties'), %old implementation - md.rifts=rifts(); - md.rifts.riftproperties=structmd.riftproperties; - md.rifts.riftstruct=structmd.rifts; - md.rifts.riftproperties=structmd.riftinfo; - end - if isfield(structmd,'bamg'), md.private.bamg=structmd.bamg; end - if isfield(structmd,'lowmem'), md.settings.lowmem=structmd.lowmem; end - if isfield(structmd,'io_gather'), md.settings.io_gather=structmd.io_gather; end - if isfield(structmd,'spcwatercolumn'), md.hydrology.spcwatercolumn=structmd.spcwatercolumn; end - if isfield(structmd,'hydro_n'), md.hydrology.n=structmd.hydro_n; end - if isfield(structmd,'hydro_p'), md.hydrology.p=structmd.hydro_p; end - if isfield(structmd,'hydro_q'), md.hydrology.q=structmd.hydro_q; end - if isfield(structmd,'hydro_CR'), md.hydrology.CR=structmd.hydro_CR; end - if isfield(structmd,'hydro_kn'), md.hydrology.kn=structmd.hydro_kn; end - if isfield(structmd,'spctemperature'), md.thermal.spctemperature=structmd.spctemperature; end - if isfield(structmd,'min_thermal_constraints'), md.thermal.penalty_threshold=structmd.min_thermal_constraints; end - if isfield(structmd,'artificial_diffusivity'), md.thermal.stabilization=structmd.artificial_diffusivity; end - if isfield(structmd,'max_nonlinear_iterations'), md.thermal.maxiter=structmd.max_nonlinear_iterations; end - if isfield(structmd,'stabilize_constraints'), md.thermal.penalty_lock=structmd.stabilize_constraints; end - if isfield(structmd,'penalty_offset'), md.thermal.penalty_factor=structmd.penalty_offset; end - if isfield(structmd,'name'), md.miscellaneous.name=structmd.name; end - if isfield(structmd,'notes'), md.miscellaneous.notes=structmd.notes; end - if isfield(structmd,'dummy'), md.miscellaneous.dummy=structmd.dummy; end - if isfield(structmd,'dt'), md.timestepping.time_step=structmd.dt; end - if isfield(structmd,'ndt'), md.timestepping.final_time=structmd.ndt; end - if isfield(structmd,'time_adapt'), md.timestepping.time_adapt=structmd.time_adapt; end - if isfield(structmd,'cfl_coefficient'), md.timestepping.cfl_coefficient=structmd.cfl_coefficient; end - if isfield(structmd,'spcthickness'), md.prognostic.spcthickness=structmd.spcthickness; end - if isfield(structmd,'artificial_diffusivity'), md.prognostic.stabilization=structmd.artificial_diffusivity; end - if isfield(structmd,'hydrostatic_adjustment'), md.prognostic.hydrostatic_adjustment=structmd.hydrostatic_adjustment; end - if isfield(structmd,'penalties'), md.prognostic.vertex_pairing=structmd.penalties; end - if isfield(structmd,'penalty_offset'), md.prognostic.penalty_factor=structmd.penalty_offset; end - if isfield(structmd,'B'), md.materials.rheology_B=structmd.B; end - if isfield(structmd,'n'), md.materials.rheology_n=structmd.n; end - if isfield(structmd,'rheology_B'), md.materials.rheology_B=structmd.rheology_B; end - if isfield(structmd,'rheology_n'), md.materials.rheology_n=structmd.rheology_n; end - if isfield(structmd,'elementoniceshelf'), md.mask.elementonfloatingice=structmd.elementoniceshelf; end - if isfield(structmd,'elementonicesheet'), md.mask.elementongroundedice=structmd.elementonicesheet; end - if isfield(structmd,'elementonwater'), md.mask.elementonwater=structmd.elementonwater; end - if isfield(structmd,'nodeoniceshelf'), md.mask.vertexonfloatingice=structmd.nodeoniceshelf; end - if isfield(structmd,'nodeonicesheet'), md.mask.vertexongroundedice=structmd.nodeonicesheet; end - if isfield(structmd,'nodeonwater'), md.mask.vertexonwater=structmd.nodeonwater; end - if isfield(structmd,'spcthickness'), md.balancethickness.spcthickness=structmd.spcthickness; end - if isfield(structmd,'artificial_diffusivity'), md.balancethickness.stabilization=structmd.artificial_diffusivity; end - if isfield(structmd,'dhdt'), md.balancethickness.thickening_rate=structmd.dhdt; end - if isfield(structmd,'ismacayealpattyn'), md.flowequation.ismacayealpattyn=structmd.ismacayealpattyn; end - if isfield(structmd,'ishutter'), md.flowequation.ishutter=structmd.ishutter; end - if isfield(structmd,'isstokes'), md.flowequation.isstokes=structmd.isstokes; end - if isfield(structmd,'elements_type'), md.flowequation.element_equation=structmd.elements_type; end - if isfield(structmd,'vertices_type'), md.flowequation.vertex_equation=structmd.vertices_type; end - if isfield(structmd,'eps_rel'), md.steadystate.reltol=structmd.eps_rel; end - if isfield(structmd,'max_steadystate_iterations'), md.steadystate.maxiter=structmd.max_steadystate_iterations; end - if isfield(structmd,'isdiagnostic'), md.transient.isdiagnostic=structmd.isdiagnostic; end - if isfield(structmd,'isprognostic'), md.transient.isprognostic=structmd.isprognostic; end - if isfield(structmd,'isthermal'), md.transient.isthermal=structmd.isthermal; end - if isfield(structmd,'control_analysis'), md.inversion.iscontrol=structmd.control_analysis; end - if isfield(structmd,'weights'), md.inversion.cost_functions_coefficients=structmd.weights; end - if isfield(structmd,'nsteps'), md.inversion.nsteps=structmd.nsteps; end - if isfield(structmd,'maxiter_per_step'), md.inversion.maxiter_per_step=structmd.maxiter_per_step; end - if isfield(structmd,'cm_min'), md.inversion.min_parameters=structmd.cm_min; end - if isfield(structmd,'cm_max'), md.inversion.max_parameters=structmd.cm_max; end - if isfield(structmd,'vx_obs'), md.inversion.vx_obs=structmd.vx_obs; end - if isfield(structmd,'vy_obs'), md.inversion.vy_obs=structmd.vy_obs; end - if isfield(structmd,'vel_obs'), md.inversion.vel_obs=structmd.vel_obs; end - if isfield(structmd,'thickness_obs'), md.inversion.thickness_obs=structmd.thickness_obs; end - if isfield(structmd,'vx'), md.initialization.vx=structmd.vx; end - if isfield(structmd,'vy'), md.initialization.vy=structmd.vy; end - if isfield(structmd,'vz'), md.initialization.vz=structmd.vz; end - if isfield(structmd,'vel'), md.initialization.vel=structmd.vel; end - if isfield(structmd,'pressure'), md.initialization.pressure=structmd.pressure; end - if isfield(structmd,'temperature'), md.initialization.temperature=structmd.temperature; end - if isfield(structmd,'waterfraction'), md.initialization.waterfraction=structmd.waterfraction; end - if isfield(structmd,'watercolumn'), md.initialization.watercolumn=structmd.watercolumn; end - if isfield(structmd,'surface'), md.geometry.surface=structmd.surface; end - if isfield(structmd,'bed'), md.geometry.bed=structmd.bed; end - if isfield(structmd,'thickness'), md.geometry.thickness=structmd.thickness; end - if isfield(structmd,'bathymetry'), md.geometry.bathymetry=structmd.bathymetry; end - if isfield(structmd,'thickness_coeff'), md.geometry.hydrostatic_ratio=structmd.thickness_coeff; end - if isfield(structmd,'connectivity'), md.mesh.average_vertex_connectivity=structmd.connectivity; end - if isfield(structmd,'extractednodes'), md.mesh.extractedvertices=structmd.extractednodes; end - if isfield(structmd,'extractedelements'), md.mesh.extractedelements=structmd.extractedelements; end - if isfield(structmd,'nodeonboundary'), md.mesh.vertexonboundary=structmd.nodeonboundary; end - if isfield(structmd,'hemisphere'), md.mesh.hemisphere=structmd.hemisphere; end - if isfield(structmd,'lat'), md.mesh.lat=structmd.lat; end - if isfield(structmd,'long'), md.mesh.long=structmd.long; end - if isfield(structmd,'segments'), md.mesh.segments=structmd.segments; end - if isfield(structmd,'segmentmarkers'), md.mesh.segmentmarkers=structmd.segmentmarkers; end - if isfield(structmd,'dim'), md.mesh.dimension=structmd.dim; end - if isfield(structmd,'numlayers'), md.mesh.numberoflayers=structmd.numlayers; end - if isfield(structmd,'numberofelements'), md.mesh.numberofelements=structmd.numberofelements; end - if isfield(structmd,'numberofvertices'), md.mesh.numberofvertices=structmd.numberofvertices; end - if isfield(structmd,'numberofnodes'), md.mesh.numberofvertices=structmd.numberofnodes; end - if isfield(structmd,'numberofedges'), md.mesh.numberofedges=structmd.numberofedges; end - if isfield(structmd,'numberofelements2d'), md.mesh.numberofelements2d=structmd.numberofelements2d; end - if isfield(structmd,'numberofnodes2d'), md.mesh.numberofvertices2d=structmd.numberofnodes2d; end - if isfield(structmd,'nodeconnectivity'), md.mesh.vertexconnectivity=structmd.nodeconnectivity; end - if isfield(structmd,'elementconnectivity'), md.mesh.elementconnectivity=structmd.elementconnectivity; end - if isfield(structmd,'uppernodes'), md.mesh.uppervertex=structmd.uppernodes; end - if isfield(structmd,'lowernodes'), md.mesh.lowervertex=structmd.lowernodes; end - if isfield(structmd,'upperelements'), md.mesh.upperelements=structmd.upperelements; end - if isfield(structmd,'lowerelements'), md.mesh.lowerelements=structmd.lowerelements; end - if isfield(structmd,'elementonbed'), md.mesh.elementonbed=structmd.elementonbed; end - if isfield(structmd,'elementonsurface'), md.mesh.elementonsurface=structmd.elementonsurface; end - if isfield(structmd,'nodeonsurface'), md.mesh.vertexonsurface=structmd.nodeonsurface; end - if isfield(structmd,'nodeonbed'), md.mesh.vertexonbed=structmd.nodeonbed; end - if isfield(structmd,'elements2d'), md.mesh.elements2d=structmd.elements2d; end - if isfield(structmd,'y2d'), md.mesh.y2d=structmd.y2d; end - if isfield(structmd,'x2d'), md.mesh.x2d=structmd.x2d; end - if isfield(structmd,'elements'), md.mesh.elements=structmd.elements; end - if isfield(structmd,'edges'), md.mesh.edges=structmd.edges; end - if isfield(structmd,'y'), md.mesh.y=structmd.y; end - if isfield(structmd,'x'), md.mesh.x=structmd.x; end - if isfield(structmd,'z'), md.mesh.z=structmd.z; end - if isfield(structmd,'mask'), md.flaim.criterion=structmd.mask; end - if isfield(structmd,'pressureload'), md.diagnostic.icefront=structmd.pressureload; end - if isfield(structmd,'diagnostic_ref'), md.diagnostic.referential=structmd.diagnostic_ref; end - if isfield(structmd,'npart'); md.qmu.numberofpartitions=structmd.npart; end - if isfield(structmd,'part'); md.qmu.partition=structmd.part; end - - %Field changes - if (isfield(structmd,'type') & ischar(structmd.type)), - if strcmpi(structmd.type,'2d'), md.mesh.dimension=2; end - if strcmpi(structmd.type,'3d'), md.mesh.dimension=3; end - end - if isnumeric(md.verbose), - md.verbose=verbose; - end - if size(md.diagnostic.icefront,2)==3 || size(md.diagnostic.icefront,2)==5, - front=md.diagnostic.icefront; - md.diagnostic.icefront=[front 1*md.mask.elementonfloatingice(front(:,end))]; - end - if isfield(structmd,'spcvelocity'), - md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1); - md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1); - md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1); - pos=find(structmd.spcvelocity(:,1)); md.diagnostic.spcvx(pos)=structmd.spcvelocity(pos,4); - pos=find(structmd.spcvelocity(:,2)); md.diagnostic.spcvy(pos)=structmd.spcvelocity(pos,5); - pos=find(structmd.spcvelocity(:,3)); md.diagnostic.spcvz(pos)=structmd.spcvelocity(pos,6); - end - if isfield(structmd,'spcvx'), - md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1); - pos=find(~isnan(structmd.spcvx)); md.diagnostic.spcvx(pos)=structmd.spcvx(pos); - end - if isfield(structmd,'spcvy'), - md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1); - pos=find(~isnan(structmd.spcvy)); md.diagnostic.spcvy(pos)=structmd.spcvy(pos); - end - if isfield(structmd,'spcvz'), - md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1); - pos=find(~isnan(structmd.spcvz)); md.diagnostic.spcvz(pos)=structmd.spcvz(pos); - end - if ~isempty(structmd.pressureload) & ismember(structmd.pressureload(end,end),[118 119 120]), - pos=find(structmd.pressureload(:,end)==120); md.diagnostic.icefront(pos,end)=0; - pos=find(structmd.pressureload(:,end)==118); md.diagnostic.icefront(pos,end)=1; - pos=find(structmd.pressureload(:,end)==119); md.diagnostic.icefront(pos,end)=2; - end - if isfield(structmd,'elements_type') & structmd.elements_type(end,end)>50, - pos=find(structmd.elements_type==59); md.flowequation.element_equation(pos,end)=0; - pos=find(structmd.elements_type==55); md.flowequation.element_equation(pos,end)=1; - pos=find(structmd.elements_type==56); md.flowequation.element_equation(pos,end)=2; - pos=find(structmd.elements_type==60); md.flowequation.element_equation(pos,end)=3; - pos=find(structmd.elements_type==62); md.flowequation.element_equation(pos,end)=4; - pos=find(structmd.elements_type==57); md.flowequation.element_equation(pos,end)=5; - pos=find(structmd.elements_type==58); md.flowequation.element_equation(pos,end)=6; - pos=find(structmd.elements_type==61); md.flowequation.element_equation(pos,end)=7; - end - if isfield(structmd,'vertices_type') & structmd.vertices_type(end,end)>50, - pos=find(structmd.vertices_type==59); md.flowequation.vertex_equation(pos,end)=0; - pos=find(structmd.vertices_type==55); md.flowequation.vertex_equation(pos,end)=1; - pos=find(structmd.vertices_type==56); md.flowequation.vertex_equation(pos,end)=2; - pos=find(structmd.vertices_type==60); md.flowequation.vertex_equation(pos,end)=3; - pos=find(structmd.vertices_type==62); md.flowequation.vertex_equation(pos,end)=4; - pos=find(structmd.vertices_type==57); md.flowequation.vertex_equation(pos,end)=5; - pos=find(structmd.vertices_type==58); md.flowequation.vertex_equation(pos,end)=6; - pos=find(structmd.vertices_type==61); md.flowequation.vertex_equation(pos,end)=7; - end - if isfield(structmd,'rheology_law') & isnumeric(structmd.rheology_law), - if (structmd.rheology_law==272), md.materials.rheology_law='None'; end - if (structmd.rheology_law==368), md.materials.rheology_law='Paterson'; end - if (structmd.rheology_law==369), md.materials.rheology_law='Arrhenius'; end - end - if isfield(structmd,'groundingline_migration') & isnumeric(structmd.groundingline_migration), - if (structmd.groundingline_migration==272), md.groundingline.migration='None'; end - if (structmd.groundingline_migration==273), md.groundingline.migration='AgressiveMigration'; end - if (structmd.groundingline_migration==274), md.groundingline.migration='SoftMigration'; end - end - if isfield(structmd,'control_type') & isnumeric(structmd.control_type), - if (structmd.control_type==143), md.inversion.control_parameters={'FrictionCoefficient'}; end - if (structmd.control_type==190), md.inversion.control_parameters={'RheologyBbar'}; end - if (structmd.control_type==147), md.inversion.control_parameters={'Thickeningrate'}; end - end - if isfield(structmd,'cm_responses') & ismember(structmd.cm_responses(end,end),[165:170 383 388 389]), - pos=find(structmd.cm_responses==166); md.inversion.cost_functions(pos)=101; - pos=find(structmd.cm_responses==167); md.inversion.cost_functions(pos)=102; - pos=find(structmd.cm_responses==168); md.inversion.cost_functions(pos)=103; - pos=find(structmd.cm_responses==169); md.inversion.cost_functions(pos)=104; - pos=find(structmd.cm_responses==170); md.inversion.cost_functions(pos)=105; - pos=find(structmd.cm_responses==165); md.inversion.cost_functions(pos)=201; - pos=find(structmd.cm_responses==389); md.inversion.cost_functions(pos)=501; - pos=find(structmd.cm_responses==388); md.inversion.cost_functions(pos)=502; - pos=find(structmd.cm_responses==382); md.inversion.cost_functions(pos)=503; - end - - if isfield(structmd,'artificial_diffusivity') & structmd.artificial_diffusivity==2, - md.thermal.stabilization=2; - md.prognostic.stabilization=1; - md.balancethickness.stabilization=1; - end - if isnumeric(md.prognostic.hydrostatic_adjustment) - if md.prognostic.hydrostatic_adjustment==269, - md.prognostic.hydrostatic_adjustment='Incremental'; - else - md.prognostic.hydrostatic_adjustment='Absolute'; - end - end - - %New fields - if ~isfield(structmd,'upperelements'); - md.mesh.upperelements=transpose(1:md.mesh.numberofelements)+md.mesh.numberofelements2d; - md.mesh.upperelements(end-md.mesh.numberofelements2d+1:end)=NaN; - end - if ~isfield(structmd,'lowerelements'); - md.mesh.lowerelements=transpose(1:md.mesh.numberofelements)-md.mesh.numberofelements2d; - md.mesh.lowerelements(1:md.mesh.numberofelements2d)=NaN; - end - - if ~isfield(structmd,'diagnostic_ref'); - md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6); - end - - end% }}} - function md = setdefaultparameters(md) % {{{ - - %initialize subclasses - md.mesh = mesh(); - md.mask = mask(); - md.constants = constants(); - md.geometry = geometry(); - md.initialization = initialization(); - md.surfaceforcings = surfaceforcings(); - md.basalforcings = basalforcings(); - md.friction = friction(); - md.rifts = rifts(); - md.timestepping = timestepping(); - md.groundingline = groundingline(); - md.materials = matice(); - md.flowequation = flowequation(); - md.debug = debug(); - md.verbose = verbose('solution',true,'qmu',true,'control',true); - md.settings = settings(); - md.solver = solver(); - if ismumps(), - md.solver = addoptions(md.solver,DiagnosticVertAnalysisEnum(),mumpsoptions()); - else - md.solver = addoptions(md.solver,DiagnosticVertAnalysisEnum(),iluasmoptions()); - end - md.cluster = generic(); - md.balancethickness = balancethickness(); - md.diagnostic = diagnostic(); - md.hydrology = hydrology(); - md.prognostic = prognostic(); - md.thermal = thermal(); - md.steadystate = steadystate(); - md.transient = transient(); - md.autodiff = autodiff(); - md.flaim = flaim(); - md.inversion = inversion(); - md.qmu = qmu(); - md.radaroverlay = radaroverlay(); - md.results = struct(); - md.miscellaneous = miscellaneous(); - md.private = private(); - end - %}}} - function disp(obj) % {{{ - disp(sprintf('%19s: %-22s -- %s','mesh' ,['[1x1 ' class(obj.mesh) ']'],'mesh properties')); - disp(sprintf('%19s: %-22s -- %s','mask' ,['[1x1 ' class(obj.mask) ']'],'defines grounded and floating elements')); - disp(sprintf('%19s: %-22s -- %s','geometry' ,['[1x1 ' class(obj.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...')); - disp(sprintf('%19s: %-22s -- %s','constants' ,['[1x1 ' class(obj.constants) ']'],'physical constants')); - disp(sprintf('%19s: %-22s -- %s','surfaceforcings' ,['[1x1 ' class(obj.surfaceforcings) ']'],'surface forcings')); - disp(sprintf('%19s: %-22s -- %s','basalforcings' ,['[1x1 ' class(obj.basalforcings) ']'],'bed forcings')); - disp(sprintf('%19s: %-22s -- %s','materials' ,['[1x1 ' class(obj.materials) ']'],'material properties')); - disp(sprintf('%19s: %-22s -- %s','friction' ,['[1x1 ' class(obj.friction) ']'],'basal friction/drag properties')); - disp(sprintf('%19s: %-22s -- %s','flowequation' ,['[1x1 ' class(obj.flowequation) ']'],'flow equations')); - disp(sprintf('%19s: %-22s -- %s','timestepping' ,['[1x1 ' class(obj.timestepping) ']'],'time stepping for transient models')); - disp(sprintf('%19s: %-22s -- %s','initialization' ,['[1x1 ' class(obj.initialization) ']'],'initial guess/state')); - disp(sprintf('%19s: %-22s -- %s','rifts' ,['[1x1 ' class(obj.rifts) ']'],'rifts properties')); - disp(sprintf('%19s: %-22s -- %s','debug' ,['[1x1 ' class(obj.debug) ']'],'debugging tools (valgrind, gprof)')); - disp(sprintf('%19s: %-22s -- %s','verbose' ,['[1x1 ' class(obj.verbose) ']'],'verbosity level in solve')); - disp(sprintf('%19s: %-22s -- %s','settings' ,['[1x1 ' class(obj.settings) ']'],'settings properties')); - disp(sprintf('%19s: %-22s -- %s','solver' ,['[1x1 ' class(obj.solver) ']'],'PETSc options for each solution')); - disp(sprintf('%19s: %-22s -- %s','cluster' ,['[1x1 ' class(obj.cluster) ']'],'cluster parameters (number of cpus...)')); - disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(obj.balancethickness) ']'],'parameters for balancethickness solution')); - disp(sprintf('%19s: %-22s -- %s','diagnostic' ,['[1x1 ' class(obj.diagnostic) ']'],'parameters for diagnostic solution')); - disp(sprintf('%19s: %-22s -- %s','groundingline' ,['[1x1 ' class(obj.groundingline) ']'],'parameters for groundingline solution')); - disp(sprintf('%19s: %-22s -- %s','hydrology' ,['[1x1 ' class(obj.hydrology) ']'],'parameters for hydrology solution')); - disp(sprintf('%19s: %-22s -- %s','prognostic' ,['[1x1 ' class(obj.prognostic) ']'],'parameters for prognostic solution')); - disp(sprintf('%19s: %-22s -- %s','thermal' ,['[1x1 ' class(obj.thermal) ']'],'parameters for thermal solution')); - disp(sprintf('%19s: %-22s -- %s','steadystate' ,['[1x1 ' class(obj.steadystate) ']'],'parameters for steadystate solution')); - disp(sprintf('%19s: %-22s -- %s','transient' ,['[1x1 ' class(obj.transient) ']'],'parameters for transient solution')); - disp(sprintf('%19s: %-22s -- %s','autodiff' ,['[1x1 ' class(obj.autodiff) ']'],'automatic differentiation parameters')); - disp(sprintf('%19s: %-22s -- %s','flaim' ,['[1x1 ' class(obj.flaim) ']'],'flaim parameters')); - disp(sprintf('%19s: %-22s -- %s','inversion' ,['[1x1 ' class(obj.inversion) ']'],'parameters for inverse methods')); - disp(sprintf('%19s: %-22s -- %s','qmu' ,['[1x1 ' class(obj.qmu) ']'],'dakota properties')); - disp(sprintf('%19s: %-22s -- %s','results' ,['[1x1 ' class(obj.results) ']'],'model results')); - disp(sprintf('%19s: %-22s -- %s','radaroverlay' ,['[1x1 ' class(obj.radaroverlay) ']'],'radar image for plot overlay')); - disp(sprintf('%19s: %-22s -- %s','miscellaneous' ,['[1x1 ' class(obj.miscellaneous) ']'],'miscellaneous fields')); - end % }}} - end + function md = structtomodel(md,structmd) % {{{ + + if ~isstruct(structmd) error('input model is not a structure'); end + + %loaded model is a struct, initialize output and recover all fields + md = structtoobj(model,structmd); + + %Old field now classes + if (isfield(structmd,'timestepping') & isnumeric(md.timestepping)), md.timestepping=timestepping(); end + if (isfield(structmd,'mask') & isnumeric(md.mask)),md.mask=mask(); end + + %Field name change + if isfield(structmd,'drag'), md.friction.coefficient=structmd.drag; end + if isfield(structmd,'p'), md.friction.p=structmd.p; end + if isfield(structmd,'q'), md.friction.q=structmd.p; end + if isfield(structmd,'melting'), md.basalforcings.melting_rate=structmd.melting; end + if isfield(structmd,'melting_rate'), md.basalforcings.melting_rate=structmd.melting_rate; end + if isfield(structmd,'accumulation'), md.surfaceforcings.mass_balance=structmd.accumulation; end + if isfield(structmd,'numberofgrids'), md.mesh.numberofvertices=structmd.numberofgrids; end + if isfield(structmd,'numberofgrids2d'), md.mesh.numberofvertices2d=structmd.numberofgrids2d; end + if isfield(structmd,'uppergrids'), md.mesh.uppervertex=structmd.uppergrids; end + if isfield(structmd,'lowergrids'), md.mesh.lowervertex=structmd.lowergrids; end + if isfield(structmd,'gridonbed'), md.mesh.vertexonbed=structmd.gridonbed; end + if isfield(structmd,'gridonsurface'), md.mesh.vertexonsurface=structmd.gridonsurface; end + if isfield(structmd,'extractedgrids'), md.mesh.extractedvertices=structmd.extractedgrids; end + if isfield(structmd,'gridoniceshelf'), md.mask.vertexonfloatingice=structmd.gridoniceshelf; end + if isfield(structmd,'gridonicesheet'), md.mask.vertexongroundedice=structmd.gridonicesheet; end + if isfield(structmd,'gridonwater'), md.mask.vertexonwater=structmd.gridonwater; end + if isfield(structmd,'gridonboundary'), md.mesh.vertexonboundary=structmd.gridonboundary; end + if isfield(structmd,'petscoptions') & ~isempty(structmd.petscoptions), md.solver=structmd.petscoptions; end + if isfield(structmd,'g'), md.constants.g=structmd.g; end + if isfield(structmd,'yts'), md.constants.yts=structmd.yts; end + if isfield(structmd,'surface_mass_balance'), md.surfaceforcings.mass_balance=structmd.surface_mass_balance; end + if isfield(structmd,'basal_melting_rate'), md.basalforcings.melting_rate=structmd.basal_melting_rate; end + if isfield(structmd,'basal_melting_rate_correction'), md.basalforcings.melting_rate_correction=structmd.basal_melting_rate_correction; end + if isfield(structmd,'geothermalflux'), md.basalforcings.geothermalflux=structmd.geothermalflux; end + if isfield(structmd,'drag'), md.friction.coefficient=structmd.drag; end + if isfield(structmd,'drag_coefficient'), md.friction.coefficient=structmd.drag_coefficient; end + if isfield(structmd,'drag_p'), md.friction.p=structmd.drag_p; end + if isfield(structmd,'drag_q'), md.friction.q=structmd.drag_q; end + if isfield(structmd,'riftproperties'), %old implementation + md.rifts=rifts(); + md.rifts.riftproperties=structmd.riftproperties; + md.rifts.riftstruct=structmd.rifts; + md.rifts.riftproperties=structmd.riftinfo; + end + if isfield(structmd,'bamg'), md.private.bamg=structmd.bamg; end + if isfield(structmd,'lowmem'), md.settings.lowmem=structmd.lowmem; end + if isfield(structmd,'io_gather'), md.settings.io_gather=structmd.io_gather; end + if isfield(structmd,'spcwatercolumn'), md.hydrology.spcwatercolumn=structmd.spcwatercolumn; end + if isfield(structmd,'hydro_n'), md.hydrology.n=structmd.hydro_n; end + if isfield(structmd,'hydro_p'), md.hydrology.p=structmd.hydro_p; end + if isfield(structmd,'hydro_q'), md.hydrology.q=structmd.hydro_q; end + if isfield(structmd,'hydro_CR'), md.hydrology.CR=structmd.hydro_CR; end + if isfield(structmd,'hydro_kn'), md.hydrology.kn=structmd.hydro_kn; end + if isfield(structmd,'spctemperature'), md.thermal.spctemperature=structmd.spctemperature; end + if isfield(structmd,'min_thermal_constraints'), md.thermal.penalty_threshold=structmd.min_thermal_constraints; end + if isfield(structmd,'artificial_diffusivity'), md.thermal.stabilization=structmd.artificial_diffusivity; end + if isfield(structmd,'max_nonlinear_iterations'), md.thermal.maxiter=structmd.max_nonlinear_iterations; end + if isfield(structmd,'stabilize_constraints'), md.thermal.penalty_lock=structmd.stabilize_constraints; end + if isfield(structmd,'penalty_offset'), md.thermal.penalty_factor=structmd.penalty_offset; end + if isfield(structmd,'name'), md.miscellaneous.name=structmd.name; end + if isfield(structmd,'notes'), md.miscellaneous.notes=structmd.notes; end + if isfield(structmd,'dummy'), md.miscellaneous.dummy=structmd.dummy; end + if isfield(structmd,'dt'), md.timestepping.time_step=structmd.dt; end + if isfield(structmd,'ndt'), md.timestepping.final_time=structmd.ndt; end + if isfield(structmd,'time_adapt'), md.timestepping.time_adapt=structmd.time_adapt; end + if isfield(structmd,'cfl_coefficient'), md.timestepping.cfl_coefficient=structmd.cfl_coefficient; end + if isfield(structmd,'spcthickness'), md.prognostic.spcthickness=structmd.spcthickness; end + if isfield(structmd,'artificial_diffusivity'), md.prognostic.stabilization=structmd.artificial_diffusivity; end + if isfield(structmd,'hydrostatic_adjustment'), md.prognostic.hydrostatic_adjustment=structmd.hydrostatic_adjustment; end + if isfield(structmd,'penalties'), md.prognostic.vertex_pairing=structmd.penalties; end + if isfield(structmd,'penalty_offset'), md.prognostic.penalty_factor=structmd.penalty_offset; end + if isfield(structmd,'B'), md.materials.rheology_B=structmd.B; end + if isfield(structmd,'n'), md.materials.rheology_n=structmd.n; end + if isfield(structmd,'rheology_B'), md.materials.rheology_B=structmd.rheology_B; end + if isfield(structmd,'rheology_n'), md.materials.rheology_n=structmd.rheology_n; end + if isfield(structmd,'elementoniceshelf'), md.mask.elementonfloatingice=structmd.elementoniceshelf; end + if isfield(structmd,'elementonicesheet'), md.mask.elementongroundedice=structmd.elementonicesheet; end + if isfield(structmd,'elementonwater'), md.mask.elementonwater=structmd.elementonwater; end + if isfield(structmd,'nodeoniceshelf'), md.mask.vertexonfloatingice=structmd.nodeoniceshelf; end + if isfield(structmd,'nodeonicesheet'), md.mask.vertexongroundedice=structmd.nodeonicesheet; end + if isfield(structmd,'nodeonwater'), md.mask.vertexonwater=structmd.nodeonwater; end + if isfield(structmd,'spcthickness'), md.balancethickness.spcthickness=structmd.spcthickness; end + if isfield(structmd,'artificial_diffusivity'), md.balancethickness.stabilization=structmd.artificial_diffusivity; end + if isfield(structmd,'dhdt'), md.balancethickness.thickening_rate=structmd.dhdt; end + if isfield(structmd,'ismacayealpattyn'), md.flowequation.ismacayealpattyn=structmd.ismacayealpattyn; end + if isfield(structmd,'ishutter'), md.flowequation.ishutter=structmd.ishutter; end + if isfield(structmd,'isstokes'), md.flowequation.isstokes=structmd.isstokes; end + if isfield(structmd,'elements_type'), md.flowequation.element_equation=structmd.elements_type; end + if isfield(structmd,'vertices_type'), md.flowequation.vertex_equation=structmd.vertices_type; end + if isfield(structmd,'eps_rel'), md.steadystate.reltol=structmd.eps_rel; end + if isfield(structmd,'max_steadystate_iterations'), md.steadystate.maxiter=structmd.max_steadystate_iterations; end + if isfield(structmd,'isdiagnostic'), md.transient.isdiagnostic=structmd.isdiagnostic; end + if isfield(structmd,'isprognostic'), md.transient.isprognostic=structmd.isprognostic; end + if isfield(structmd,'isthermal'), md.transient.isthermal=structmd.isthermal; end + if isfield(structmd,'control_analysis'), md.inversion.iscontrol=structmd.control_analysis; end + if isfield(structmd,'weights'), md.inversion.cost_functions_coefficients=structmd.weights; end + if isfield(structmd,'nsteps'), md.inversion.nsteps=structmd.nsteps; end + if isfield(structmd,'maxiter_per_step'), md.inversion.maxiter_per_step=structmd.maxiter_per_step; end + if isfield(structmd,'cm_min'), md.inversion.min_parameters=structmd.cm_min; end + if isfield(structmd,'cm_max'), md.inversion.max_parameters=structmd.cm_max; end + if isfield(structmd,'vx_obs'), md.inversion.vx_obs=structmd.vx_obs; end + if isfield(structmd,'vy_obs'), md.inversion.vy_obs=structmd.vy_obs; end + if isfield(structmd,'vel_obs'), md.inversion.vel_obs=structmd.vel_obs; end + if isfield(structmd,'thickness_obs'), md.inversion.thickness_obs=structmd.thickness_obs; end + if isfield(structmd,'vx'), md.initialization.vx=structmd.vx; end + if isfield(structmd,'vy'), md.initialization.vy=structmd.vy; end + if isfield(structmd,'vz'), md.initialization.vz=structmd.vz; end + if isfield(structmd,'vel'), md.initialization.vel=structmd.vel; end + if isfield(structmd,'pressure'), md.initialization.pressure=structmd.pressure; end + if isfield(structmd,'temperature'), md.initialization.temperature=structmd.temperature; end + if isfield(structmd,'waterfraction'), md.initialization.waterfraction=structmd.waterfraction; end + if isfield(structmd,'watercolumn'), md.initialization.watercolumn=structmd.watercolumn; end + if isfield(structmd,'surface'), md.geometry.surface=structmd.surface; end + if isfield(structmd,'bed'), md.geometry.bed=structmd.bed; end + if isfield(structmd,'thickness'), md.geometry.thickness=structmd.thickness; end + if isfield(structmd,'bathymetry'), md.geometry.bathymetry=structmd.bathymetry; end + if isfield(structmd,'thickness_coeff'), md.geometry.hydrostatic_ratio=structmd.thickness_coeff; end + if isfield(structmd,'connectivity'), md.mesh.average_vertex_connectivity=structmd.connectivity; end + if isfield(structmd,'extractednodes'), md.mesh.extractedvertices=structmd.extractednodes; end + if isfield(structmd,'extractedelements'), md.mesh.extractedelements=structmd.extractedelements; end + if isfield(structmd,'nodeonboundary'), md.mesh.vertexonboundary=structmd.nodeonboundary; end + if isfield(structmd,'hemisphere'), md.mesh.hemisphere=structmd.hemisphere; end + if isfield(structmd,'lat'), md.mesh.lat=structmd.lat; end + if isfield(structmd,'long'), md.mesh.long=structmd.long; end + if isfield(structmd,'segments'), md.mesh.segments=structmd.segments; end + if isfield(structmd,'segmentmarkers'), md.mesh.segmentmarkers=structmd.segmentmarkers; end + if isfield(structmd,'dim'), md.mesh.dimension=structmd.dim; end + if isfield(structmd,'numlayers'), md.mesh.numberoflayers=structmd.numlayers; end + if isfield(structmd,'numberofelements'), md.mesh.numberofelements=structmd.numberofelements; end + if isfield(structmd,'numberofvertices'), md.mesh.numberofvertices=structmd.numberofvertices; end + if isfield(structmd,'numberofnodes'), md.mesh.numberofvertices=structmd.numberofnodes; end + if isfield(structmd,'numberofedges'), md.mesh.numberofedges=structmd.numberofedges; end + if isfield(structmd,'numberofelements2d'), md.mesh.numberofelements2d=structmd.numberofelements2d; end + if isfield(structmd,'numberofnodes2d'), md.mesh.numberofvertices2d=structmd.numberofnodes2d; end + if isfield(structmd,'nodeconnectivity'), md.mesh.vertexconnectivity=structmd.nodeconnectivity; end + if isfield(structmd,'elementconnectivity'), md.mesh.elementconnectivity=structmd.elementconnectivity; end + if isfield(structmd,'uppernodes'), md.mesh.uppervertex=structmd.uppernodes; end + if isfield(structmd,'lowernodes'), md.mesh.lowervertex=structmd.lowernodes; end + if isfield(structmd,'upperelements'), md.mesh.upperelements=structmd.upperelements; end + if isfield(structmd,'lowerelements'), md.mesh.lowerelements=structmd.lowerelements; end + if isfield(structmd,'elementonbed'), md.mesh.elementonbed=structmd.elementonbed; end + if isfield(structmd,'elementonsurface'), md.mesh.elementonsurface=structmd.elementonsurface; end + if isfield(structmd,'nodeonsurface'), md.mesh.vertexonsurface=structmd.nodeonsurface; end + if isfield(structmd,'nodeonbed'), md.mesh.vertexonbed=structmd.nodeonbed; end + if isfield(structmd,'elements2d'), md.mesh.elements2d=structmd.elements2d; end + if isfield(structmd,'y2d'), md.mesh.y2d=structmd.y2d; end + if isfield(structmd,'x2d'), md.mesh.x2d=structmd.x2d; end + if isfield(structmd,'elements'), md.mesh.elements=structmd.elements; end + if isfield(structmd,'edges'), md.mesh.edges=structmd.edges; end + if isfield(structmd,'y'), md.mesh.y=structmd.y; end + if isfield(structmd,'x'), md.mesh.x=structmd.x; end + if isfield(structmd,'z'), md.mesh.z=structmd.z; end + if isfield(structmd,'mask'), md.flaim.criterion=structmd.mask; end + if isfield(structmd,'pressureload'), md.diagnostic.icefront=structmd.pressureload; end + if isfield(structmd,'diagnostic_ref'), md.diagnostic.referential=structmd.diagnostic_ref; end + if isfield(structmd,'npart'); md.qmu.numberofpartitions=structmd.npart; end + if isfield(structmd,'part'); md.qmu.partition=structmd.part; end + + %Field changes + if (isfield(structmd,'type') & ischar(structmd.type)), + if strcmpi(structmd.type,'2d'), md.mesh.dimension=2; end + if strcmpi(structmd.type,'3d'), md.mesh.dimension=3; end + end + if isnumeric(md.verbose), + md.verbose=verbose; + end + if size(md.diagnostic.icefront,2)==3 || size(md.diagnostic.icefront,2)==5, + front=md.diagnostic.icefront; + md.diagnostic.icefront=[front 1*md.mask.elementonfloatingice(front(:,end))]; + end + if isfield(structmd,'spcvelocity'), + md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1); + md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1); + md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1); + pos=find(structmd.spcvelocity(:,1)); md.diagnostic.spcvx(pos)=structmd.spcvelocity(pos,4); + pos=find(structmd.spcvelocity(:,2)); md.diagnostic.spcvy(pos)=structmd.spcvelocity(pos,5); + pos=find(structmd.spcvelocity(:,3)); md.diagnostic.spcvz(pos)=structmd.spcvelocity(pos,6); + end + if isfield(structmd,'spcvx'), + md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1); + pos=find(~isnan(structmd.spcvx)); md.diagnostic.spcvx(pos)=structmd.spcvx(pos); + end + if isfield(structmd,'spcvy'), + md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1); + pos=find(~isnan(structmd.spcvy)); md.diagnostic.spcvy(pos)=structmd.spcvy(pos); + end + if isfield(structmd,'spcvz'), + md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1); + pos=find(~isnan(structmd.spcvz)); md.diagnostic.spcvz(pos)=structmd.spcvz(pos); + end + if ~isempty(structmd.pressureload) & ismember(structmd.pressureload(end,end),[118 119 120]), + pos=find(structmd.pressureload(:,end)==120); md.diagnostic.icefront(pos,end)=0; + pos=find(structmd.pressureload(:,end)==118); md.diagnostic.icefront(pos,end)=1; + pos=find(structmd.pressureload(:,end)==119); md.diagnostic.icefront(pos,end)=2; + end + if isfield(structmd,'elements_type') & structmd.elements_type(end,end)>50, + pos=find(structmd.elements_type==59); md.flowequation.element_equation(pos,end)=0; + pos=find(structmd.elements_type==55); md.flowequation.element_equation(pos,end)=1; + pos=find(structmd.elements_type==56); md.flowequation.element_equation(pos,end)=2; + pos=find(structmd.elements_type==60); md.flowequation.element_equation(pos,end)=3; + pos=find(structmd.elements_type==62); md.flowequation.element_equation(pos,end)=4; + pos=find(structmd.elements_type==57); md.flowequation.element_equation(pos,end)=5; + pos=find(structmd.elements_type==58); md.flowequation.element_equation(pos,end)=6; + pos=find(structmd.elements_type==61); md.flowequation.element_equation(pos,end)=7; + end + if isfield(structmd,'vertices_type') & structmd.vertices_type(end,end)>50, + pos=find(structmd.vertices_type==59); md.flowequation.vertex_equation(pos,end)=0; + pos=find(structmd.vertices_type==55); md.flowequation.vertex_equation(pos,end)=1; + pos=find(structmd.vertices_type==56); md.flowequation.vertex_equation(pos,end)=2; + pos=find(structmd.vertices_type==60); md.flowequation.vertex_equation(pos,end)=3; + pos=find(structmd.vertices_type==62); md.flowequation.vertex_equation(pos,end)=4; + pos=find(structmd.vertices_type==57); md.flowequation.vertex_equation(pos,end)=5; + pos=find(structmd.vertices_type==58); md.flowequation.vertex_equation(pos,end)=6; + pos=find(structmd.vertices_type==61); md.flowequation.vertex_equation(pos,end)=7; + end + if isfield(structmd,'rheology_law') & isnumeric(structmd.rheology_law), + if (structmd.rheology_law==272), md.materials.rheology_law='None'; end + if (structmd.rheology_law==368), md.materials.rheology_law='Paterson'; end + if (structmd.rheology_law==369), md.materials.rheology_law='Arrhenius'; end + end + if isfield(structmd,'groundingline_migration') & isnumeric(structmd.groundingline_migration), + if (structmd.groundingline_migration==272), md.groundingline.migration='None'; end + if (structmd.groundingline_migration==273), md.groundingline.migration='AgressiveMigration'; end + if (structmd.groundingline_migration==274), md.groundingline.migration='SoftMigration'; end + end + if isfield(structmd,'control_type') & isnumeric(structmd.control_type), + if (structmd.control_type==143), md.inversion.control_parameters={'FrictionCoefficient'}; end + if (structmd.control_type==190), md.inversion.control_parameters={'RheologyBbar'}; end + if (structmd.control_type==147), md.inversion.control_parameters={'Thickeningrate'}; end + end + if isfield(structmd,'cm_responses') & ismember(structmd.cm_responses(end,end),[165:170 383 388 389]), + pos=find(structmd.cm_responses==166); md.inversion.cost_functions(pos)=101; + pos=find(structmd.cm_responses==167); md.inversion.cost_functions(pos)=102; + pos=find(structmd.cm_responses==168); md.inversion.cost_functions(pos)=103; + pos=find(structmd.cm_responses==169); md.inversion.cost_functions(pos)=104; + pos=find(structmd.cm_responses==170); md.inversion.cost_functions(pos)=105; + pos=find(structmd.cm_responses==165); md.inversion.cost_functions(pos)=201; + pos=find(structmd.cm_responses==389); md.inversion.cost_functions(pos)=501; + pos=find(structmd.cm_responses==388); md.inversion.cost_functions(pos)=502; + pos=find(structmd.cm_responses==382); md.inversion.cost_functions(pos)=503; + end + + if isfield(structmd,'artificial_diffusivity') & structmd.artificial_diffusivity==2, + md.thermal.stabilization=2; + md.prognostic.stabilization=1; + md.balancethickness.stabilization=1; + end + if isnumeric(md.prognostic.hydrostatic_adjustment) + if md.prognostic.hydrostatic_adjustment==269, + md.prognostic.hydrostatic_adjustment='Incremental'; + else + md.prognostic.hydrostatic_adjustment='Absolute'; + end + end + + %New fields + if ~isfield(structmd,'upperelements'); + md.mesh.upperelements=transpose(1:md.mesh.numberofelements)+md.mesh.numberofelements2d; + md.mesh.upperelements(end-md.mesh.numberofelements2d+1:end)=NaN; + end + if ~isfield(structmd,'lowerelements'); + md.mesh.lowerelements=transpose(1:md.mesh.numberofelements)-md.mesh.numberofelements2d; + md.mesh.lowerelements(1:md.mesh.numberofelements2d)=NaN; + end + + if ~isfield(structmd,'diagnostic_ref'); + md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6); + end + + end% }}} + function md = setdefaultparameters(md) % {{{ + + %initialize subclasses + md.mesh = mesh(); + md.mask = mask(); + md.constants = constants(); + md.geometry = geometry(); + md.initialization = initialization(); + md.surfaceforcings = surfaceforcings(); + md.basalforcings = basalforcings(); + md.friction = friction(); + md.rifts = rifts(); + md.timestepping = timestepping(); + md.groundingline = groundingline(); + md.materials = matice(); + md.flowequation = flowequation(); + md.debug = debug(); + md.verbose = verbose('solution',true,'qmu',true,'control',true); + md.settings = settings(); + md.solver = solver(); + if ismumps(), + md.solver = addoptions(md.solver,DiagnosticVertAnalysisEnum(),mumpsoptions()); + else + md.solver = addoptions(md.solver,DiagnosticVertAnalysisEnum(),iluasmoptions()); + end + md.cluster = generic(); + md.balancethickness = balancethickness(); + md.diagnostic = diagnostic(); + md.hydrology = hydrology(); + md.prognostic = prognostic(); + md.thermal = thermal(); + md.steadystate = steadystate(); + md.transient = transient(); + md.autodiff = autodiff(); + md.flaim = flaim(); + md.inversion = inversion(); + md.qmu = qmu(); + md.radaroverlay = radaroverlay(); + md.results = struct(); + md.miscellaneous = miscellaneous(); + md.private = private(); + end + %}}} + function disp(obj) % {{{ + disp(sprintf('%19s: %-22s -- %s','mesh' ,['[1x1 ' class(obj.mesh) ']'],'mesh properties')); + disp(sprintf('%19s: %-22s -- %s','mask' ,['[1x1 ' class(obj.mask) ']'],'defines grounded and floating elements')); + disp(sprintf('%19s: %-22s -- %s','geometry' ,['[1x1 ' class(obj.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...')); + disp(sprintf('%19s: %-22s -- %s','constants' ,['[1x1 ' class(obj.constants) ']'],'physical constants')); + disp(sprintf('%19s: %-22s -- %s','surfaceforcings' ,['[1x1 ' class(obj.surfaceforcings) ']'],'surface forcings')); + disp(sprintf('%19s: %-22s -- %s','basalforcings' ,['[1x1 ' class(obj.basalforcings) ']'],'bed forcings')); + disp(sprintf('%19s: %-22s -- %s','materials' ,['[1x1 ' class(obj.materials) ']'],'material properties')); + disp(sprintf('%19s: %-22s -- %s','friction' ,['[1x1 ' class(obj.friction) ']'],'basal friction/drag properties')); + disp(sprintf('%19s: %-22s -- %s','flowequation' ,['[1x1 ' class(obj.flowequation) ']'],'flow equations')); + disp(sprintf('%19s: %-22s -- %s','timestepping' ,['[1x1 ' class(obj.timestepping) ']'],'time stepping for transient models')); + disp(sprintf('%19s: %-22s -- %s','initialization' ,['[1x1 ' class(obj.initialization) ']'],'initial guess/state')); + disp(sprintf('%19s: %-22s -- %s','rifts' ,['[1x1 ' class(obj.rifts) ']'],'rifts properties')); + disp(sprintf('%19s: %-22s -- %s','debug' ,['[1x1 ' class(obj.debug) ']'],'debugging tools (valgrind, gprof)')); + disp(sprintf('%19s: %-22s -- %s','verbose' ,['[1x1 ' class(obj.verbose) ']'],'verbosity level in solve')); + disp(sprintf('%19s: %-22s -- %s','settings' ,['[1x1 ' class(obj.settings) ']'],'settings properties')); + disp(sprintf('%19s: %-22s -- %s','solver' ,['[1x1 ' class(obj.solver) ']'],'PETSc options for each solution')); + disp(sprintf('%19s: %-22s -- %s','cluster' ,['[1x1 ' class(obj.cluster) ']'],'cluster parameters (number of cpus...)')); + disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(obj.balancethickness) ']'],'parameters for balancethickness solution')); + disp(sprintf('%19s: %-22s -- %s','diagnostic' ,['[1x1 ' class(obj.diagnostic) ']'],'parameters for diagnostic solution')); + disp(sprintf('%19s: %-22s -- %s','groundingline' ,['[1x1 ' class(obj.groundingline) ']'],'parameters for groundingline solution')); + disp(sprintf('%19s: %-22s -- %s','hydrology' ,['[1x1 ' class(obj.hydrology) ']'],'parameters for hydrology solution')); + disp(sprintf('%19s: %-22s -- %s','prognostic' ,['[1x1 ' class(obj.prognostic) ']'],'parameters for prognostic solution')); + disp(sprintf('%19s: %-22s -- %s','thermal' ,['[1x1 ' class(obj.thermal) ']'],'parameters for thermal solution')); + disp(sprintf('%19s: %-22s -- %s','steadystate' ,['[1x1 ' class(obj.steadystate) ']'],'parameters for steadystate solution')); + disp(sprintf('%19s: %-22s -- %s','transient' ,['[1x1 ' class(obj.transient) ']'],'parameters for transient solution')); + disp(sprintf('%19s: %-22s -- %s','autodiff' ,['[1x1 ' class(obj.autodiff) ']'],'automatic differentiation parameters')); + disp(sprintf('%19s: %-22s -- %s','flaim' ,['[1x1 ' class(obj.flaim) ']'],'flaim parameters')); + disp(sprintf('%19s: %-22s -- %s','inversion' ,['[1x1 ' class(obj.inversion) ']'],'parameters for inverse methods')); + disp(sprintf('%19s: %-22s -- %s','qmu' ,['[1x1 ' class(obj.qmu) ']'],'dakota properties')); + disp(sprintf('%19s: %-22s -- %s','results' ,['[1x1 ' class(obj.results) ']'],'model results')); + disp(sprintf('%19s: %-22s -- %s','radaroverlay' ,['[1x1 ' class(obj.radaroverlay) ']'],'radar image for plot overlay')); + disp(sprintf('%19s: %-22s -- %s','miscellaneous' ,['[1x1 ' class(obj.miscellaneous) ']'],'miscellaneous fields')); + end % }}} + end end Index: /issm/trunk-jpl/src/m/classes/model/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model/model.py (revision 13691) +++ /issm/trunk-jpl/src/m/classes/model/model.py (revision 13692) @@ -170,4 +170,294 @@ # }}} + """ + function md2 = extract(md,area) % {{{ + %extract - extract a model according to an Argus contour or flag list + % + % This routine extracts a submodel from a bigger model with respect to a given contour + % md must be followed by the corresponding exp file or flags list + % It can either be a domain file (argus type, .exp extension), or an array of element flags. + % If user wants every element outside the domain to be + % extract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp'); + % an empty string '' will be considered as an empty domain + % a string 'all' will be considered as the entire domain + % add an argument 0 if you do not want the elements to be checked (faster) + % + % Usage: + % md2=extract(md,area); + % + % Examples: + % md2=extract(md,'Domain.exp'); + % md2=extract(md,md.mask.elementonfloatingice); + % + % See also: EXTRUDE, COLLAPSE + + %copy model + md1=md; + + %some checks + if ((nargin~=2) | (nargout~=1)), + help extract + error('extract error message: bad usage'); + end + + %get check option + if (nargin==3 & varargin{1}==0), + checkoutline=0; + else + checkoutline=1; + end + + %get elements that are inside area + flag_elem=FlagElements(md1,area); + if ~any(flag_elem), + error('extracted model is empty'); + end + + %kick out all elements with 3 dirichlets + spc_elem=find(~flag_elem); + spc_node=sort(unique(md1.mesh.elements(spc_elem,:))); + flag=ones(md1.mesh.numberofvertices,1); + flag(spc_node)=0; + pos=find(sum(flag(md1.mesh.elements),2)==0); + flag_elem(pos)=0; + + %extracted elements and nodes lists + pos_elem=find(flag_elem); + pos_node=sort(unique(md1.mesh.elements(pos_elem,:))); + + %keep track of some fields + numberofvertices1=md1.mesh.numberofvertices; + numberofelements1=md1.mesh.numberofelements; + numberofvertices2=length(pos_node); + numberofelements2=length(pos_elem); + flag_node=zeros(numberofvertices1,1); + flag_node(pos_node)=1; + + %Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) + Pelem=zeros(numberofelements1,1); + Pelem(pos_elem)=[1:numberofelements2]'; + Pnode=zeros(numberofvertices1,1); + Pnode(pos_node)=[1:numberofvertices2]'; + + %renumber the elements (some node won't exist anymore) + elements_1=md1.mesh.elements; + elements_2=elements_1(pos_elem,:); + elements_2(:,1)=Pnode(elements_2(:,1)); + elements_2(:,2)=Pnode(elements_2(:,2)); + elements_2(:,3)=Pnode(elements_2(:,3)); + if md1.mesh.dimension==3, + elements_2(:,4)=Pnode(elements_2(:,4)); + elements_2(:,5)=Pnode(elements_2(:,5)); + elements_2(:,6)=Pnode(elements_2(:,6)); + end + + %OK, now create the new model ! + + %take every fields from model + md2=md1; + + %automatically modify fields + + %loop over model fields + model_fields=fields(md1); + for i=1:length(model_fields), + %get field + field=md1.(model_fields{i}); + fieldsize=size(field); + if isobject(field), %recursive call + object_fields=fields(md1.(model_fields{i})); + for j=1:length(object_fields), + %get field + field=md1.(model_fields{i}).(object_fields{j}); + fieldsize=size(field); + %size = number of nodes * n + if fieldsize(1)==numberofvertices1 + md2.(model_fields{i}).(object_fields{j})=field(pos_node,:); + elseif (fieldsize(1)==numberofvertices1+1) + md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)]; + %size = number of elements * n + elseif fieldsize(1)==numberofelements1 + md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:); + end + end + else + %size = number of nodes * n + if fieldsize(1)==numberofvertices1 + md2.(model_fields{i})=field(pos_node,:); + elseif (fieldsize(1)==numberofvertices1+1) + md2.(model_fields{i})=[field(pos_node,:); field(end,:)]; + %size = number of elements * n + elseif fieldsize(1)==numberofelements1 + md2.(model_fields{i})=field(pos_elem,:); + end + end + end + + %modify some specific fields + + %Mesh + md2.mesh.numberofelements=numberofelements2; + md2.mesh.numberofvertices=numberofvertices2; + md2.mesh.elements=elements_2; + + %mesh.uppervertex mesh.lowervertex + if md1.mesh.dimension==3 + md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node); + pos=find(~isnan(md2.mesh.uppervertex)); + md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos)); + + md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node); + pos=find(~isnan(md2.mesh.lowervertex)); + md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos)); + + md2.mesh.upperelements=md1.mesh.upperelements(pos_elem); + pos=find(~isnan(md2.mesh.upperelements)); + md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos)); + + md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem); + pos=find(~isnan(md2.mesh.lowerelements)); + md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos)); + end + + %Initial 2d mesh + if md1.mesh.dimension==3 + flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d); + pos_elem_2d=find(flag_elem_2d); + flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d); + pos_node_2d=find(flag_node_2d); + + md2.mesh.numberofelements2d=length(pos_elem_2d); + md2.mesh.numberofvertices2d=length(pos_node_2d); + md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:); + md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1)); + md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2)); + md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3)); + + md2.mesh.x2d=md1.mesh.x(pos_node_2d); + md2.mesh.y2d=md1.mesh.y(pos_node_2d); + end + + %Edges + if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs... + %renumber first two columns + pos=find(md2.mesh.edges(:,4)~=-1); + md2.mesh.edges(: ,1)=Pnode(md2.mesh.edges(:,1)); + md2.mesh.edges(: ,2)=Pnode(md2.mesh.edges(:,2)); + md2.mesh.edges(: ,3)=Pelem(md2.mesh.edges(:,3)); + md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4)); + %remove edges when the 2 vertices are not in the domain. + md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:); + %Replace all zeros by -1 in the last two columns; + pos=find(md2.mesh.edges(:,3)==0); + md2.mesh.edges(pos,3)=-1; + pos=find(md2.mesh.edges(:,4)==0); + md2.mesh.edges(pos,4)=-1; + %Invert -1 on the third column with last column (Also invert first two columns!!) + pos=find(md2.mesh.edges(:,3)==-1); + md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4); + md2.mesh.edges(pos,4)=-1; + values=md2.mesh.edges(pos,2); + md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1); + md2.mesh.edges(pos,1)=values; + %Finally remove edges that do not belong to any element + pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1); + md2.mesh.edges(pos,:)=[]; + end + + %Penalties + if ~isnan(md2.diagnostic.vertex_pairing), + for i=1:size(md1.diagnostic.vertex_pairing,1); + md2.diagnostic.vertex_pairing(i,:)=Pnode(md1.diagnostic.vertex_pairing(i,:)); + end + md2.diagnostic.vertex_pairing=md2.diagnostic.vertex_pairing(find(md2.diagnostic.vertex_pairing(:,1)),:); + end + if ~isnan(md2.prognostic.vertex_pairing), + for i=1:size(md1.prognostic.vertex_pairing,1); + md2.prognostic.vertex_pairing(i,:)=Pnode(md1.prognostic.vertex_pairing(i,:)); + end + md2.prognostic.vertex_pairing=md2.prognostic.vertex_pairing(find(md2.prognostic.vertex_pairing(:,1)),:); + end + + %recreate segments + if md1.mesh.dimension==2 + md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); + md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); + md2.mesh.segments=contourenvelope(md2); + md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; + else + %First do the connectivity for the contourenvelope in 2d + md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d); + md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity); + md2.mesh.segments=contourenvelope(md2); + md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; + md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1); + %Then do it for 3d as usual + md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); + md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); + end + + %Boundary conditions: Dirichlets on new boundary + %Catch the elements that have not been extracted + orphans_elem=find(~flag_elem); + orphans_node=unique(md1.mesh.elements(orphans_elem,:))'; + %Figure out which node are on the boundary between md2 and md1 + nodestoflag1=intersect(orphans_node,pos_node); + nodestoflag2=Pnode(nodestoflag1); + if numel(md1.diagnostic.spcvx)>1 & numel(md1.diagnostic.spcvy)>2 & numel(md1.diagnostic.spcvz)>2, + if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1 + md2.diagnostic.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2); + md2.diagnostic.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2); + else + md2.diagnostic.spcvx(nodestoflag2)=NaN; + md2.diagnostic.spcvy(nodestoflag2)=NaN; + disp(' ') + disp('!! extract warning: spc values should be checked !!') + disp(' ') + end + %put 0 for vz + md2.diagnostic.spcvz(nodestoflag2)=0; + end + if ~isnan(md1.thermal.spctemperature), + md2.thermal.spctemperature(nodestoflag2,1)=1; + end + + %Diagnostic + if ~isnan(md2.diagnostic.icefront) + md2.diagnostic.icefront(:,1)=Pnode(md1.diagnostic.icefront(:,1)); + md2.diagnostic.icefront(:,2)=Pnode(md1.diagnostic.icefront(:,2)); + md2.diagnostic.icefront(:,end-1)=Pelem(md1.diagnostic.icefront(:,end-1)); + if md1.mesh.dimension==3 + md2.diagnostic.icefront(:,3)=Pnode(md1.diagnostic.icefront(:,3)); + md2.diagnostic.icefront(:,4)=Pnode(md1.diagnostic.icefront(:,4)); + end + md2.diagnostic.icefront=md2.diagnostic.icefront(find(md2.diagnostic.icefront(:,1) & md2.diagnostic.icefront(:,2) & md2.diagnostic.icefront(:,end)),:); + end + + %Results fields + if isstruct(md1.results), + md2.results=struct(); + solutionfields=fields(md1.results); + for i=1:length(solutionfields), + %get subfields + solutionsubfields=fields(md1.results.(solutionfields{i})); + for j=1:length(solutionsubfields), + field=md1.results.(solutionfields{i}).(solutionsubfields{j}); + if length(field)==numberofvertices1, + md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_node); + elseif length(field)==numberofelements1, + md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_elem); + else + md2.results.(solutionfields{i}).(solutionsubfields{j})=field; + end + end + end + end + + %Keep track of pos_node and pos_elem + md2.mesh.extractedvertices=pos_node; + md2.mesh.extractedelements=pos_elem; + end % }}} + """ + def extrude(md,*args): # {{{ """