0001 function pe=CreatePVector(triaelem,grids,materials,inputs,analysis_type);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 if strcmpi(analysis_type,'diagnostic_horiz'),
0013
0014 pe=CreatePVectorHoriz(triaelem,grids,materials,inputs);
0015
0016 elseif strcmpi(analysis_type,'diagnostic_vert'),
0017
0018 pe=CreatePVectorVert(triaelem,grids,materials,inputs);
0019
0020 elseif strcmpi(analysis_type,'temperature'),
0021
0022 pe=CreatePVectorThermal(triaelem,grids,materials,inputs);
0023
0024 elseif strcmpi(analysis_type,'prognostic'),
0025
0026 pe=CreatePVectorPrognostic(triaelem,grids,materials,inputs);
0027
0028 elseif strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y') | strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
0029
0030 pe=CreatePVectorSlopeCompute(triaelem,grids,materials,inputs,analysis_type);
0031
0032 end
0033 end
0034
0035 function pe=CreatePVectorVert(triaelem,grids,materials,inputs)
0036 error('CreatePVector error message: Vertical analysis not implemented yet');
0037 end
0038
0039 function pe=CreatePVectorThermal(triaelem,grids,materials,inputs)
0040 error('CreatePVector error message: Thermal analysis not implemented yet');
0041 end
0042
0043 function pe=CreatePVectorSlopeCompute(triaelem,grids,materials,inputs,analysis_type);
0044
0045
0046 numgrids=3;
0047 DOFPERGRID=1;
0048 numdof=numgrids*DOFPERGRID;
0049
0050
0051 pe=elemvector(numdof);
0052
0053
0054 xyz_list=getgriddata(triaelem,grids);
0055
0056
0057 xyz_list=xyz_list(1:3,:);
0058
0059
0060 if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y') ,
0061 [param param_is_present]=recover_input(inputs,'bed');
0062 elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
0063 [param param_is_present]=recover_input(inputs,'surface');
0064 end
0065
0066 param_list=zeros(numgrids,1);
0067
0068
0069 for i=1:numgrids,
0070 doflist=grids(triaelem.g(i)).grid.doflist;
0071 if(param_is_present) param_list(i)=param(doflist(1));end;
0072 dof=doflist(1);
0073 pe.row_indices(i)=dof;
0074 end
0075
0076 if param_is_present,
0077 param=param_list(1:3);
0078 else
0079 if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y') ,
0080 param=triaelem.b(1:3);
0081 elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
0082 param=triaelem.s(1:3);
0083 end
0084 end
0085
0086
0087 [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0088
0089
0090 for ig=1:num_gauss2D,
0091
0092
0093 gauss_weight=gauss_weights(ig);
0094 gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0095
0096
0097 dparam=GetParameterDerivativeValue(triaelem,param,xyz_list,gauss_coord);
0098 dparamdx=dparam(1);
0099 dparamdy=dparam(2);
0100
0101
0102 Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
0103
0104
0105 L=GetL(triaelem,gauss_coord,DOFPERGRID);
0106
0107
0108 if strcmp(analysis_type,'bed_slope_compute_x') | strcmp(analysis_type,'surface_slope_compute_x'),
0109 pe_g_gaussian=Jdettria*gauss_weight*dparamdx*L';
0110 elseif strcmp(analysis_type,'bed_slope_compute_y')| strcmp(analysis_type,'surface_slope_compute_y'),
0111 pe_g_gaussian=Jdettria*gauss_weight*dparamdy*L';
0112 end
0113
0114
0115 pe.terms=pe.terms+pe_g_gaussian;
0116 end
0117
0118 end
0119
0120
0121 function pe=CreatePVectorHoriz(triaelem,grids,materials,inputs);
0122
0123 global element_debug element_debugid
0124
0125
0126 if triaelem.acceleration==1,
0127 pe={};
0128 end
0129
0130
0131 numgrids=3;
0132 NDOF2=2;
0133
0134
0135 pe=elemvector(numgrids*NDOF2);
0136
0137
0138 matice=materials(triaelem.matid).material;
0139 matpar=materials(end).constants;
0140
0141
0142 gravity=matpar.g;
0143 rho_ice=matpar.rho_ice;
0144
0145
0146 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0147 [surface_param surface_is_present]=recover_input(inputs,'surface');
0148 [bed_param bed_is_present]=recover_input(inputs,'bed');
0149 [basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
0150
0151
0152 xyz_list=getgriddata(triaelem,grids);
0153
0154
0155 K_list=zeros(numgrids,1);
0156 thickness_list=zeros(numgrids,1);
0157 surface_list=zeros(numgrids,1);
0158 bed_list=zeros(numgrids,1);
0159
0160
0161 for i=1:numgrids,
0162 doflist=grids(triaelem.g(i)).grid.doflist;
0163 for j=1:NDOF2,
0164 pe.row_indices((i-1)*NDOF2+j)=doflist(j);
0165 end
0166 dof=doflist(1);
0167 if(thickness_is_present),
0168 thickness_list(i)=thickness_param(dof);
0169 else
0170 thickness_list(i)=triaelem.h(i);
0171 end
0172 if(surface_is_present),
0173 surface_list(i)=surface_param(dof);
0174 else
0175 surface_list(i)=triaelem.s(i);
0176 end;
0177 if(bed_is_present),
0178 bed_list(i)=bed_param(dof);
0179 else
0180 bed_list(i)=triaelem.b(i);
0181 end
0182 if(basal_drag_is_present),
0183 K_list(i)=basal_drag_param(dof);
0184 else
0185 K_list(i)=triaelem.k(i);
0186 end
0187 end
0188
0189 if (element_debug & triaelem.id==element_debugid),
0190 disp(sprintf('gravity %g\n',gravity));
0191 disp(sprintf('rho_ice %g\n',rho_ice));
0192 disp(sprintf('thickness_list (%g,%g,%g)\n',thickness_list(1),thickness_list(2),thickness_list(3)));
0193 disp(sprintf('surface_list (%g,%g,%g)\n',surface_list(1),surface_list(2),surface_list(3)));
0194 disp(sprintf('bed_list (%g,%g,%g)\n',bed_list(1),bed_list(2),bed_list(3)));
0195 disp(sprintf('K_list (%g,%g,%g)\n',K_list(1),K_list(2),K_list(3)));
0196 end
0197
0198
0199 [num_gauss,first_gauss_area_coord,second_gauss_area_coord,third_gauss_area_coord,gauss_weights]=GaussTria(2);
0200
0201 if (element_debug & triaelem.id==element_debugid),
0202 disp(sprintf(' gaussian points: '));
0203 for i=1:num_gauss,
0204 disp(sprintf(' %g %g %g : %g',first_gauss_area_coord(i),second_gauss_area_coord(i),third_gauss_area_coord(i),gauss_weights(i)));
0205 end
0206 end
0207
0208
0209 for ig=1:num_gauss,
0210
0211 gauss_weight=gauss_weights(ig);
0212 gauss_l1l2l3=[first_gauss_area_coord(ig) second_gauss_area_coord(ig) third_gauss_area_coord(ig)];
0213
0214
0215 thickness=GetParameterValue(triaelem,thickness_list,gauss_l1l2l3);
0216
0217 slopevector=GetParameterDerivativeValue(triaelem,surface_list,xyz_list,gauss_l1l2l3);
0218
0219
0220
0221 if (~triaelem.shelf) & (triaelem.friction_type==1),
0222 plastic_stress=GetParameterValue(triaelem,K_list,gauss_l1l2l3);
0223 end
0224
0225
0226 Jdet=GetJacobianDeterminant(triaelem,xyz_list,gauss_l1l2l3);
0227
0228
0229 l1l2l3=GetNodalFunctions(triaelem,gauss_l1l2l3);
0230
0231
0232 driving_stress_baseline=rho_ice*gravity*thickness;
0233
0234 if (element_debug & triaelem.id==element_debugid),
0235 disp(sprintf(' gaussian %i\n',ig));
0236 disp(sprintf(' thickness %g\n',thickness));
0237 disp(sprintf(' slope(%g,%g)\n',slopevector(1),slopevector(2)));
0238 disp(sprintf(' Jdet %g\n',Jdet));
0239 disp(sprintf(' gaussweigth %g\n',gauss_weight));
0240 disp(sprintf(' l1l2l3 (%g,%g,%g)\n',l1l2l3(1),l1l2l3(2),l1l2l3(3)));
0241 if(triaelem.friction_type==1), disp(sprintf(' plastic_stress(%g)\n',plastic_stress));end;
0242 end
0243
0244
0245 pe_gaussian=zeros(numgrids*NDOF2,1);
0246
0247 if(~triaelem.shelf) & (triaelem.friction_type==1),
0248 for i=1:numgrids,
0249 for j=1:NDOF2,
0250 pe_gaussian((i-1)*NDOF2+j)=(-driving_stress_baseline*slopevector(j)-plastic_stress)*Jdet*gauss_weight*l1l2l3(i);
0251 end
0252 end
0253 else
0254 for i=1:numgrids,
0255 for j=1:NDOF2,
0256 pe_gaussian((i-1)*NDOF2+j)=-driving_stress_baseline*slopevector(j)*Jdet*gauss_weight*l1l2l3(i);
0257 end
0258 end
0259 end
0260
0261
0262 for i=1:pe.nrows,
0263 pe.terms(i)=pe.terms(i)+pe_gaussian(i);
0264 end
0265
0266 end
0267
0268 if (element_debug & triaelem.id==element_debugid),
0269 disp(sprintf(' pe_g->terms\n',ig));
0270 pe.terms
0271 disp(sprintf(' pe_g->row_indices\n',ig));
0272 pe.row_indices
0273 end
0274
0275 end
0276
0277 function pe=CreatePVectorPrognostic(triaelem,grids,materials,inputs);
0278
0279
0280 numgrids=3;
0281 DOFPERGRID=1;
0282 numdof=numgrids*DOFPERGRID;
0283
0284
0285 pe=elemvector(numdof);
0286
0287
0288 xyz_list=getgriddata(triaelem,grids);
0289
0290
0291 [accumulation_param accumulation_is_present]=recover_input(inputs,'accumulation');
0292 [melting_param melting_is_present]=recover_input(inputs,'melting');
0293 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0294 [dt dt_is_present]=recover_input(inputs,'dt');
0295
0296
0297 if ( ~accumulation_is_present | ~melting_is_present | ~thickness_is_present | ~dt_is_present),
0298 error('CreateKMatrixPrognostic error message: missing input parameters!');
0299 end
0300
0301
0302 thickness_list=zeros(numgrids,1);
0303 accumulation_list=zeros(numgrids,1);
0304 melting_list=zeros(numgrids,1);
0305
0306
0307 for i=1:numgrids,
0308 doflist=grids(triaelem.g(i)).grid.doflist;
0309 dof=doflist(1);
0310
0311 pe.row_indices(i)=dof;
0312
0313
0314 thickness_list(i)=thickness_param(dof);
0315 melting_list(i)=melting_param(dof);
0316 accumulation_list(i)=accumulation_param(dof);
0317 end
0318
0319
0320 [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0321
0322 for ig=1:num_gauss2D,
0323
0324
0325 gauss_weight=gauss_weights(ig);
0326 gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0327
0328
0329 melting=GetParameterValue(triaelem,melting_list,gauss_coord);
0330
0331
0332 accumulation=GetParameterValue(triaelem,accumulation_list,gauss_coord);
0333
0334
0335 thickness=GetParameterValue(triaelem,thickness_list,gauss_coord);
0336
0337
0338 Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
0339
0340
0341 L=GetL(triaelem,gauss_coord,DOFPERGRID);
0342
0343
0344 pe_g_gaussian=Jdettria*gauss_weight*(thickness+dt*(accumulation-melting))*L';
0345
0346
0347 pe.terms=pe.terms+pe_g_gaussian;
0348 end
0349
0350 end