0001 function pe=CreatePVector(pentaelem,grids,materials,inputs,analysis_type);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 if strcmpi(analysis_type,'diagnostic_horiz'),
0015
0016 pe=CreatePVectorHoriz(pentaelem,grids,materials,inputs,analysis_type);
0017
0018 elseif strcmpi(analysis_type,'diagnostic_vert'),
0019
0020 pe=CreatePVectorVert(pentaelem,grids,materials,inputs);
0021
0022 elseif strcmpi(analysis_type,'diagnostic_basevert'),
0023
0024 pe=CreatePVectorBaseVert(pentaelem,grids,materials,inputs);
0025
0026 elseif strcmpi(analysis_type,'prognostic'),
0027
0028 pe=CreatePVectorPrognostic(pentaelem,grids,materials,inputs);
0029
0030 elseif (strcmpi(analysis_type,'thermalsteady') | strcmpi(analysis_type,'thermaltransient')),
0031
0032 pe=CreatePVectorThermal(pentaelem,grids,materials,inputs);
0033
0034 elseif strcmpi(analysis_type,'melting')
0035
0036 pe=elemvector(0);
0037
0038 elseif strcmpi(analysis_type,'diagnostic_stokes'),
0039
0040 pe=CreatePVectorStokes(pentaelem,grids,materials,inputs);
0041
0042 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'),
0043
0044 pe=CreatePVectorSlopeCompute(pentaelem,grids,materials,inputs,analysis_type);
0045
0046
0047 end
0048 end
0049
0050 function pe=CreatePVectorStokes(pentaelem,grids,materials,inputs,analysis_type);
0051
0052
0053 numgrids=6;
0054 numgrids2d=3;
0055 NDOF4=4;
0056
0057
0058 pe=elemvector(numgrids*NDOF4);
0059
0060
0061 matice=materials(pentaelem.matid).material;
0062 matpar=materials(end).constants;
0063
0064
0065 gravity=matpar.g;
0066 viscosity_overshoot=matpar.viscosity_overshoot;
0067 rho_ice=matpar.rho_ice;
0068 rho_water=matpar.rho_water;
0069
0070
0071 conditioning=pentaelem.reconditioning_number;
0072
0073
0074 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0075 [temperature_param temperature_is_present]=recover_input(inputs,'temperature');
0076 [surface_param surface_is_present]=recover_input(inputs,'surface');
0077 [bed_param bed_is_present]=recover_input(inputs,'bed');
0078 [melting_param melting_is_present]=recover_input(inputs,'melting');
0079 [basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
0080 [flow_law_param flow_law_is_present]=recover_input(inputs,'B');
0081 [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0082
0083
0084 xyz_list=getgriddata(pentaelem,grids);
0085
0086
0087 if velocity_is_present,
0088 vxvyvz_list=zeros(numgrids,3);
0089 end
0090 B_list=zeros(numgrids,1);
0091 basal_drag_list=zeros(numgrids,1);
0092 melting_list=zeros(numgrids,1);
0093 thickness_list=zeros(numgrids,1);
0094 surface_list=zeros(numgrids,1);
0095 bed_list=zeros(numgrids,1);
0096 temperature_list=zeros(numgrids,1);
0097
0098
0099 for i=1:numgrids,
0100
0101 doflist=grids(pentaelem.g(i)).grid.doflist;
0102 for j=1:NDOF4,
0103 dof=doflist(j);
0104 pe.row_indices((i-1)*NDOF4+j)=dof;
0105 if j<4,
0106 if velocity_is_present, vxvyvz_list(i,j)=velocity_param(dof); end;
0107 end
0108 end
0109 dof=doflist(1);
0110 if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
0111 if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
0112 if(temperature_is_present) temperature_list(i)=temperature_param(dof);end;
0113 if(surface_is_present) surface_list(i)=surface_param(dof);end;
0114 if(bed_is_present) bed_list(i)=bed_param(dof);end;
0115 if(melting_is_present) melting_list(i)=melting_param(dof);end;
0116 if(basal_drag_is_present) basal_drag_list(i)=basal_drag_param(dof);end;
0117
0118 end
0119
0120
0121 area_order=5;
0122 num_vert_gauss=5;
0123 [num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);
0124
0125
0126
0127
0128
0129 pe_temp=zeros(27,1);
0130 Kebubble_temp=zeros(27,3);
0131
0132
0133 for igarea=1:num_area_gauss,
0134 for igvert=1:num_vert_gauss,
0135
0136 gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
0137 gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
0138
0139
0140
0141 if (~pentaelem.shelf) & (pentaelem.friction_type==1),
0142 if basal_drag_is_present,
0143 plastic_stress=GetParameterValueStokes(pentaelem,basal_drag_list,gauss_coord);
0144 else
0145 plastic_stress=GetParameterValueStokes(pentaelem,pentaelem.k,gauss_coord);
0146 end
0147 end
0148
0149
0150 Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
0151
0152
0153
0154 l1l7=GetNodalFunctionsStokes(pentaelem,gauss_coord);
0155
0156
0157 pe_gaussian=zeros(numgrids*NDOF4+3,1);
0158
0159 if(~pentaelem.shelf) & (pentaelem.friction_type==1),
0160 for i=1:numgrids+1,
0161 for j=1:2,
0162 pe_gaussian((i-1)*NDOF4+j)=-plastic_stress*Jdet*gauss_weight*l1l7(i);
0163 end
0164 pe_gaussian((i-1)*NDOF4+3)=-rho_ice*gravity*Jdet*gauss_weight*l1l7(i);
0165 end
0166 else
0167 for i=1:numgrids+1,
0168 pe_gaussian((i-1)*NDOF4+3)=-rho_ice*gravity*Jdet*gauss_weight*l1l7(i);
0169 end
0170 end
0171
0172
0173 pe_temp=pe_temp+pe_gaussian;
0174
0175
0176
0177
0178 if velocity_is_present,
0179 epsilon=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord);
0180 else
0181 epsilon=zeros(6,1);
0182 end
0183
0184
0185
0186
0187
0188
0189
0190 if temperature_is_present,
0191 temperature=GetParameterValueStokes(pentaelem,temperature_list,gauss_coord);
0192 matice.B=paterson(temperature);
0193 end
0194
0195
0196
0197 if flow_law_is_present,
0198 B_param=GetParameterValueStokes(pentaelem,B_list,gauss_coord);
0199 matice.B=B_param; clear B_param.
0200 end
0201
0202
0203 viscosity=GetViscosity3d(matice,epsilon);
0204
0205
0206
0207 B=GetBStokes(pentaelem,xyz_list,gauss_coord);
0208 Bprime=GetBprimeStokes(pentaelem,xyz_list,gauss_coord);
0209 Bprime_bubble=Bprime(:,25:27);
0210
0211
0212
0213 D=gauss_weight*Jdet*diag([viscosity*ones(6,1);-conditioning*ones(2,1)]);
0214
0215
0216 Ke_gg_gaussian=B'*D*Bprime_bubble;
0217
0218
0219 Kebubble_temp=Kebubble_temp+Ke_gg_gaussian;
0220
0221 end
0222 end
0223
0224
0225
0226 if (pentaelem.onbed==1 & pentaelem.shelf==1),
0227
0228
0229 if thickness_is_present,
0230 thickness=thickness_list(1:3);
0231 else
0232 thickness=pentaelem.h(1:3);
0233 end
0234 if bed_is_present,
0235 bed=bed_list(1:3);
0236 else
0237 bed=pentaelem.b(1:3);
0238 end
0239 if basal_drag_is_present,
0240 basal_drag=basal_drag_list(1:3);
0241 else
0242 basal_drag=pentaelem.k(1:3);
0243 end
0244
0245
0246 xyz_list_tria=xyz_list(1:3,:);
0247
0248
0249 [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0250
0251
0252 n=zeros(3,1);
0253 normal=cross((xyz_list_tria(1,:)-xyz_list_tria(3,:)),(xyz_list_tria(3,:)-xyz_list_tria(2,:)));
0254 n(1)=1/norm(normal)*normal(1);
0255 n(2)=1/norm(normal)*normal(2);
0256 n(3)=1/norm(normal)*normal(3);
0257
0258 pe_gg_water=zeros(27,1);
0259
0260 for ig=1:num_gauss2D,
0261
0262
0263 gauss_weight=gauss_weights(ig);
0264 gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0265
0266
0267 Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord);
0268
0269
0270 bed_gg=GetParameterValue(triaelem,bed,gauss_coord);
0271
0272
0273 L=GetL(triaelem,gauss_coord,1);
0274
0275
0276 pressure=gravity*rho_water*bed_gg;
0277
0278
0279 for i=1:numgrids2d,
0280 for j=1:3,
0281 pe_gg_water((i-1)*NDOF4+j)=pe_gg_water((i-1)*NDOF4+j)+pressure*gauss_weight*Jdettria*L(i)*n(j);
0282 end
0283 end
0284 end
0285
0286
0287 pe_temp=pe_temp+pe_gg_water;
0288
0289 end
0290
0291
0292 pe_reduced=ReduceVectorStokes(pentaelem,Kebubble_temp, pe_temp);
0293
0294 pe.terms=pe.terms+pe_reduced;
0295
0296 end
0297
0298
0299 function pe=CreatePVectorSlopeCompute(pentaelem,grids,materials,inputs,analysis_type);
0300
0301
0302
0303 if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y'),
0304 if pentaelem.onbed~=1,
0305 pe=elemvector(0);
0306 return;
0307 end
0308 elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
0309 if pentaelem.onsurface~=1,
0310 pe=elemvector(0);
0311 return;
0312 end
0313 end
0314
0315
0316 numgrids=3;
0317 DOFPERGRID=1;
0318 numdof=numgrids*DOFPERGRID;
0319
0320
0321 pe=elemvector(numdof);
0322
0323
0324 xyz_list=getgriddata(pentaelem,grids);
0325
0326
0327 if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y'),
0328 xyz_list=xyz_list(1:3,:);
0329 [param param_is_present]=recover_input(inputs,'bed');
0330 elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
0331 xyz_list=xyz_list(4:6,:);
0332 [param param_is_present]=recover_input(inputs,'surface');
0333 end
0334
0335 param_list=zeros(numgrids,1);
0336
0337
0338 for i=1:numgrids,
0339 if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y'),
0340 doflist=grids(pentaelem.g(i)).grid.doflist;
0341 elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
0342 doflist=grids(pentaelem.g(i+3)).grid.doflist;
0343 end
0344 if(param_is_present) param_list(i)=param(doflist(1));end;
0345 dof=doflist(1);
0346 pe.row_indices(i)=dof;
0347 end
0348
0349 if param_is_present,
0350 param=param_list(1:3);
0351 else
0352 if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y'),
0353 param=pentaelem.b(1:3);
0354 elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
0355 param=pentaelem.s(4:6);
0356 end
0357 end
0358
0359
0360 [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0361
0362
0363 for ig=1:num_gauss2D,
0364
0365
0366 gauss_weight=gauss_weights(ig);
0367 gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0368
0369
0370 dparam=GetParameterDerivativeValue(triaelem,param,xyz_list,gauss_coord);
0371 dparamdx=dparam(1);
0372 dparamdy=dparam(2);
0373
0374
0375 Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
0376
0377
0378 L=GetL(triaelem,gauss_coord,DOFPERGRID);
0379
0380
0381 if strcmp(analysis_type,'bed_slope_compute_x') | strcmp(analysis_type,'surface_slope_compute_x'),
0382 pe_g_gaussian=Jdettria*gauss_weight*dparamdx*L';
0383 elseif strcmp(analysis_type,'bed_slope_compute_y') | strcmp(analysis_type,'surface_slope_compute_y'),
0384 pe_g_gaussian=Jdettria*gauss_weight*dparamdy*L';
0385 end
0386
0387
0388 pe.terms=pe.terms+pe_g_gaussian;
0389 end
0390
0391 end
0392
0393
0394
0395 function pe=CreatePVectorBaseVert(pentaelem,grids,materials,inputs);
0396
0397
0398
0399 if pentaelem.onbed~=1,
0400 pe=elemvector(0);
0401 return;
0402 end
0403
0404
0405 numgrids=3;
0406 DOFPERGRID=1;
0407 numdof=numgrids*DOFPERGRID;
0408
0409
0410 matice=materials(pentaelem.matid).material;
0411 matpar=materials(end).constants;
0412 rho_ice=matpar.rho_ice;
0413 rho_water=matpar.rho_water;
0414 di=rho_ice/rho_water;
0415
0416
0417 pe=elemvector(numdof);
0418
0419
0420 xyz_list=getgriddata(pentaelem,grids);
0421
0422
0423 xyz_list=xyz_list(1:3,:);
0424
0425
0426 [bed_param bed_is_present]=recover_input(inputs,'bed');
0427 [melting_param melting_is_present]=recover_input(inputs,'melting');
0428 [accumulation_param accumulation_is_present]=recover_input(inputs,'accumulation');
0429 [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0430 [velocity_average_param velocity_average_is_present]=recover_input(inputs,'velocity_average');
0431 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0432 if (~velocity_is_present | ~velocity_average_is_present),
0433 error('CreatePVector error message: no input velocity!');
0434 end
0435
0436
0437 vxvy_list=zeros(numgrids,2);
0438 vxvy_average_list=zeros(numgrids,2);
0439 bed_list=zeros(numgrids,1);
0440 melting_list=zeros(numgrids,1);
0441 accumulation_list=zeros(numgrids,1);
0442 thickness_list=zeros(numgrids,1);
0443
0444
0445 for i=1:numgrids,
0446 doflist=grids(pentaelem.g(i)).grid.doflist;
0447 for j=1:2,
0448 dof=doflist(j);
0449 vxvy_list(i,j)=velocity_param(dof);
0450 vxvy_average_list(i,j)=velocity_average_param(dof);
0451 end
0452 if(bed_is_present) bed_list(i)=bed_param(doflist(1));end;
0453 if(melting_is_present) melting_list(i)=melting_param(doflist(1));end;
0454 if(accumulation_is_present) accumulation_list(i)=accumulation_param(doflist(1));end;
0455 if(thickness_is_present) thickness_list(i)=thickness_param(doflist(1));end;
0456 dof=doflist(3);
0457 pe.row_indices(i)=dof;
0458 end
0459
0460 if bed_is_present, bed=bed_list(1:3); else bed=pentaelem.b(1:3); end
0461 if thickness_is_present, thickness=thickness_list(1:3); else thickness=pentaelem.h(1:3); end
0462 if melting_is_present, melting=melting_list(1:3); else melting=pentaelem.melting(1:3); end
0463 if accumulation_is_present, accumulation=accumulation_list(1:3); else accumulation=pentaelem.accumulation(1:3); end
0464
0465
0466 [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0467
0468
0469 for ig=1:num_gauss2D,
0470
0471
0472 gauss_weight=gauss_weights(ig);
0473 gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0474
0475
0476 melting_gauss=GetParameterValue(triaelem,melting,gauss_coord);
0477
0478
0479 vx=GetParameterValue(triaelem,vxvy_list(:,1),gauss_coord);
0480 vy=GetParameterValue(triaelem,vxvy_list(:,2),gauss_coord);
0481
0482
0483 db=GetParameterDerivativeValue(triaelem,bed,xyz_list,gauss_coord);
0484 dbdx=db(1);
0485 dbdy=db(2);
0486
0487
0488 Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
0489
0490
0491 L=GetL(triaelem,gauss_coord,DOFPERGRID);
0492
0493
0494 pe_g_basevert_gaussian=Jdettria*gauss_weight*(vx*dbdx+vy*dbdy-melting_gauss)*L';
0495
0496
0497 pe.terms=pe.terms+pe_g_basevert_gaussian;
0498 end
0499
0500
0501
0502 if pentaelem.shelf,
0503
0504
0505 HU=thickness.*vxvy_average_list(:,1);
0506 HV=thickness.*vxvy_average_list(:,2);
0507
0508 for ig=1:num_gauss2D,
0509
0510
0511 gauss_weight=gauss_weights(ig);
0512 gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0513
0514
0515 dHU=GetParameterDerivativeValue(triaelem,HU,xyz_list,gauss_coord);
0516 dHV=GetParameterDerivativeValue(triaelem,HV,xyz_list,gauss_coord);
0517 divHVel=dHU(1)+dHV(2);
0518
0519
0520 melting_gauss=GetParameterValue(triaelem,melting,gauss_coord);
0521 accumulation_gauss=GetParameterValue(triaelem,accumulation,gauss_coord);
0522
0523
0524 L=GetL(triaelem,gauss_coord,DOFPERGRID);
0525
0526
0527 Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
0528
0529
0530 pe_g_basevert_iceshelf_gaussian=Jdettria*gauss_weight*(-di*(-divHVel+accumulation_gauss-melting_gauss))*L';
0531
0532
0533 pe.terms=pe.terms+pe_g_basevert_iceshelf_gaussian;
0534 end
0535 end
0536 end
0537
0538 function pe=CreatePVectorVert(pentaelem,grids,materials,inputs);
0539
0540
0541 numgrids=6;
0542 NDOF1=1;
0543 NDOF2=2;
0544
0545
0546 pe=elemvector(numgrids*NDOF1);
0547
0548
0549 [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0550 if ~velocity_is_present,
0551 error('CreatePVector error message: no input velocity!');
0552 end
0553
0554
0555 vxvy_list=zeros(numgrids,2);
0556
0557
0558 xyz_list=getgriddata(pentaelem,grids);
0559
0560
0561 for i=1:numgrids,
0562 doflist=grids(pentaelem.g(i)).grid.doflist;
0563 for j=1:NDOF2,
0564 dof=doflist(j);
0565 vxvy_list(i,j)=velocity_param(dof);
0566 end
0567 dof=doflist(3);
0568 pe.row_indices(i)=dof;
0569 end
0570
0571
0572 area_order=2;
0573 num_vert_gauss=2;
0574 [num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);
0575
0576
0577 for igarea=1:num_area_gauss,
0578 for igvert=1:num_vert_gauss,
0579
0580 gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
0581 gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
0582
0583
0584 dudx=GetParameterDerivativeValue(pentaelem,vxvy_list(:,1),xyz_list,gauss_coord); dudx=dudx(1);
0585 dvdy=GetParameterDerivativeValue(pentaelem,vxvy_list(:,2),xyz_list,gauss_coord); dvdy=dvdy(2);
0586
0587
0588 Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
0589
0590
0591
0592 l1l6=GetNodalFunctions(pentaelem,gauss_coord);
0593
0594 pe.terms=pe.terms+l1l6'*(dudx+dvdy)*Jdet*gauss_weight;
0595
0596 end
0597 end
0598
0599 end
0600
0601 function pe=CreatePVectorThermal(pentaelem,grids,materials,inputs)
0602
0603
0604 NDOF1=1;
0605 numgrids=6;
0606
0607
0608 pe=elemvector(numgrids*NDOF1);
0609
0610
0611 matice=materials(pentaelem.matid).material;
0612 matpar=materials(end).constants;
0613
0614 rho_ice=matpar.rho_ice;
0615 rho_water=matpar.rho_water;
0616 gravity=matpar.g;
0617 heatcapacity=matpar.heatcapacity;
0618 beta=matpar.beta;
0619 meltingpoint=matpar.meltingpoint;
0620
0621
0622 xyz_list=getgriddata(pentaelem,grids);
0623
0624
0625 [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0626 [flow_law_param flow_law_is_present]=recover_input(inputs,'B');
0627 [basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
0628 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0629 [surface_param surface_is_present]=recover_input(inputs,'surface');
0630 [bed_param bed_is_present]=recover_input(inputs,'bed');
0631 [temperature_param temperature_is_present]=recover_input(inputs,'temperature');
0632 [dt dt_is_present]=recover_input(inputs,'dt');
0633 [pressure_param pressure_is_present]=recover_input(inputs,'pressure');
0634
0635
0636
0637 if ~velocity_is_present,
0638 error('CreatePVectorThermal error message: input velocity not present!');
0639 end
0640 if ~dt_is_present & ~pentaelem.thermal_steadystate,
0641 error('CreatePVectorThermal error message: transient thermal solution requires present of time step as parameter!');
0642 end
0643 if ~temperature_is_present & ~pentaelem.thermal_steadystate,
0644 error('CreatePVectorThermal error message: transient thermal solution requires temperature at previous iteration!');
0645 end
0646 if ~pressure_is_present,
0647 error('CreatePVectorThermal error message: thermal solution requires pressure');
0648 end
0649
0650
0651
0652
0653 vxvyvz_list=zeros(numgrids,3);
0654 B_list=zeros(numgrids,1);
0655 K_list=zeros(numgrids,1);
0656 thickness_list=zeros(numgrids,1);
0657 surface_list=zeros(numgrids,1);
0658 bed_list=zeros(numgrids,1);
0659 temperature_list=zeros(numgrids,1);
0660 pressure_list=zeros(numgrids,1);
0661
0662
0663
0664 for i=1:numgrids,
0665 doflist=grids(pentaelem.g(i)).grid.doflist;
0666 for j=1:3,
0667 dof=doflist(j);
0668 vxvyvz_list(i,j)=velocity_param(dof);
0669 end
0670 dof=doflist(1);
0671 pe.row_indices(i)=dof;
0672 if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
0673 if(basal_drag_is_present), K_list(i)=basal_drag_param(dof);end;
0674 if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
0675 if(surface_is_present) surface_list(i)=surface_param(dof);end;
0676 if(bed_is_present) bed_list(i)=bed_param(dof);end;
0677 if(temperature_is_present) temperature_list(i)=temperature_param(dof);end;
0678 if (pressure_is_present) pressure_list(i)=pressure_param(dof); end;
0679 end
0680 vxvy_list= vxvyvz_list(:,1:2);
0681
0682
0683
0684
0685
0686
0687 area_order=2;
0688 num_vert_gauss=3;
0689 [num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);
0690
0691
0692 for igarea=1:num_area_gauss,
0693 for igvert=1:num_vert_gauss,
0694
0695 gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
0696 gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
0697
0698
0699 if temperature_is_present,
0700 temperature=GetParameterValue(pentaelem,temperature_list,gauss_coord);
0701 matice.B=paterson(temperature);
0702 end
0703
0704
0705 epsilon=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord);
0706 viscosity=GetViscosity3d(matice,epsilon);
0707
0708
0709 epsilon_matrix=[epsilon(1) epsilon(4) epsilon(5)
0710 epsilon(4) epsilon(2) epsilon(6)
0711 epsilon(5) epsilon(6) epsilon(3)];
0712
0713 epsilon_e=effective_value(epsilon_matrix);
0714
0715 phi=2*viscosity*epsilon_e^2;
0716
0717
0718 Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
0719
0720
0721 l1l6=GetNodalFunctions(pentaelem,gauss_coord);
0722
0723
0724 pe_gaussian_deformational=zeros(numgrids*NDOF1,1);
0725
0726 if pentaelem.thermal_steadystate,
0727 pe_gaussian_deformational=phi/(rho_ice*heatcapacity)*Jdet*gauss_weight*l1l6';
0728 else
0729 pe_gaussian_deformational=dt*phi/(rho_ice*heatcapacity)*Jdet*gauss_weight*l1l6';
0730 end
0731
0732
0733 if pentaelem.thermal_steadystate,
0734 pe_gaussian_transient=0;
0735 else
0736 pe_gaussian_transient=temperature*Jdet*gauss_weight*l1l6';
0737 end
0738
0739
0740 pe.terms=pe.terms+pe_gaussian_deformational+pe_gaussian_transient;
0741
0742 end
0743 end
0744
0745
0746 if (pentaelem.onbed & pentaelem.shelf),
0747
0748
0749 mixed_layer_capacity=matpar.mixed_layer_capacity;
0750 thermal_exchange_velocity=matpar.thermal_exchange_velocity;
0751
0752
0753 pressure_tria=pressure_list(1:3);
0754
0755
0756 [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0757
0758 for ig=1:num_gauss2D,
0759
0760
0761 gauss_weight=gauss_weights(ig);
0762 gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0763
0764
0765 xyz_list_tria=xyz_list(1:3,:);
0766 Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord);
0767
0768
0769 l1l2l3=GetNodalFunctions(triaelem,gauss_coord);
0770
0771
0772 pressure_g=GetParameterValue(triaelem,pressure_tria,gauss_coord);
0773 t_pmp=meltingpoint-beta*pressure_g;
0774
0775
0776 if pentaelem.thermal_steadystate,
0777 pe_ocean_gaussian=gauss_weight*Jdettria*rho_water*mixed_layer_capacity*thermal_exchange_velocity*t_pmp/(heatcapacity*rho_ice)*l1l2l3';
0778 else
0779 pe_ocean_gaussian=dt*gauss_weight*Jdettria*rho_water*mixed_layer_capacity*thermal_exchange_velocity*t_pmp/(heatcapacity*rho_ice)*l1l2l3';
0780 end
0781
0782
0783 pe.terms(1:3)=pe.terms(1:3)+pe_ocean_gaussian;
0784 end
0785 end
0786
0787
0788 if (pentaelem.onbed & ~pentaelem.shelf),
0789
0790
0791 geothermalflux=pentaelem.geothermalflux(1:3,1);
0792
0793
0794 basal_friction=zeros(3,1);
0795 if (pentaelem.friction_type==2),
0796
0797
0798 frictionparameters=struct();
0799 frictionparameters.element_type='3d';
0800 frictionparameters.rho_ice=rho_ice;
0801 frictionparameters.rho_water=rho_water;
0802 frictionparameters.g=gravity;
0803 frictionparameters.p=pentaelem.p;
0804 frictionparameters.q=pentaelem.q;
0805 frictionparameters.velocities=vxvy_list(1:3,:);
0806
0807 if thickness_is_present,
0808 frictionparameters.h=thickness_list(1:3,:);
0809 else
0810 frictionparameters.h=pentaelem.h(1:3,:);
0811 end
0812 if(bed_is_present),
0813 frictionparameters.b=bed_list(1:3,:);
0814 else
0815 frictionparameters.b=pentaelem.b(1:3,:);
0816 end
0817 if basal_drag_is_present,
0818 frictionparameters.k=K_list(1:3,:);
0819 else
0820 frictionparameters.k=pentaelem.k(1:3,:);
0821 end
0822
0823 alpha2=Getalpha2(frictionparameters);
0824
0825
0826 u_b2=vxvy_list(1:3,1).^2+vxvy_list(1:3,2).^2;
0827
0828
0829 basal_friction=alpha2.*u_b2;
0830 end
0831
0832
0833 [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0834
0835 for ig=1:num_gauss2D,
0836
0837
0838 gauss_weight=gauss_weights(ig);
0839 gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0840
0841
0842 xyz_list_tria=xyz_list(1:3,:);
0843 Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord);
0844
0845
0846 G_g=GetParameterValue(triaelem,geothermalflux,gauss_coord);
0847 basal_friction_g=GetParameterValue(triaelem,basal_friction,gauss_coord);
0848
0849
0850 l1l2l3=GetNodalFunctions(triaelem,gauss_coord);
0851
0852
0853 if pentaelem.thermal_steadystate,
0854 pe_geo_gaussian=gauss_weight*Jdettria*(basal_friction_g+G_g)/(heatcapacity*rho_ice)*l1l2l3';
0855 else
0856 pe_geo_gaussian=dt*gauss_weight*Jdettria*(basal_friction_g+G_g)/(heatcapacity*rho_ice)*l1l2l3';
0857 end
0858
0859
0860 pe.terms(1:3)=pe.terms(1:3)+pe_geo_gaussian;
0861 end
0862 end
0863
0864
0865 end
0866
0867
0868 function pe=CreatePVectorHoriz(pentaelem,grids,materials,inputs,analysis_type);
0869
0870
0871
0872
0873
0874 if pentaelem.collapse & ~pentaelem.onbed
0875
0876 pe=elemvector(0);
0877 return;
0878
0879 elseif pentaelem.collapse & pentaelem.onbed
0880
0881
0882 element=triaelem;
0883 element.type='triaelem';
0884 element.id=NaN;
0885 element.matid=pentaelem.matid;
0886 element.g=pentaelem.g(1:3);
0887 element.h=pentaelem.h(1:3);
0888 element.s=pentaelem.s(1:3);
0889 element.b=pentaelem.b(1:3);
0890 element.friction_type=pentaelem.friction_type;
0891 element.k=pentaelem.k(1:3);
0892 element.p=pentaelem.p;
0893 element.q=pentaelem.q;
0894 element.shelf=pentaelem.shelf;
0895 element.meanvel=pentaelem.meanvel;
0896 element.epsvel=pentaelem.epsvel;
0897 element.acceleration=pentaelem.acceleration;
0898
0899
0900 pe=CreatePVector(element,grids,materials,inputs,analysis_type);
0901
0902
0903
0904 return;
0905
0906 else
0907
0908
0909
0910
0911
0912 numgrids=6;
0913 numgrids2d=3;
0914 NDOF2=2;
0915
0916
0917 pe=elemvector(numgrids*NDOF2);
0918
0919
0920 matice=materials(pentaelem.matid).material;
0921 matpar=materials(end).constants;
0922
0923
0924 gravity=matpar.g;
0925 rho_ice=matpar.rho_ice;
0926 rho_water=matpar.rho_water;
0927
0928
0929 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0930 [surface_param surface_is_present]=recover_input(inputs,'surface');
0931 [bed_param bed_is_present]=recover_input(inputs,'bed');
0932 [melting_param melting_is_present]=recover_input(inputs,'melting');
0933 [basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
0934
0935
0936 xyz_list=getgriddata(pentaelem,grids);
0937
0938
0939 basal_drag_list=zeros(numgrids,1);
0940 melting_list=zeros(numgrids,1);
0941 thickness_list=zeros(numgrids,1);
0942 surface_list=zeros(numgrids,1);
0943 bed_list=zeros(numgrids,1);
0944 temperature_list=zeros(numgrids,1);
0945
0946
0947 for i=1:numgrids,
0948
0949 doflist=grids(pentaelem.g(i)).grid.doflist;
0950 for j=1:NDOF2,
0951 pe.row_indices((i-1)*NDOF2+j)=doflist(j);
0952 end
0953 dof=doflist(1);
0954 if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
0955 if(surface_is_present) surface_list(i)=surface_param(dof);end;
0956 if(bed_is_present) bed_list(i)=bed_param(dof);end;
0957 if(melting_is_present) melting_list(i)=melting_param(dof);end;
0958 if(basal_drag_is_present) basal_drag_list(i)=basal_drag_param(dof);end;
0959
0960 end
0961
0962
0963
0964
0965
0966
0967
0968
0969
0970
0971
0972
0973
0974
0975
0976
0977
0978
0979
0980
0981
0982
0983
0984
0985
0986
0987
0988
0989
0990
0991
0992
0993
0994
0995
0996
0997
0998
0999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035 area_order=2;
1036 num_vert_gauss=3;
1037 [num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);
1038
1039
1040
1041
1042
1043
1044 for igarea=1:num_area_gauss,
1045 for igvert=1:num_vert_gauss,
1046
1047 gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
1048 gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
1049
1050
1051 if thickness_is_present,
1052 thickness=GetParameterValue(pentaelem,thickness_list,gauss_coord);
1053 else
1054 thickness=GetParameterValue(pentaelem,pentaelem.h,gauss_coord);
1055 end
1056
1057
1058 if surface_is_present,
1059 slope=GetParameterDerivativeValue(pentaelem,surface_list,xyz_list,gauss_coord);
1060 else
1061 slope=GetParameterDerivativeValue(pentaelem,pentaelem.s,xyz_list,gauss_coord);
1062 end
1063
1064
1065
1066
1067 if (~pentaelem.shelf) & (pentaelem.friction_type==1),
1068 if basal_drag_is_present,
1069 plastic_stress=GetParameterValue(pentaelem,basal_drag_list,gauss_coord);
1070 else
1071 plastic_stress=GetParameterValue(pentaelem,pentaelem.k,gauss_coord);
1072 end
1073 end
1074
1075
1076 Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
1077
1078
1079
1080 l1l6=GetNodalFunctions(pentaelem,gauss_coord);
1081
1082
1083 driving_stress_baseline=rho_ice*gravity;
1084
1085
1086 pe_gaussian=zeros(numgrids*NDOF2,1);
1087
1088 if(~pentaelem.shelf) & (pentaelem.friction_type==1),
1089 for i=1:numgrids,
1090 for j=1:NDOF2,
1091 pe_gaussian((i-1)*NDOF2+j)=(-driving_stress_baseline*slope(j)-plastic_stress)*Jdet*gauss_weight*l1l6(i);
1092 end
1093 end
1094 else
1095 for i=1:numgrids,
1096 for j=1:NDOF2,
1097 pe_gaussian((i-1)*NDOF2+j)=-driving_stress_baseline*slope(j)*Jdet*gauss_weight*l1l6(i);
1098 end
1099 end
1100 end
1101
1102
1103 for i=1:pe.nrows,
1104 pe.terms(i)=pe.terms(i)+pe_gaussian(i);
1105 end
1106
1107 end
1108 end
1109 end
1110 end
1111
1112
1113
1114
1115 function pe=CreatePVectorPrognostic(pentaelem,grids,materials,inputs);
1116
1117
1118
1119 if pentaelem.onbed~=1,
1120 pe=elemvector(0);
1121 return;
1122 end
1123
1124
1125 numgrids=3;
1126 DOFPERGRID=1;
1127 numdof=numgrids*DOFPERGRID;
1128
1129
1130 pe=elemvector(numdof);
1131
1132
1133 xyz_list=getgriddata(pentaelem,grids);
1134
1135
1136 xyz_list=xyz_list(1:3,:);
1137
1138
1139 [accumulation_param accumulation_is_present]=recover_input(inputs,'accumulation');
1140 [melting_param melting_is_present]=recover_input(inputs,'melting');
1141 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
1142 [dt dt_is_present]=recover_input(inputs,'dt');
1143
1144
1145 if ( ~accumulation_is_present | ~melting_is_present | ~thickness_is_present | ~dt_is_present),
1146 error('CreateKMatrixPrognostic error message: missing input parameters!');
1147 end
1148
1149
1150 thickness_list=zeros(numgrids,1);
1151 accumulation_list=zeros(numgrids,1);
1152 melting_list=zeros(numgrids,1);
1153
1154
1155 for i=1:numgrids,
1156 doflist=grids(pentaelem.g(i)).grid.doflist;
1157 dof=doflist(1);
1158
1159 pe.row_indices(i)=dof;
1160
1161
1162 thickness_list(i)=thickness_param(dof);
1163 melting_list(i)=melting_param(dof);
1164 accumulation_list(i)=accumulation_param(dof);
1165 end
1166
1167
1168 [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
1169
1170 for ig=1:num_gauss2D,
1171
1172
1173 gauss_weight=gauss_weights(ig);
1174 gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
1175
1176
1177 melting=GetParameterValue(triaelem,melting_list,gauss_coord);
1178
1179
1180 accumulation=GetParameterValue(triaelem,accumulation_list,gauss_coord);
1181
1182
1183 thickness=GetParameterValue(triaelem,thickness_list,gauss_coord);
1184
1185
1186 Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
1187
1188
1189 L=GetL(triaelem,gauss_coord,DOFPERGRID);
1190
1191
1192 pe_g_gaussian=Jdettria*gauss_weight*(thickness+dt*(accumulation-melting))*L';
1193
1194
1195 pe.terms=pe.terms+pe_g_gaussian;
1196 end
1197
1198 end