0001 function Ke=CreateKMatrix(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 Ke=CreateKMatrixHoriz(triaelem,grids,materials,inputs);
0015
0016 elseif strcmpi(analysis_type,'diagnostic_vert'),
0017
0018 Ke=CreateKMatrixVert(triaelem,grids,materials,inputs);
0019
0020 elseif strcmpi(analysis_type,'temperature'),
0021
0022 Ke=CreateKMatrixThermal(triaelem,grids,materials,inputs);
0023
0024 elseif strcmpi(analysis_type,'prognostic'),
0025
0026 Ke=CreateKMatrixPrognostic(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 Ke=CreateKMatrixSlopeCompute(triaelem,grids,materials,inputs);
0031
0032 end
0033 end
0034
0035
0036 function Ke=CreateKMatrixVert(triaelem,grids,materials,inputs)
0037 error('CreateKMatrix error message: Vertical analysis not implemented yet');
0038 end
0039
0040
0041 function Ke=CreateKMatrixPrognostic(triaelem,grids,materials,inputs)
0042
0043
0044 numgrids=3;
0045 DOFPERGRID=1;
0046 numdof=numgrids*DOFPERGRID;
0047
0048
0049 Ke=elemmatrix(numdof);
0050
0051
0052 xyz_list=getgriddata(triaelem,grids);
0053
0054
0055 [dt dt_is_present]=recover_input(inputs,'dt');
0056 [velocity_average_param velocity_average_is_present]=recover_input(inputs,'velocity_average');
0057
0058
0059 if ((~dt_is_present) | (~velocity_average_is_present)),
0060 error('CreateKMatrixPrognostic error message: missing input parameters!');
0061 end
0062
0063
0064 vxvy_list=zeros(numgrids,2);
0065
0066
0067 for i=1:numgrids,
0068 doflist=grids(triaelem.g(i)).grid.doflist;
0069 Ke.row_indices(i)=doflist(1);
0070 vxvy_list(i,:)=velocity_average_param(doflist(1:2));
0071 end
0072
0073
0074 if triaelem.artificial_diffusivity,
0075
0076 Jdettria=GetJacobianDeterminant(triaelem,xyz_list,[1/3 1/3 1/3]);
0077
0078
0079 v_gauss=1/3*(vxvy_list(1,:)+vxvy_list(2,:)+vxvy_list(3,:));
0080 K=Jdettria^(1/2)/2*[abs(v_gauss(1)) 0 ; 0 abs(v_gauss(2))];
0081 end
0082
0083
0084 [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0085
0086 for ig=1:num_gauss2D,
0087
0088
0089 gauss_weight=gauss_weights(ig);
0090 gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0091
0092
0093 Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
0094
0095
0096 L=GetL(triaelem,gauss_coord,DOFPERGRID);
0097
0098 DL_scalar=gauss_weight*Jdettria;
0099
0100
0101 Ke_gg_gaussian=L'*DL_scalar*L;
0102
0103
0104 B=GetB_prog(triaelem,xyz_list,gauss_coord);
0105
0106
0107 Bprime=GetBprime_prog(triaelem,xyz_list,gauss_coord);
0108
0109
0110 u_g=GetParameterValue(triaelem,vxvy_list(:,1),gauss_coord);
0111 v_g=GetParameterValue(triaelem,vxvy_list(:,2),gauss_coord);
0112 du_g=GetParameterDerivativeValue(triaelem,vxvy_list(:,1),xyz_list,gauss_coord);
0113 dv_g=GetParameterDerivativeValue(triaelem,vxvy_list(:,2),xyz_list,gauss_coord);
0114 dux_g=du_g(1);
0115 dvy_g=dv_g(2);
0116
0117 DL_scalar=dt*gauss_weight*Jdettria;
0118
0119
0120 DL=DL_scalar*[dux_g, 0; 0, dvy_g];
0121 DLprime=DL_scalar*[u_g, 0; 0, v_g];
0122
0123
0124 Ke_gg_thickness=B'*DL*B+B'*DLprime*Bprime;
0125
0126
0127 Ke.terms=Ke.terms+Ke_gg_gaussian+Ke_gg_thickness;
0128
0129
0130 if triaelem.artificial_diffusivity,
0131
0132
0133 Ke_gg_gaussian=Bprime'*DL_scalar*K*Bprime;
0134
0135
0136 Ke.terms=Ke.terms+Ke_gg_gaussian;
0137 end
0138 end
0139 end
0140
0141
0142 function Ke=CreateKMatrixThermal(triaelem,grids,materials,inputs)
0143 error('CreateKMatrix error message: Thermal analysis not implemented yet');
0144 end
0145
0146
0147 function Ke=CreateKMatrixSlopeCompute(triaelem,grids,materials,inputs)
0148
0149
0150 numgrids=3;
0151 DOFPERGRID=1;
0152 numdof=numgrids*DOFPERGRID;
0153
0154
0155 Ke=elemmatrix(numdof);
0156
0157
0158 xyz_list=getgriddata(triaelem,grids);
0159
0160
0161 xyz_list=xyz_list(1:3,:);
0162
0163
0164 for i=1:numgrids,
0165 doflist=grids(triaelem.g(i)).grid.doflist;
0166 dof=doflist(1);
0167 Ke.row_indices(i)=dof;
0168 end
0169
0170
0171 [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0172
0173 for ig=1:num_gauss2D,
0174
0175
0176 gauss_weight=gauss_weights(ig);
0177 gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0178
0179
0180 Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
0181
0182
0183 L=GetL(triaelem,gauss_coord,DOFPERGRID);
0184
0185 DL_scalar=gauss_weight*Jdettria;
0186
0187
0188 Ke_gg_gaussian=L'*DL_scalar*L;
0189
0190
0191 Ke.terms=Ke.terms+Ke_gg_gaussian;
0192 end
0193
0194 end
0195
0196
0197
0198 function Ke=CreateKMatrixHoriz(triaelem,grids,materials,inputs);
0199
0200 global element_debug element_debugid
0201
0202
0203
0204 if triaelem.acceleration==1,
0205 Ke=elemmatrix(0);
0206 return;
0207 end
0208
0209
0210 numgrids=3;
0211 DOFPERGRID=2;
0212 numdof=numgrids*DOFPERGRID;
0213
0214
0215 MAXSLOPE=.06;
0216 MOUNTAINKEXPONENT=10;
0217
0218
0219 Ke=elemmatrix(numdof);
0220
0221
0222 matice=materials(triaelem.matid).material;
0223 matpar=materials(end).constants;
0224
0225
0226 gravity=matpar.g;
0227 viscosity_overshoot=matpar.viscosity_overshoot;
0228 rho_ice=matpar.rho_ice;
0229 rho_water=matpar.rho_water;
0230
0231
0232 [basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
0233 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0234 [surface_param surface_is_present]=recover_input(inputs,'surface');
0235 [bed_param bed_is_present]=recover_input(inputs,'bed');
0236 [temperature_average_param temperature_average_is_present]=recover_input(inputs,'temperature_average');
0237 [flow_law_param flow_law_is_present]=recover_input(inputs,'B');
0238 [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0239 [oldvelocity_param oldvelocity_is_present]=recover_input(inputs,'oldvelocity');
0240
0241
0242 vxvy_list=zeros(numgrids,2);
0243 oldvxvy_list=zeros(numgrids,2);
0244
0245 B_list=zeros(numgrids,1);
0246 K_list=zeros(numgrids,1);
0247 thickness_list=zeros(numgrids,1);
0248 surface_list=zeros(numgrids,1);
0249 bed_list=zeros(numgrids,1);
0250 temperature_average_list=zeros(numgrids,1);
0251
0252
0253
0254 xyz_list=getgriddata(triaelem,grids);
0255
0256
0257
0258 for i=1:numgrids,
0259 doflist=grids(triaelem.g(i)).grid.doflist;
0260 for j=1:DOFPERGRID,
0261 dof=doflist(j);
0262 Ke.row_indices((i-1)*DOFPERGRID+j)=dof;
0263 if velocity_is_present, vxvy_list(i,j)=velocity_param(dof); end;
0264 if oldvelocity_is_present, oldvxvy_list(i,j)=oldvelocity_param(dof); end;
0265 end
0266
0267 dof=doflist(1);
0268 if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
0269 if(basal_drag_is_present),
0270 K_list(i)=basal_drag_param(dof);
0271 else
0272 K_list(i)=triaelem.k(i);
0273 end
0274 if(thickness_is_present),
0275 thickness_list(i)=thickness_param(dof);
0276 else
0277 thickness_list(i)=triaelem.h(i);
0278 end
0279 if(surface_is_present),
0280 surface_list(i)=surface_param(dof);
0281 else
0282 surface_list(i)=triaelem.s(i);
0283 end
0284 if(bed_is_present),
0285 bed_list(i)=bed_param(dof);
0286 else
0287 bed_list(i)=triaelem.b(i);
0288 end
0289 if(temperature_average_is_present) temperature_average_list(i)=temperature_average_param(dof);end;
0290 end
0291 if (element_debug & triaelem.id==element_debugid),
0292 disp(sprintf('El id %i TriaElemnet input list before gaussian loop: \n',element_debugid));
0293 disp(sprintf(' rho_ice: %g ',rho_ice));
0294 disp(sprintf(' gravity: %g',gravity));
0295 disp(sprintf(' rho_water: %g',rho_water));
0296 disp(sprintf(' Velocity: '));
0297 for i=1:3,
0298 disp(sprintf(' grid %i [%g,%g]',i,vxvy_list(i,1),vxvy_list(i,2)));
0299 end
0300 disp(sprintf(' B [%g %g %g ]',B_list(1),B_list(2),B_list(3)));
0301 disp(sprintf(' K [%g %g %g ]',K_list(1),K_list(2),K_list(3)));
0302 disp(sprintf(' thickness [%g %g %g]',thickness_list(1),thickness_list(2),thickness_list(3)));
0303 disp(sprintf(' surface [%g %g %g ]',surface_list(1),surface_list(2),surface_list(3)));
0304 disp(sprintf(' bed [%g %g %g]',bed_list(1),bed_list(2),bed_list(3)));
0305 if(temperature_average_is_present)disp(sprintf(' temperature_average [%g %g %g]',temperature_average_list(1),temperature_average_list(2),temperature_average_list(3)));end;
0306 end
0307
0308
0309
0310 alpha2_list=zeros(3,1);
0311
0312 if (~triaelem.shelf & triaelem.friction_type==2)
0313
0314
0315 frictionparameters=struct();
0316 frictionparameters.element_type='2d';
0317 frictionparameters.rho_ice=rho_ice;
0318 frictionparameters.rho_water=rho_water;
0319 frictionparameters.g=gravity;
0320 frictionparameters.p=triaelem.p;
0321 frictionparameters.q=triaelem.q;
0322
0323 if velocity_is_present,
0324 frictionparameters.velocities=vxvy_list;
0325 else
0326 frictionparameters.velocities=zeros(3,2);
0327 end
0328 frictionparameters.h=thickness_list;
0329 frictionparameters.b=bed_list;
0330 frictionparameters.k=K_list;
0331
0332 alpha2=Getalpha2(frictionparameters);
0333 if (element_debug & triaelem.id==element_debugid),
0334 disp(sprintf(' alpha2_list (%g %g %g )',alpha2_list(1),alpha2_list(2),alpha2_list(3)));
0335 end
0336 end
0337
0338
0339
0340
0341 [num_gauss,first_gauss_area_coord,second_gauss_area_coord,third_gauss_area_coord,gauss_weights]=GaussTria(2);
0342
0343 if (element_debug & triaelem.id==element_debugid),
0344 disp(sprintf(' gaussian points: '));
0345 for i=1:num_gauss,
0346 disp(sprintf(' %g %g %g : %g',first_gauss_area_coord(i),second_gauss_area_coord(i),third_gauss_area_coord(i),gauss_weights(i)));
0347 end
0348 end
0349
0350
0351 for ig=1:num_gauss,
0352
0353 gauss_weight=gauss_weights(ig);
0354 gauss_l1l2l3=[first_gauss_area_coord(ig) second_gauss_area_coord(ig) third_gauss_area_coord(ig)];
0355
0356
0357 if thickness_is_present,
0358 thickness=GetParameterValue(triaelem,thickness_list,gauss_l1l2l3);
0359 else
0360 thickness=GetParameterValue(triaelem,triaelem.h,gauss_l1l2l3);
0361 end
0362
0363
0364
0365
0366 slopevector=zeros(2,1);
0367 if(~triaelem.shelf),
0368
0369 if surface_is_present,
0370 slopevector=GetParameterDerivativeValue(triaelem,surface_list,xyz_list,gauss_l1l2l3);
0371 else
0372 slopevector=GetParameterDerivativeValue(triaelem,triaelem.s,xyz_list,gauss_l1l2l3);
0373 end
0374 slope_magnitude=sqrt(slopevector(1)^2+slopevector(2)^2);
0375
0376 if (slope_magnitude>MAXSLOPE),
0377 K=10^(-MOUNTAINKEXPONENT);
0378 end
0379 end
0380
0381
0382 if velocity_is_present,
0383 epsilon=GetStrainRate(triaelem,vxvy_list,xyz_list,gauss_l1l2l3);
0384 else
0385 epsilon=zeros(3,1);
0386 end
0387 if oldvelocity_is_present,
0388 oldepsilon=GetStrainRate(triaelem,oldvxvy_list,xyz_list,gauss_l1l2l3);
0389 else
0390 oldepsilon=zeros(3,1);
0391 end
0392
0393
0394
0395
0396 if temperature_average_is_present,
0397 temperature_average=GetParameterValue(triaelem,temperature_average_list,gauss_l1l2l3);
0398 matice.B=paterson(temperature_average);
0399 end
0400
0401
0402
0403 if flow_law_is_present,
0404 B_param=GetParameterValue(triaelem,B_list,gauss_l1l2l3);
0405 matice.B=B_param; clear B_param.
0406 end
0407
0408
0409
0410 newviscosity=GetViscosity2d(matice,epsilon);
0411 oldviscosity=GetViscosity2d(matice,oldepsilon);
0412 viscosity=newviscosity+viscosity_overshoot*(newviscosity-oldviscosity);
0413
0414
0415 Jdet=GetJacobianDeterminant(triaelem,xyz_list,gauss_l1l2l3);
0416
0417
0418
0419 D=viscosity*thickness*gauss_weight*Jdet*diag(ones(numgrids,1));
0420
0421 if (element_debug & triaelem.id==element_debugid),
0422 disp(sprintf(' gaussian loop %i\n',ig));
0423 disp(sprintf(' thickness %g',thickness));
0424 disp(sprintf(' slope (%g,%g)',slopevector(1),slopevector(2)));
0425 disp(sprintf(' alpha2_list (%g,%g,%g)',alpha2_list(1),alpha2_list(2),alpha2_list(3)));
0426 disp(sprintf(' epsilon (%g,%g,%g)',epsilon(1),epsilon(2),epsilon(3)));
0427 disp(sprintf(' Matice: '));
0428 matice
0429 disp(sprintf('\n viscosity: %g ',viscosity));
0430 disp(sprintf(' jacobian: %g ',Jdet));
0431 disp(sprintf(' gauss_weight: %g ',gauss_weight));
0432 end
0433
0434
0435
0436 B=GetB(triaelem,xyz_list,gauss_l1l2l3);
0437 Bprime=GetBprime(triaelem,xyz_list,gauss_l1l2l3);
0438
0439
0440 L=zeros(2,DOFPERGRID*numgrids);
0441 if (triaelem.friction_type==2 && triaelem.shelf==0),
0442 L=GetL(triaelem,gauss_l1l2l3,DOFPERGRID);
0443 end
0444
0445
0446
0447 Ke_gg_gaussian=B'*D*Bprime;
0448
0449
0450 Ke.terms=Ke.terms+Ke_gg_gaussian;
0451
0452
0453 alpha2_g=0;
0454 if (~triaelem.shelf) & (triaelem.friction_type==2),
0455
0456
0457 alpha2_g=GetParameterValue(triaelem,alpha2,gauss_l1l2l3);
0458
0459 if velocity_is_present
0460 DL_scalar=alpha2_g*gauss_weight*Jdet;
0461 else
0462 DL_scalar=0;
0463 end
0464
0465 DL=diag(DL_scalar*ones(2,1));
0466
0467
0468 Ke_gg_drag_gaussian=L'*DL*L;
0469
0470
0471 Ke.terms=Ke.terms+Ke_gg_drag_gaussian;
0472 end
0473
0474 if (element_debug & triaelem.id==element_debugid),
0475 disp(sprintf(' alpha2 %g\n',alpha2_g));
0476 disp(sprintf('B:\n'));
0477 B
0478 disp(sprintf('Bprime:\n'));
0479 Bprime
0480 disp(sprintf('L:\n'));
0481 L
0482 end
0483
0484 end
0485
0486 if (element_debug & triaelem.id==element_debugid),
0487 disp(sprintf('Ke_gg->terms:\n'));
0488 Ke.terms
0489 disp(sprintf('Ke_gg->row_indices:\n'));
0490 Ke.row_indices'
0491 end
0492
0493 end
0494