source: issm/trunk/src/m/model/collapse.m@ 9681

Last change on this file since 9681 was 9681, checked in by Mathieu Morlighem, 14 years ago

Added inversion

File size: 4.0 KB
RevLine 
[1]1function md=collapse(md)
2%COLLAPSE - collapses a 3d mesh into a 2d mesh
3%
4% This routine collapses a 3d model into a 2d model
5% and collapses all the fileds of the 3d model by
6% taking their depth-averaged values
7%
8% Usage:
9% md=collapse(md)
10%
11% See also: EXTRUDE, MODELEXTRACT
12
13%Check that the model is really a 3d model
[3994]14if ~md.dim==3,
[1]15 error('collapse error message: only 3d mesh can be collapsed')
16end
17
18%Start with changing alle the fields from the 3d mesh
19
[8298]20%drag is limited to nodes that are on the bedrock.
[9610]21md.friction.coefficient=project2d(md,md.friction.coefficient,1);
[1]22
23%p and q (same deal, except for element that are on the bedrock: )
[9610]24md.friction.p=project2d(md,md.friction.p,1);
25md.friction.q=project2d(md,md.friction.q,1);
[1]26
27%observations
[9681]28if ~isnan(md.inversion.vx_obs), md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.numlayers); end;
29if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.numlayers); end;
30if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.numlayers); end;
[9607]31if ~isnan(md.surfaceforcings.mass_balance),
32 md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.numlayers);
33end;
[9646]34if ~isnan(md.balancethickness.thickening_rate), md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.numlayers); end;
[1]35
36%results
37if ~isnan(md.vx),md.vx=DepthAverage(md,md.vx);end;
38if ~isnan(md.vy),md.vy=DepthAverage(md,md.vy);end;
39if ~isnan(md.vz),md.vz=DepthAverage(md,md.vz);end;
40if ~isnan(md.vel),md.vel=DepthAverage(md,md.vel);end;
41
42%bedinfo and surface info
43md.elementonbed=ones(md.numberofelements2d,1);
44md.elementonsurface=ones(md.numberofelements2d,1);
[8298]45md.nodeonbed=ones(md.numberofnodes2d,1);
46md.nodeonsurface=ones(md.numberofnodes2d,1);
[1]47
48%elementstype
[9661]49if ~isnan(md.flowequation.element_equation)
50 md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1);
[1]51end
52
53%boundary conditions
[9679]54md.diagnostic.spcvx=project2d(md,md.diagnostic.spcvx,md.numlayers);
55md.diagnostic.spcvy=project2d(md,md.diagnostic.spcvy,md.numlayers);
56md.diagnostic.spcvz=project2d(md,md.diagnostic.spcvz,md.numlayers);
[9642]57md.prognostic.spcthickness=project2d(md,md.prognostic.spcthickness,md.numlayers);
[9632]58md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.numlayers);
[1]59
60%Extrusion of Neumann BC
[9679]61if ~isnan(md.diagnostic.icefront),
62 numberofneumann2d=size(md.diagnostic.icefront,1)/md.numlayers;
63 md.diagnostic.icefront=[md.diagnostic.icefront(1:numberofneumann2d,1:2) md.diagnostic.icefront(1:numberofneumann2d,5:6)]; %Add two columns on the first layer
[6412]64end
[1]65
66%materials
[9636]67md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B);
68md.materials.rheology_n=project2d(md,md.materials.rheology_n,1);
[1]69
70%special for thermal modeling:
[9612]71md.basalforcings.melting_rate=project2d(md,md.basalforcings.melting_rate,1);
72md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); %bedrock only gets geothermal flux
[1]73
74%update of connectivity matrix
75md.connectivity=25;
76
77%Collapse the mesh
[8298]78nodes2d=md.numberofnodes2d;
[1]79elements2d=md.numberofelements2d;
80
81%parameters
82md.surface=project2d(md,md.surface,1);
83md.thickness=project2d(md,md.thickness,1);
84md.bed=project2d(md,md.bed,1);
[8298]85md.nodeonboundary=project2d(md,md.nodeonboundary,1);
[9641]86md.mask.elementonfloatingice=project2d(md,md.mask.elementonfloatingice,1);
87md.mask.vertexonfloatingice=project2d(md,md.mask.vertexonfloatingice,1);
88md.mask.elementongroundedice=project2d(md,md.mask.elementongroundedice,1);
89md.mask.vertexongroundedice=project2d(md,md.mask.vertexongroundedice,1);
[1]90
91%Initialize with the 2d mesh
92md.x=md.x2d;
93md.y=md.y2d;
[9451]94md.z=zeros(size(md.x2d));
[8298]95md.numberofnodes=md.numberofnodes2d;
[1]96md.numberofelements=md.numberofelements2d;
97md.elements=md.elements2d;
98
[8298]99%Keep a trace of lower and upper nodes
100md.lowernodes=NaN;
101md.uppernodes=NaN;
[1]102
103%Remove old mesh
104md.x2d=NaN;
105md.y2d=NaN;
106md.elements2d=NaN;
107md.numberofelements2d=md.numberofelements;
[8298]108md.numberofnodes2d=md.numberofnodes;
[1]109md.numlayers=0;
110
111%Update mesh type
[3994]112md.dim=2;
Note: See TracBrowser for help on using the repository browser.