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