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