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.dim==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.numlayers); end; if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.numlayers); end; if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.numlayers); end; if ~isnan(md.surfaceforcings.mass_balance), md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.numlayers); end; if ~isnan(md.balancethickness.thickening_rate), md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.numlayers); 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; %bedinfo and surface info md.elementonbed=ones(md.numberofelements2d,1); md.elementonsurface=ones(md.numberofelements2d,1); md.nodeonbed=ones(md.numberofnodes2d,1); md.nodeonsurface=ones(md.numberofnodes2d,1); %elementstype if ~isnan(md.flowequation.element_equation) md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1); end %boundary conditions md.diagnostic.spcvx=project2d(md,md.diagnostic.spcvx,md.numlayers); md.diagnostic.spcvy=project2d(md,md.diagnostic.spcvy,md.numlayers); md.diagnostic.spcvz=project2d(md,md.diagnostic.spcvz,md.numlayers); md.prognostic.spcthickness=project2d(md,md.prognostic.spcthickness,md.numlayers); md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.numlayers); %Extrusion of Neumann BC if ~isnan(md.diagnostic.icefront), numberofneumann2d=size(md.diagnostic.icefront,1)/md.numlayers; 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); %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.connectivity=25; %Collapse the mesh nodes2d=md.numberofnodes2d; elements2d=md.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.nodeonboundary=project2d(md,md.nodeonboundary,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); %Initialize with the 2d mesh md.x=md.x2d; md.y=md.y2d; md.z=zeros(size(md.x2d)); md.numberofnodes=md.numberofnodes2d; md.numberofelements=md.numberofelements2d; md.elements=md.elements2d; %Keep a trace of lower and upper nodes md.lowernodes=NaN; md.uppernodes=NaN; %Remove old mesh md.x2d=NaN; md.y2d=NaN; md.elements2d=NaN; md.numberofelements2d=md.numberofelements; md.numberofnodes2d=md.numberofnodes; md.numlayers=0; %Update mesh type md.dim=2;