0001 function stress=Stress(pentaelem,grids,materials,inputs)
0002
0003
0004
0005
0006
0007
0008
0009
0010 stress=zeros(6,1);
0011 volume=0;
0012
0013
0014 NDOF1=1;
0015 numgrids=6;
0016
0017
0018 matice=materials(pentaelem.matid).material;
0019 matpar=materials(end).constants;
0020
0021 gravity=matpar.g;
0022 rho_ice=matpar.rho_ice;
0023 B=matice.B;
0024
0025
0026 xyz_list=getgriddata(pentaelem,grids);
0027
0028
0029 [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0030 [surface_param surface_is_present]=recover_input(inputs,'surface');
0031 [flow_law_param flow_law_is_present]=recover_input(inputs,'B');
0032 [elementonstokes_param elementonstokes_is_present]=recover_input(inputs,'elementonstokes');
0033
0034
0035
0036 if ~velocity_is_present &~elementonstokes_is_present
0037 error('Stress error message: missing inputs!');
0038 end
0039
0040
0041 vxvyvz_list=zeros(numgrids,3);
0042 pressure=zeros(numgrids,1);
0043
0044
0045 for i=1:numgrids,
0046 doflist=grids(pentaelem.g(i)).grid.doflist;
0047 for j=1:3,
0048 dof=doflist(j);
0049 vxvyvz_list(i,j)=velocity_param(dof);
0050 end
0051 pressure(i)=velocity_param(doflist(4));
0052 dof=doflist(1);
0053 if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
0054 if(surface_is_present) surface_list(i)=surface_param(dof);end;
0055 end
0056
0057
0058 if flow_law_is_present,
0059 B_param=GetParameterValue(pentaelem,B_list,gauss_coord);
0060 matice.B=B_param; clear B_param.
0061 end
0062
0063
0064 area_order=2;
0065 num_vert_gauss=2;
0066 [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);
0067
0068
0069 for igarea=1:num_area_gauss,
0070 for igvert=1:num_vert_gauss,
0071
0072
0073 gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
0074 gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
0075
0076
0077 strainrate_g=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord);
0078 viscosity=GetViscosity3d(matice,strainrate_g);
0079 deviatoricstress=viscosity*strainrate_g;
0080
0081
0082 if elementonstokes_param(pentaelem.id)==stokesenum();
0083 pressure_g=GetParameterValue(pentaelem,pressure,gauss_coord);
0084 else
0085 z_av=mean(xyz_list(:,3));
0086 if surface_is_present,
0087 s_av=mean(surface_list);
0088 else
0089 s_av=mean(pentaelem.s);
0090 end
0091 pressure_g=deviatoricstress(3)+rho_ice*gravity*(s_av-z_av);
0092 end
0093
0094
0095 stress_g=deviatoricstress;
0096 stress_g([1;2;3],1)=stress_g(1:3,1)-pressure_g;
0097
0098
0099 Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
0100
0101 stress=stress+stress_g*Jdet*gauss_weight;
0102 volume=volume+Jdet*gauss_weight;
0103
0104 end
0105 end
0106
0107
0108 stress=stress/volume;