0001 function pe=CreatePVector(load,elements,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(load,elements,grids,materials,inputs);
0015
0016 elseif strcmpi(analysis_type,'diagnostic_stokes'),
0017
0018 pe=CreatePVectorStokes(load,elements,grids,materials,inputs);
0019
0020 else
0021 error('CreatePVector error message: only horizontal diagnostic should use icefront objects');
0022
0023 end
0024 end
0025
0026 function pe=CreatePVectorStokes(load,elements,grids,materials,inputs);
0027
0028
0029 numgrids=4;
0030 NDOF4=4;
0031
0032 numgridstria=4;
0033 NDOF4=4;
0034
0035
0036
0037 pe=elemvector(numgrids*NDOF4);
0038
0039
0040 pentaelem=elements(load.eid).element;
0041 matpar=materials(end).constants;
0042
0043
0044 gravity=matpar.g;
0045 rho_ice=matpar.rho_ice;
0046 rho_water=matpar.rho_water;
0047
0048
0049 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0050 [bed_param bed_is_present]=recover_input(inputs,'bed');
0051
0052
0053 thickness_list=zeros(6,1);
0054 bed_list=zeros(6,1);
0055
0056
0057 xyz_list=getgriddata(pentaelem,grids);
0058
0059
0060 for i=1:6,
0061 if pentaelem.g(i)==load.g(1),
0062 grid1=i;
0063 end
0064 if pentaelem.g(i)==load.g(2),
0065 grid2=i;
0066 end
0067 if pentaelem.g(i)==load.g(3),
0068 grid3=i;
0069 end
0070 if pentaelem.g(i)==load.g(4),
0071 grid4=i;
0072 end
0073 end
0074 quad_grids=[grid1 grid2 grid3 grid4];
0075
0076
0077 for i=1:numgrids,
0078 doflist=grids(load.g(i)).grid.doflist;
0079 for j=1:NDOF4,
0080 pe.row_indices((i-1)*NDOF4+j)=doflist(j);
0081 end
0082 dof=doflist(1);
0083 end
0084
0085 for i=1:6,
0086 doflist=grids(pentaelem.g(i)).grid.doflist; dof=doflist(1);
0087 if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
0088 if(bed_is_present) bed_list(i)=bed_param(dof);end;
0089 end
0090
0091
0092
0093 if thickness_is_present,
0094 thickness=thickness_list(quad_grids);
0095 bed=bed_list(quad_grids);
0096 else
0097 thickness=pentaelem.h(quad_grids);
0098 bed=pentaelem.b(quad_grids);
0099 end
0100
0101
0102 xyz_list_quad=xyz_list(quad_grids,:);
0103
0104
0105 x5=mean(xyz_list_quad(:,1));
0106 y5=mean(xyz_list_quad(:,2));
0107 z5=mean(xyz_list_quad(:,3));
0108
0109 xyz_list_quad=[xyz_list_quad; x5 y5 z5];
0110
0111
0112 V1=cross(xyz_list_quad(1,:)-xyz_list_quad(5,:),xyz_list_quad(2,:)-xyz_list_quad(5,:));
0113 normal1=1/norm(V1)*V1';
0114
0115 V2=cross(xyz_list_quad(2,:)-xyz_list_quad(5,:),xyz_list_quad(3,:)-xyz_list_quad(5,:));
0116 normal2=1/norm(V2)*V2';
0117
0118 V3=cross(xyz_list_quad(3,:)-xyz_list_quad(5,:),xyz_list_quad(4,:)-xyz_list_quad(5,:));
0119 normal3=1/norm(V3)*V3';
0120
0121 V4=cross(xyz_list_quad(4,:)-xyz_list_quad(5,:),xyz_list_quad(1,:)-xyz_list_quad(5,:));
0122 normal4=1/norm(V4)*V4';
0123
0124
0125 normal=[normal1 normal2 normal3 normal4];
0126
0127
0128 if pentaelem.shelf==1,
0129
0130 pe.terms=QuadPressureLoadStokes(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,1);
0131 else
0132
0133 pe.terms=QuadPressureLoadStokes(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,0);
0134 end
0135
0136 end
0137
0138 function pe=QuadPressureLoadStokes(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,fill);
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164 nx=normal(1,:);
0165 ny=normal(2,:);
0166 nz=normal(3,:);
0167
0168
0169 [num_gauss,first_gauss_area,second_gauss_area,third_gauss_area,gauss_weights]=GaussTria(2);
0170
0171
0172 s=thickness+bed;
0173
0174
0175 s(5)=mean(s);
0176
0177
0178 x=xyz_list_quad(:,1);
0179 y=xyz_list_quad(:,2);
0180 z=xyz_list_quad(:,3);
0181
0182
0183
0184
0185 tria=triaelem;
0186
0187 l12=(x(2)-x(1))^2+(y(2)-y(1))^2+(z(2)-z(1))^2;
0188 l14=(x(4)-x(1))^2+(y(4)-y(1))^2+(z(4)-z(1))^2;
0189 l15=(x(5)-x(1))^2+(y(5)-y(1))^2+(z(5)-z(1))^2;
0190 l23=(x(3)-x(2))^2+(y(3)-y(2))^2+(z(3)-z(2))^2;
0191 l25=(x(5)-x(2))^2+(y(5)-y(2))^2+(z(5)-z(2))^2;
0192 l34=(x(4)-x(3))^2+(y(4)-y(3))^2+(z(4)-z(3))^2;
0193 l35=(x(5)-x(3))^2+(y(5)-y(3))^2+(z(5)-z(3))^2;
0194 l45=(x(5)-x(4))^2+(y(5)-y(4))^2+(z(5)-z(4))^2;
0195 cos_theta_triangle1=(l15+l12-l25)/(2*sqrt(l12*l15));
0196 cos_theta_triangle2=(l25+l23-l35)/(2*sqrt(l23*l25));
0197 cos_theta_triangle3=(l35+l34-l45)/(2*sqrt(l34*l35));
0198 cos_theta_triangle4=(l45+l14-l15)/(2*sqrt(l14*l45));
0199
0200
0201 x1tria1=0;
0202 y1tria1=0;
0203 x2tria1=sqrt(l12);
0204 y2tria1=0;
0205 x3tria1=cos_theta_triangle1*sqrt(l15);
0206 y3tria1=sqrt(1-cos_theta_triangle1^2)*sqrt(l15);
0207
0208 xyz_tria1=[x1tria1 y1tria1; x2tria1 y2tria1; x3tria1 y3tria1];
0209
0210
0211 x1tria2=0;
0212 y1tria2=0;
0213 x2tria2=sqrt(l23);
0214 y2tria2=0;
0215 x3tria2=cos_theta_triangle2*sqrt(l25);
0216 y3tria2=sqrt(1-cos_theta_triangle2^2)*sqrt(l25);
0217
0218 xyz_tria2=[x1tria2 y1tria2; x2tria2 y2tria2; x3tria2 y3tria2];
0219
0220
0221 x1tria3=0;
0222 y1tria3=0;
0223 x2tria3=sqrt(l34);
0224 y2tria3=0;
0225 x3tria3=cos_theta_triangle3*sqrt(l35);
0226 y3tria3=sqrt(1-cos_theta_triangle3^2)*sqrt(l35);
0227
0228 xyz_tria3=[x1tria3 y1tria3; x2tria3 y2tria3; x3tria3 y3tria3];
0229
0230
0231 x1tria4=0;
0232 y1tria4=0;
0233 x2tria4=sqrt(l14);
0234 y2tria4=0;
0235 x3tria4=cos_theta_triangle4*sqrt(l45);
0236 y3tria4=sqrt(1-cos_theta_triangle4^2)*sqrt(l45);
0237
0238 xyz_tria4=[x1tria4 y1tria4; x2tria4 y2tria4; x3tria4 y3tria4];
0239
0240
0241
0242 xyz_tria=[xyz_tria1; xyz_tria2; xyz_tria3; xyz_tria4];
0243
0244
0245 pe=zeros(16,1);
0246
0247
0248 for ig=1:num_gauss,
0249
0250
0251 gauss_coord=[first_gauss_area(ig) second_gauss_area(ig) third_gauss_area(ig)];
0252 l1l2l3_tria=GetNodalFunctions(tria,gauss_coord);
0253
0254
0255 r_tria=gauss_coord(2)-gauss_coord(1);
0256 s_tria=-3/sqrt(3)*(gauss_coord(1)+gauss_coord(2)-2/3);
0257
0258
0259 r_quad(1,1)=r_tria;
0260 s_quad(1,1)=1/sqrt(3)*s_tria-2/3;
0261 r_quad(2,1)=-1/sqrt(3)*s_tria+2/3;
0262 s_quad(2,1)=r_tria;
0263 r_quad(3,1)=-r_tria;
0264 s_quad(3,1)=-1/sqrt(3)*s_tria+2/3;
0265 r_quad(4,1)=1/sqrt(3)*s_tria-2/3;
0266 s_quad(4,1)=-r_tria;
0267
0268
0269 l1l4_tria=1/4*[(r_quad-1).*(s_quad-1) -(r_quad+1).*(s_quad-1) (r_quad+1).*(s_quad+1) -(r_quad-1).*(s_quad+1)];
0270
0271
0272 for i=1:4
0273 complete_list=zeros(3,3);
0274 complete_list(1:3,1:2)=xyz_tria(3*i-2:3*i,:);
0275 J(i)=GetJacobianDeterminant(tria,complete_list,l1l2l3_tria);
0276 end
0277
0278
0279 z_g(1)=z(1)*l1l2l3_tria(1)+z(2)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0280 z_g(2)=z(2)*l1l2l3_tria(1)+z(3)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0281 z_g(3)=z(3)*l1l2l3_tria(1)+z(4)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0282 z_g(4)=z(4)*l1l2l3_tria(1)+z(1)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0283
0284
0285 for i=1:4
0286
0287
0288 air_pressure_tria=0;
0289
0290
0291 if fill==1,
0292 water_level_above_g_tria=min(0,z_g(i));
0293 water_pressure_tria=rho_water*gravity*water_level_above_g_tria;
0294 elseif fill==0,
0295 water_pressure_tria=0;
0296 else
0297 error('QuadPressureLoad error message: unknow fill type for icefront boundary condition');
0298 end
0299
0300
0301
0302 for j=1:4
0303 pressure_tria(j) = water_pressure_tria + air_pressure_tria;
0304 end
0305
0306 pe=pe+J(i)*gauss_weights(ig)*...
0307 [ pressure_tria(1)*l1l4_tria(i,1)*nx(i)
0308 pressure_tria(1)*l1l4_tria(i,1)*ny(i)
0309 pressure_tria(1)*l1l4_tria(i,1)*nz(i)
0310 0
0311 pressure_tria(2)*l1l4_tria(i,2)*nx(i)
0312 pressure_tria(2)*l1l4_tria(i,2)*ny(i)
0313 pressure_tria(2)*l1l4_tria(i,2)*nz(i)
0314 0
0315 pressure_tria(3)*l1l4_tria(i,3)*nx(i)
0316 pressure_tria(3)*l1l4_tria(i,3)*ny(i)
0317 pressure_tria(3)*l1l4_tria(i,3)*nz(i)
0318 0
0319 pressure_tria(4)*l1l4_tria(i,4)*nx(i)
0320 pressure_tria(4)*l1l4_tria(i,4)*ny(i)
0321 pressure_tria(4)*l1l4_tria(i,4)*nz(i)
0322 0 ];
0323
0324 end
0325
0326 end
0327
0328 end
0329
0330 function pe=CreatePVectorHoriz(load,elements,grids,materials,inputs);
0331
0332 global element_debug element_debugid
0333
0334
0335 if strcmpi(load.type,'segment'),
0336
0337 if strcmpi(elements(load.eid).element.type,'pentaelem'),
0338 if ~elements(load.eid).element.onbed
0339 pe=elemvector(0);
0340 return;
0341 end
0342 end
0343
0344
0345 numgrids=2;
0346 NDOF2=2;
0347
0348
0349 pe=elemvector(numgrids*NDOF2);
0350
0351
0352 triaelem=elements(load.eid).element;
0353 matpar=materials(end).constants;
0354
0355 if (element_debug & triaelem.id==element_debugid),
0356 disp(sprintf('%s%i','Inside icefront load connected to element ',triaelem.id));
0357 end
0358
0359
0360 gravity=matpar.g;
0361 rho_ice=matpar.rho_ice;
0362 rho_water=matpar.rho_water;
0363
0364
0365 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0366 [bed_param bed_is_present]=recover_input(inputs,'bed');
0367
0368
0369 xyz_list=getgriddata(triaelem,grids);
0370
0371
0372 thickness_list=zeros(3,1);
0373 bed_list=zeros(3,1);
0374
0375
0376 for i=1:3,
0377 if triaelem.g(i)==load.g(1),
0378 grid1=i;
0379 end
0380 if triaelem.g(i)==load.g(2),
0381 grid2=i;
0382 end
0383 end
0384
0385
0386 for i=1:numgrids,
0387 doflist=grids(load.g(i)).grid.doflist;
0388 for j=1:NDOF2,
0389 pe.row_indices((i-1)*NDOF2+j)=doflist(j);
0390 end
0391 end
0392
0393 for i=1:3,
0394 doflist=grids(triaelem.g(i)).grid.doflist; dof=doflist(1);
0395 if(thickness_is_present),
0396 thickness_list(i)=thickness_param(dof);
0397 else
0398 thickness_list(i)=triaelem.h(i);
0399 end
0400 if(bed_is_present),
0401 bed_list(i)=bed_param(dof);
0402 else
0403 bed_list(i)=triaelem.b(i);
0404 end;
0405 end
0406
0407
0408 if thickness_is_present,
0409 thickness=[thickness_list(grid1) thickness_list(grid2)];
0410 bed=[bed_list(grid1) bed_list(grid2)];
0411 else
0412 thickness=[triaelem.h(grid1) triaelem.h(grid2)];
0413 bed=[triaelem.b(grid1) triaelem.b(grid2)];
0414 end
0415
0416
0417 x1=xyz_list(grid1,1);
0418 y1=xyz_list(grid1,2);
0419 x2=xyz_list(grid2,1);
0420 y2=xyz_list(grid2,2);
0421
0422
0423 normal=zeros(2,1);
0424 normal(1)=cos(atan2(x1-x2,y2-y1));
0425 normal(2)=sin(atan2(x1-x2,y2-y1));
0426 length=sqrt((x2-x1)^2+(y2-y1)^2);
0427
0428
0429 if triaelem.shelf==1,
0430
0431 fill=1;
0432 else
0433
0434 fill=0;
0435 end
0436 pe.terms=SegmentPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,length,fill);
0437
0438 if (element_debug & triaelem.id==element_debugid),
0439
0440 disp(sprintf('\nicefront load'));
0441 disp(sprintf('grids %i %i',grid1,grid2));
0442 disp(sprintf('rho_water %g',rho_water));
0443 disp(sprintf('rho_ice %g',rho_ice));
0444 disp(sprintf('gravity %g',gravity));
0445 disp(sprintf('thickness (%g,%g)',thickness(1),thickness(2)));
0446 disp(sprintf('bed (%g,%g)',bed(1),bed(2)));
0447 disp(sprintf('normal (%g,%g)',normal(1),normal(2)));
0448 disp(sprintf('length %g',length));
0449 disp(sprintf('fill %i',fill));
0450
0451 disp(sprintf(' pe_g->terms\n'));
0452 pe.terms
0453 end
0454
0455 elseif strcmpi(load.type,'quad'),
0456
0457
0458 numgrids=4;
0459 NDOF2=2;
0460
0461 numgridstria=4;
0462 NDOF2=2;
0463
0464
0465
0466 pe=elemvector(numgrids*NDOF2);
0467
0468
0469 pentaelem=elements(load.eid).element;
0470 matpar=materials(end).constants;
0471
0472
0473 gravity=matpar.g;
0474 rho_ice=matpar.rho_ice;
0475 rho_water=matpar.rho_water;
0476
0477
0478 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0479 [bed_param bed_is_present]=recover_input(inputs,'bed');
0480
0481
0482 thickness_list=zeros(6,1);
0483 bed_list=zeros(6,1);
0484
0485
0486 xyz_list=getgriddata(pentaelem,grids);
0487
0488
0489 for i=1:6,
0490 if pentaelem.g(i)==load.g(1),
0491 grid1=i;
0492 end
0493 if pentaelem.g(i)==load.g(2),
0494 grid2=i;
0495 end
0496 if pentaelem.g(i)==load.g(3),
0497 grid3=i;
0498 end
0499 if pentaelem.g(i)==load.g(4),
0500 grid4=i;
0501 end
0502 end
0503 quad_grids=[grid1 grid2 grid3 grid4];
0504
0505
0506 for i=1:numgrids,
0507 doflist=grids(load.g(i)).grid.doflist;
0508 for j=1:NDOF2,
0509 pe.row_indices((i-1)*NDOF2+j)=doflist(j);
0510 end
0511 dof=doflist(1);
0512 end
0513
0514 for i=1:6,
0515 doflist=grids(pentaelem.g(i)).grid.doflist; dof=doflist(1);
0516 if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
0517 if(bed_is_present) bed_list(i)=bed_param(dof);end;
0518 end
0519
0520
0521
0522 if thickness_is_present,
0523 thickness=thickness_list(quad_grids);
0524 bed=bed_list(quad_grids);
0525 else
0526 thickness=pentaelem.h(quad_grids);
0527 bed=pentaelem.b(quad_grids);
0528 end
0529
0530
0531 xyz_list_quad=xyz_list(quad_grids,:);
0532
0533
0534 x5=mean(xyz_list_quad(:,1));
0535 y5=mean(xyz_list_quad(:,2));
0536 z5=mean(xyz_list_quad(:,3));
0537
0538 xyz_list_quad=[xyz_list_quad; x5 y5 z5];
0539
0540
0541 V1=cross(xyz_list_quad(1,:)-xyz_list_quad(5,:),xyz_list_quad(2,:)-xyz_list_quad(5,:));
0542 normal1=1/norm(V1)*V1';
0543
0544 V2=cross(xyz_list_quad(2,:)-xyz_list_quad(5,:),xyz_list_quad(3,:)-xyz_list_quad(5,:));
0545 normal2=1/norm(V2)*V2';
0546
0547 V3=cross(xyz_list_quad(3,:)-xyz_list_quad(5,:),xyz_list_quad(4,:)-xyz_list_quad(5,:));
0548 normal3=1/norm(V3)*V3';
0549
0550 V4=cross(xyz_list_quad(4,:)-xyz_list_quad(5,:),xyz_list_quad(1,:)-xyz_list_quad(5,:));
0551 normal4=1/norm(V4)*V4';
0552
0553
0554 normal=[normal1 normal2 normal3 normal4];
0555
0556
0557 if pentaelem.shelf==1,
0558
0559 pe.terms=QuadPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,1);
0560 else
0561
0562 pe.terms=QuadPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,0);
0563 end
0564
0565 if (element_debug & pentaelem.id==element_debugid),
0566 disp(sprintf(' pe_g->terms\n'));
0567 pe.terms
0568 end
0569
0570 else
0571 error('unsupported element type')
0572 end
0573
0574 end
0575
0576 function pe=QuadPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,fill);
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586
0587
0588
0589
0590
0591
0592
0593
0594
0595
0596
0597
0598
0599
0600
0601
0602 nx=normal(1,:);
0603 ny=normal(2,:);
0604
0605
0606 [num_gauss,first_gauss_area,second_gauss_area,third_gauss_area,gauss_weights]=GaussTria(2);
0607
0608
0609 s=thickness+bed;
0610
0611
0612 s(5)=mean(s);
0613
0614
0615 x=xyz_list_quad(:,1);
0616 y=xyz_list_quad(:,2);
0617 z=xyz_list_quad(:,3);
0618
0619
0620
0621
0622 tria=triaelem;
0623
0624 l12=(x(2)-x(1))^2+(y(2)-y(1))^2+(z(2)-z(1))^2;
0625 l14=(x(4)-x(1))^2+(y(4)-y(1))^2+(z(4)-z(1))^2;
0626 l15=(x(5)-x(1))^2+(y(5)-y(1))^2+(z(5)-z(1))^2;
0627 l23=(x(3)-x(2))^2+(y(3)-y(2))^2+(z(3)-z(2))^2;
0628 l25=(x(5)-x(2))^2+(y(5)-y(2))^2+(z(5)-z(2))^2;
0629 l34=(x(4)-x(3))^2+(y(4)-y(3))^2+(z(4)-z(3))^2;
0630 l35=(x(5)-x(3))^2+(y(5)-y(3))^2+(z(5)-z(3))^2;
0631 l45=(x(5)-x(4))^2+(y(5)-y(4))^2+(z(5)-z(4))^2;
0632 cos_theta_triangle1=(l15+l12-l25)/(2*sqrt(l12*l15));
0633 cos_theta_triangle2=(l25+l23-l35)/(2*sqrt(l23*l25));
0634 cos_theta_triangle3=(l35+l34-l45)/(2*sqrt(l34*l35));
0635 cos_theta_triangle4=(l45+l14-l15)/(2*sqrt(l14*l45));
0636
0637
0638 x1tria1=0;
0639 y1tria1=0;
0640 x2tria1=sqrt(l12);
0641 y2tria1=0;
0642 x3tria1=cos_theta_triangle1*sqrt(l15);
0643 y3tria1=sqrt(1-cos_theta_triangle1^2)*sqrt(l15);
0644
0645 xyz_tria1=[x1tria1 y1tria1; x2tria1 y2tria1; x3tria1 y3tria1];
0646
0647
0648 x1tria2=0;
0649 y1tria2=0;
0650 x2tria2=sqrt(l23);
0651 y2tria2=0;
0652 x3tria2=cos_theta_triangle2*sqrt(l25);
0653 y3tria2=sqrt(1-cos_theta_triangle2^2)*sqrt(l25);
0654
0655 xyz_tria2=[x1tria2 y1tria2; x2tria2 y2tria2; x3tria2 y3tria2];
0656
0657
0658 x1tria3=0;
0659 y1tria3=0;
0660 x2tria3=sqrt(l34);
0661 y2tria3=0;
0662 x3tria3=cos_theta_triangle3*sqrt(l35);
0663 y3tria3=sqrt(1-cos_theta_triangle3^2)*sqrt(l35);
0664
0665 xyz_tria3=[x1tria3 y1tria3; x2tria3 y2tria3; x3tria3 y3tria3];
0666
0667
0668 x1tria4=0;
0669 y1tria4=0;
0670 x2tria4=sqrt(l14);
0671 y2tria4=0;
0672 x3tria4=cos_theta_triangle4*sqrt(l45);
0673 y3tria4=sqrt(1-cos_theta_triangle4^2)*sqrt(l45);
0674
0675 xyz_tria4=[x1tria4 y1tria4; x2tria4 y2tria4; x3tria4 y3tria4];
0676
0677
0678
0679 xyz_tria=[xyz_tria1; xyz_tria2; xyz_tria3; xyz_tria4];
0680
0681
0682 pe=zeros(8,1);
0683
0684
0685 for ig=1:num_gauss,
0686
0687
0688 gauss_coord=[first_gauss_area(ig) second_gauss_area(ig) third_gauss_area(ig)];
0689 l1l2l3_tria=GetNodalFunctions(tria,gauss_coord);
0690
0691
0692 r_tria=gauss_coord(2)-gauss_coord(1);
0693 s_tria=-3/sqrt(3)*(gauss_coord(1)+gauss_coord(2)-2/3);
0694
0695
0696 r_quad(1,1)=r_tria;
0697 s_quad(1,1)=1/sqrt(3)*s_tria-2/3;
0698 r_quad(2,1)=-1/sqrt(3)*s_tria+2/3;
0699 s_quad(2,1)=r_tria;
0700 r_quad(3,1)=-r_tria;
0701 s_quad(3,1)=-1/sqrt(3)*s_tria+2/3;
0702 r_quad(4,1)=1/sqrt(3)*s_tria-2/3;
0703 s_quad(4,1)=-r_tria;
0704
0705
0706 l1l4_tria=1/4*[(r_quad-1).*(s_quad-1) -(r_quad+1).*(s_quad-1) (r_quad+1).*(s_quad+1) -(r_quad-1).*(s_quad+1)];
0707
0708
0709 for i=1:4
0710 complete_list=zeros(3,3);
0711 complete_list(1:3,1:2)=xyz_tria(3*i-2:3*i,:);
0712 J(i)=GetJacobianDeterminant(tria,complete_list,l1l2l3_tria);
0713 end
0714
0715
0716 z_g(1)=z(1)*l1l2l3_tria(1)+z(2)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0717 z_g(2)=z(2)*l1l2l3_tria(1)+z(3)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0718 z_g(3)=z(3)*l1l2l3_tria(1)+z(4)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0719 z_g(4)=z(4)*l1l2l3_tria(1)+z(1)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0720
0721
0722 for i=1:4
0723
0724
0725
0726 for j=1:4
0727 ice_pressure_tria(j)=gravity*rho_ice*(s(j)-z_g(i));
0728 end
0729 air_pressure_tria=0;
0730
0731
0732 if fill==1,
0733 water_level_above_g_tria=min(0,z_g(i));
0734 water_pressure_tria=rho_water*gravity*water_level_above_g_tria;
0735
0736 elseif fill==0,
0737 water_pressure_tria=0;
0738
0739 else
0740 error('QuadPressureLoad error message: unknow fill type for icefront boundary condition');
0741 end
0742
0743
0744
0745 for j=1:4
0746 pressure_tria(j) = ice_pressure_tria(j) + water_pressure_tria + air_pressure_tria;
0747 end
0748
0749 pe=pe+J(i)*gauss_weights(ig)*...
0750 [ pressure_tria(1)*l1l4_tria(i,1)*nx(i)
0751 pressure_tria(1)*l1l4_tria(i,1)*ny(i)
0752 pressure_tria(2)*l1l4_tria(i,2)*nx(i)
0753 pressure_tria(2)*l1l4_tria(i,2)*ny(i)
0754 pressure_tria(3)*l1l4_tria(i,3)*nx(i)
0755 pressure_tria(3)*l1l4_tria(i,3)*ny(i)
0756 pressure_tria(4)*l1l4_tria(i,4)*nx(i)
0757 pressure_tria(4)*l1l4_tria(i,4)*ny(i)];
0758
0759 end
0760
0761 end
0762
0763 end
0764
0765
0766
0767 function pe=SegmentPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,length,fill);
0768
0769
0770 nx=normal(1,1);
0771 ny=normal(2,1);
0772
0773
0774 num_gauss=3;
0775 [segment_gauss_coord, gauss_weights]=GaussSegment(num_gauss);
0776
0777
0778 h1=thickness(1);
0779 b1=bed(1);
0780 h2=thickness(2);
0781 b2=bed(2);
0782
0783
0784 Jdet=1/2*length;
0785
0786
0787 pe=zeros(4,1);
0788
0789
0790 for ig=1:num_gauss,
0791
0792 thickness=h1*(1+segment_gauss_coord(ig))/2+h2*(1-segment_gauss_coord(ig))/2;
0793 bed=b1*(1+segment_gauss_coord(ig))/2+b2*(1-segment_gauss_coord(ig))/2;
0794
0795 if fill==1,
0796 ice_pressure=1.0/2.0*gravity*rho_ice*thickness^2;
0797 air_pressure=0;
0798
0799
0800 surface_under_water=min(0,thickness+bed);
0801 base_under_water=min(0,bed);
0802 water_pressure=1.0/2.0*gravity*rho_water*(surface_under_water^2 - base_under_water^2);
0803
0804 elseif fill==0,
0805 ice_pressure=1.0/2.0*gravity*rho_ice*thickness^2;
0806 air_pressure=0;
0807 water_pressure=0;
0808
0809 else
0810 error('SegmentPressureLoad error message: unknow fill type for icefront boundary condition');
0811 end
0812
0813 pressure = ice_pressure + water_pressure + air_pressure;
0814
0815
0816 pe=pe+pressure*Jdet*gauss_weights(ig)*...
0817 [ (1+segment_gauss_coord(ig))/2*nx
0818 (1+segment_gauss_coord(ig))/2*ny
0819 (1-segment_gauss_coord(ig))/2*nx
0820 (1-segment_gauss_coord(ig))/2*ny];
0821
0822 end
0823
0824 end