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

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

removed verticestype2d and elementstype2d
elementstype and verticestype are not enums anymore

File size: 3.6 KB
Line 
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
14if ~md.dim==3,
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
20%drag is limited to nodes that are on the bedrock.
21md.drag_coefficient=project2d(md,md.drag_coefficient,1);
22
23%p and q (same deal, except for element that are on the bedrock: )
24md.drag_p=project2d(md,md.drag_p,1);
25md.drag_q=project2d(md,md.drag_q,1);
26
27%observations
28if ~isnan(md.vx_obs), md.vx_obs=project2d(md,md.vx_obs,md.numlayers); end;
29if ~isnan(md.vy_obs), md.vy_obs=project2d(md,md.vy_obs,md.numlayers); end;
30if ~isnan(md.vel_obs), md.vel_obs=project2d(md,md.vel_obs,md.numlayers); end;
31if ~isnan(md.surface_mass_balance), md.surface_mass_balance=project2d(md,md.surface_mass_balance,md.numlayers); end;
32if ~isnan(md.dhdt), md.dhdt=project2d(md,md.dhdt,md.numlayers); end;
33
34%results
35if ~isnan(md.vx),md.vx=DepthAverage(md,md.vx);end;
36if ~isnan(md.vy),md.vy=DepthAverage(md,md.vy);end;
37if ~isnan(md.vz),md.vz=DepthAverage(md,md.vz);end;
38if ~isnan(md.vel),md.vel=DepthAverage(md,md.vel);end;
39
40%bedinfo and surface info
41md.elementonbed=ones(md.numberofelements2d,1);
42md.elementonsurface=ones(md.numberofelements2d,1);
43md.nodeonbed=ones(md.numberofnodes2d,1);
44md.nodeonsurface=ones(md.numberofnodes2d,1);
45
46%elementstype
47if ~isnan(md.elements_type)
48 md.elements_type=project2d(md,md.elements_type,1);
49end
50md.nodeonhutter=project2d(md,md.nodeonhutter,1);
51md.nodeonmacayeal=project2d(md,md.nodeonmacayeal,1);
52md.nodeonpattyn=project2d(md,md.nodeonpattyn,1);
53md.nodeonstokes=project2d(md,md.nodeonstokes,1);
54
55%boundary conditions
56md.spcvx=project2d(md,md.spcvx,md.numlayers);
57md.spcvy=project2d(md,md.spcvy,md.numlayers);
58md.spcvz=project2d(md,md.spcvz,md.numlayers);
59md.spcthickness=project2d(md,md.spcthickness,md.numlayers);
60md.spctemperature=project2d(md,md.spctemperature,md.numlayers);
61
62%Extrusion of Neumann BC
63if ~isnan(md.pressureload),
64 numberofneumann2d=size(md.pressureload,1)/md.numlayers;
65 md.pressureload=[md.pressureload(1:numberofneumann2d,1:2) md.pressureload(1:numberofneumann2d,5:6)]; %Add two columns on the first layer
66end
67
68%materials
69md.rheology_B=DepthAverage(md,md.rheology_B);
70md.rheology_n=project2d(md,md.rheology_n,1);
71
72%special for thermal modeling:
73md.basal_melting_rate=project2d(md,md.basal_melting_rate,1);
74md.geothermalflux=project2d(md,md.geothermalflux,1); %bedrock only gets geothermal flux
75
76%update of connectivity matrix
77md.connectivity=25;
78
79%Collapse the mesh
80nodes2d=md.numberofnodes2d;
81elements2d=md.numberofelements2d;
82
83%parameters
84md.surface=project2d(md,md.surface,1);
85md.thickness=project2d(md,md.thickness,1);
86md.bed=project2d(md,md.bed,1);
87md.nodeonboundary=project2d(md,md.nodeonboundary,1);
88md.elementoniceshelf=project2d(md,md.elementoniceshelf,1);
89md.nodeoniceshelf=project2d(md,md.nodeoniceshelf,1);
90md.elementonicesheet=project2d(md,md.elementonicesheet,1);
91md.nodeonicesheet=project2d(md,md.nodeonicesheet,1);
92
93%Initialize with the 2d mesh
94md.x=md.x2d;
95md.y=md.y2d;
96md.z=zeros(size(md.x2d));
97md.numberofnodes=md.numberofnodes2d;
98md.numberofelements=md.numberofelements2d;
99md.elements=md.elements2d;
100
101%Keep a trace of lower and upper nodes
102md.lowernodes=NaN;
103md.uppernodes=NaN;
104
105%Remove old mesh
106md.x2d=NaN;
107md.y2d=NaN;
108md.elements2d=NaN;
109md.numberofelements2d=md.numberofelements;
110md.numberofnodes2d=md.numberofnodes;
111md.numlayers=0;
112
113%Update mesh type
114md.dim=2;
Note: See TracBrowser for help on using the repository browser.