0001 function deviatoricstress=DeviatoricStress(pentaelem,grids,materials,inputs)
0002
0003
0004
0005
0006
0007
0008 deviatoricstress=zeros(6,1);
0009 volume=0;
0010
0011
0012 NDOF1=1;
0013 numgrids=6;
0014
0015
0016 matice=materials(pentaelem.matid).material;
0017 B=matice.B;
0018
0019
0020 xyz_list=getgriddata(pentaelem,grids);
0021
0022
0023 [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0024 [flow_law_param flow_law_is_present]=recover_input(inputs,'B');
0025
0026
0027
0028 if ~velocity_is_present,
0029 error('DeviatoricStress error message: input velocity not present!');
0030 end
0031
0032
0033 vxvyvz_list=zeros(numgrids,3);
0034
0035
0036 for i=1:numgrids,
0037 doflist=grids(pentaelem.g(i)).grid.doflist;
0038 for j=1:3,
0039 dof=doflist(j);
0040 vxvyvz_list(i,j)=velocity_param(dof);
0041 end
0042 dof=doflist(1);
0043 if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
0044 end
0045
0046
0047 if flow_law_is_present,
0048 B_param=GetParameterValue(pentaelem,B_list,gauss_coord);
0049 matice.B=B_param; clear B_param.
0050 end
0051
0052
0053 area_order=2;
0054 num_vert_gauss=2;
0055 [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);
0056
0057
0058 for igarea=1:num_area_gauss,
0059 for igvert=1:num_vert_gauss,
0060
0061
0062 gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
0063 gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
0064
0065
0066 strainrate_g=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord);
0067 viscosity=GetViscosity3d(matice,strainrate_g);
0068 deviatoricstress_g=viscosity*strainrate_g;
0069
0070
0071 Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
0072
0073 deviatoricstress=deviatoricstress+deviatoricstress_g*Jdet*gauss_weight;
0074 volume=volume+Jdet*gauss_weight;
0075
0076 end
0077 end
0078
0079
0080 deviatoricstress=deviatoricstress/volume;