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

Last change on this file since 6275 was 6275, checked in by seroussi, 14 years ago

minor forgot to collapse dhdt

File size: 3.9 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 grids 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
28md.vx_obs=project2d(md,md.vx_obs,md.numlayers);
29md.vy_obs=project2d(md,md.vy_obs,md.numlayers);
30md.vel_obs=project2d(md,md.vel_obs,md.numlayers);
31md.accumulation_rate=project2d(md,md.accumulation_rate,md.numlayers);
32md.dhdt=project2d(md,md.dhdt,md.numlayers);
33md.firn_layer=project2d(md,md.firn_layer,md.numlayers);
34
35%results
36if ~isnan(md.vx),md.vx=DepthAverage(md,md.vx);end;
37if ~isnan(md.vy),md.vy=DepthAverage(md,md.vy);end;
38if ~isnan(md.vz),md.vz=DepthAverage(md,md.vz);end;
39if ~isnan(md.vel),md.vel=DepthAverage(md,md.vel);end;
40if ~isnan(md.surface_slopex),md.surface_slopex=project2d(md,md.surface_slopex,md.numlayers);end;
41if ~isnan(md.surface_slopey),md.surface_slopey=project2d(md,md.surface_slopey,md.numlayers);end;
42if ~isnan(md.bed_slopex),md.bed_slopex=project2d(md,md.bed_slopex,1);end;
43if ~isnan(md.bed_slopey),md.bed_slopey=project2d(md,md.bed_slopey,1);end;
44
45%bedinfo and surface info
46md.elementonbed=ones(md.numberofelements2d,1);
47md.elementonsurface=ones(md.numberofelements2d,1);
48md.gridonbed=ones(md.numberofgrids2d,1);
49md.gridonsurface=ones(md.numberofgrids2d,1);
50
51%elementstype
52if ~isnan(md.elements_type2d)
53 md.elements_type=md.elements_type2d;
54elseif ~isnan(md.elements_type)
55 md.elements_type=project2d(md,md.elements_type,1);
56end
57md.gridonhutter=project2d(md,md.gridonhutter,1);
58md.gridonmacayeal=project2d(md,md.gridonmacayeal,1);
59md.gridonpattyn=project2d(md,md.gridonpattyn,1);
60md.gridonstokes=project2d(md,md.gridonstokes,1);
61
62%boundary conditions
63md.spcvelocity=project2d(md,md.spcvelocity,md.numlayers);
64md.spcthickness=project2d(md,md.spcthickness,md.numlayers);
65md.spctemperature=project2d(md,md.spctemperature,md.numlayers);
66
67%Extrusion of Neumann BC
68%in 2d, segmentonnumann is: [grid1 grid2 element]
69numberofneumann2d=size(md.pressureload,1)/md.numlayers;
70md.pressureload=[md.pressureload(1:numberofneumann2d,1:2) md.pressureload(1:numberofneumann2d,5:6)]; %Add two columns on the first layer
71
72%materials
73md.rheology_B=DepthAverage(md,md.rheology_B);
74md.rheology_n=project2d(md,md.rheology_n,1);
75
76%special for thermal modeling:
77md.melting_rate=project2d(md,md.melting_rate,1);
78md.observed_temperature=DepthAverage(md,md.observed_temperature);
79md.geothermalflux=project2d(md,md.geothermalflux,1); %bedrock only gets geothermal flux
80
81%update of connectivity matrix
82md.connectivity=25;
83
84%Collapse the mesh
85grids2d=md.numberofgrids2d;
86elements2d=md.numberofelements2d;
87
88%parameters
89md.surface=project2d(md,md.surface,1);
90md.thickness=project2d(md,md.thickness,1);
91md.bed=project2d(md,md.bed,1);
92md.gridonboundary=project2d(md,md.gridonboundary,1);
93md.elementoniceshelf=project2d(md,md.elementoniceshelf,1);
94md.gridoniceshelf=project2d(md,md.gridoniceshelf,1);
95md.elementonicesheet=project2d(md,md.elementonicesheet,1);
96md.gridonicesheet=project2d(md,md.gridonicesheet,1);
97
98%Initialize with the 2d mesh
99md.x=md.x2d;
100md.y=md.y2d;
101md.z=md.z2d;
102md.numberofgrids=md.numberofgrids2d;
103md.numberofelements=md.numberofelements2d;
104md.elements=md.elements2d;
105
106%Keep a trace of lower and upper grids
107md.lowergrids=NaN;
108md.uppergrids=NaN;
109
110%Remove old mesh
111md.x2d=NaN;
112md.y2d=NaN;
113md.z2d=NaN;
114md.elements2d=NaN;
115md.elements_type2d=md.elements_type;
116md.numberofelements2d=md.numberofelements;
117md.numberofgrids2d=md.numberofgrids;
118md.numlayers=0;
119
120%Update mesh type
121md.dim=2;
Note: See TracBrowser for help on using the repository browser.