[1] | 1 | function 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] | 14 | if ~md.dim==3,
|
---|
[1] | 15 | error('collapse error message: only 3d mesh can be collapsed')
|
---|
| 16 | end
|
---|
| 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] | 21 | md.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] | 24 | md.friction.p=project2d(md,md.friction.p,1);
|
---|
| 25 | md.friction.q=project2d(md,md.friction.q,1);
|
---|
[1] | 26 |
|
---|
| 27 | %observations
|
---|
[6335] | 28 | if ~isnan(md.vx_obs), md.vx_obs=project2d(md,md.vx_obs,md.numlayers); end;
|
---|
| 29 | if ~isnan(md.vy_obs), md.vy_obs=project2d(md,md.vy_obs,md.numlayers); end;
|
---|
| 30 | if ~isnan(md.vel_obs), md.vel_obs=project2d(md,md.vel_obs,md.numlayers); end;
|
---|
[9607] | 31 | if ~isnan(md.surfaceforcings.mass_balance),
|
---|
| 32 | md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.numlayers);
|
---|
| 33 | end;
|
---|
[6335] | 34 | if ~isnan(md.dhdt), md.dhdt=project2d(md,md.dhdt,md.numlayers); end;
|
---|
[1] | 35 |
|
---|
| 36 | %results
|
---|
| 37 | if ~isnan(md.vx),md.vx=DepthAverage(md,md.vx);end;
|
---|
| 38 | if ~isnan(md.vy),md.vy=DepthAverage(md,md.vy);end;
|
---|
| 39 | if ~isnan(md.vz),md.vz=DepthAverage(md,md.vz);end;
|
---|
| 40 | if ~isnan(md.vel),md.vel=DepthAverage(md,md.vel);end;
|
---|
| 41 |
|
---|
| 42 | %bedinfo and surface info
|
---|
| 43 | md.elementonbed=ones(md.numberofelements2d,1);
|
---|
| 44 | md.elementonsurface=ones(md.numberofelements2d,1);
|
---|
[8298] | 45 | md.nodeonbed=ones(md.numberofnodes2d,1);
|
---|
| 46 | md.nodeonsurface=ones(md.numberofnodes2d,1);
|
---|
[1] | 47 |
|
---|
| 48 | %elementstype
|
---|
[9532] | 49 | if ~isnan(md.elements_type)
|
---|
[5898] | 50 | md.elements_type=project2d(md,md.elements_type,1);
|
---|
[1] | 51 | end
|
---|
[8298] | 52 | md.nodeonhutter=project2d(md,md.nodeonhutter,1);
|
---|
| 53 | md.nodeonmacayeal=project2d(md,md.nodeonmacayeal,1);
|
---|
| 54 | md.nodeonpattyn=project2d(md,md.nodeonpattyn,1);
|
---|
| 55 | md.nodeonstokes=project2d(md,md.nodeonstokes,1);
|
---|
[1] | 56 |
|
---|
| 57 | %boundary conditions
|
---|
[8823] | 58 | md.spcvx=project2d(md,md.spcvx,md.numlayers);
|
---|
| 59 | md.spcvy=project2d(md,md.spcvy,md.numlayers);
|
---|
| 60 | md.spcvz=project2d(md,md.spcvz,md.numlayers);
|
---|
[1755] | 61 | md.spcthickness=project2d(md,md.spcthickness,md.numlayers);
|
---|
| 62 | md.spctemperature=project2d(md,md.spctemperature,md.numlayers);
|
---|
[1] | 63 |
|
---|
| 64 | %Extrusion of Neumann BC
|
---|
[6412] | 65 | if ~isnan(md.pressureload),
|
---|
| 66 | numberofneumann2d=size(md.pressureload,1)/md.numlayers;
|
---|
| 67 | md.pressureload=[md.pressureload(1:numberofneumann2d,1:2) md.pressureload(1:numberofneumann2d,5:6)]; %Add two columns on the first layer
|
---|
| 68 | end
|
---|
[1] | 69 |
|
---|
| 70 | %materials
|
---|
[3760] | 71 | md.rheology_B=DepthAverage(md,md.rheology_B);
|
---|
| 72 | md.rheology_n=project2d(md,md.rheology_n,1);
|
---|
[1] | 73 |
|
---|
| 74 | %special for thermal modeling:
|
---|
[9612] | 75 | md.basalforcings.melting_rate=project2d(md,md.basalforcings.melting_rate,1);
|
---|
| 76 | md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); %bedrock only gets geothermal flux
|
---|
[1] | 77 |
|
---|
| 78 | %update of connectivity matrix
|
---|
| 79 | md.connectivity=25;
|
---|
| 80 |
|
---|
| 81 | %Collapse the mesh
|
---|
[8298] | 82 | nodes2d=md.numberofnodes2d;
|
---|
[1] | 83 | elements2d=md.numberofelements2d;
|
---|
| 84 |
|
---|
| 85 | %parameters
|
---|
| 86 | md.surface=project2d(md,md.surface,1);
|
---|
| 87 | md.thickness=project2d(md,md.thickness,1);
|
---|
| 88 | md.bed=project2d(md,md.bed,1);
|
---|
[8298] | 89 | md.nodeonboundary=project2d(md,md.nodeonboundary,1);
|
---|
[1] | 90 | md.elementoniceshelf=project2d(md,md.elementoniceshelf,1);
|
---|
[8298] | 91 | md.nodeoniceshelf=project2d(md,md.nodeoniceshelf,1);
|
---|
[1] | 92 | md.elementonicesheet=project2d(md,md.elementonicesheet,1);
|
---|
[8298] | 93 | md.nodeonicesheet=project2d(md,md.nodeonicesheet,1);
|
---|
[1] | 94 |
|
---|
| 95 | %Initialize with the 2d mesh
|
---|
| 96 | md.x=md.x2d;
|
---|
| 97 | md.y=md.y2d;
|
---|
[9451] | 98 | md.z=zeros(size(md.x2d));
|
---|
[8298] | 99 | md.numberofnodes=md.numberofnodes2d;
|
---|
[1] | 100 | md.numberofelements=md.numberofelements2d;
|
---|
| 101 | md.elements=md.elements2d;
|
---|
| 102 |
|
---|
[8298] | 103 | %Keep a trace of lower and upper nodes
|
---|
| 104 | md.lowernodes=NaN;
|
---|
| 105 | md.uppernodes=NaN;
|
---|
[1] | 106 |
|
---|
| 107 | %Remove old mesh
|
---|
| 108 | md.x2d=NaN;
|
---|
| 109 | md.y2d=NaN;
|
---|
| 110 | md.elements2d=NaN;
|
---|
| 111 | md.numberofelements2d=md.numberofelements;
|
---|
[8298] | 112 | md.numberofnodes2d=md.numberofnodes;
|
---|
[1] | 113 | md.numlayers=0;
|
---|
| 114 |
|
---|
| 115 | %Update mesh type
|
---|
[3994] | 116 | md.dim=2;
|
---|