0001 function deviatoricstress=DeviatoricStress(triaelem,grids,materials,inputs)
0002
0003
0004
0005
0006
0007
0008 deviatoricstress=zeros(3,1);
0009 volume=0;
0010
0011
0012 NDOF1=1;
0013 numgrids=3;
0014
0015
0016 matice=materials(triaelem.matid).material;
0017 B=matice.B;
0018
0019
0020 xyz_list=getgriddata(triaelem,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 vxvy_list=zeros(numgrids,2);
0034
0035
0036 for i=1:numgrids,
0037 doflist=grids(triaelem.g(i)).grid.doflist;
0038 for j=1:2,
0039 dof=doflist(j);
0040 vxvy_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(triaelem,B_list,gauss_coord);
0049 matice.B=B_param; clear B_param.
0050 end
0051
0052
0053 [num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights]=GaussTria(2);
0054
0055
0056 for igarea=1:num_area_gauss,
0057
0058
0059 gauss_weight=area_gauss_weights(igarea);
0060 gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea)];
0061
0062
0063 deviatoricstress_g=zeros(3,1);
0064 epsilon=GetStrainRate(triaelem,vxvy_list,xyz_list,gauss_coord);
0065 viscosity=GetViscosity2d(matice,epsilon);
0066 deviatoricstress_g=viscosity*epsilon;
0067
0068
0069 Jdet=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
0070
0071 deviatoricstress=deviatoricstress+deviatoricstress_g*Jdet*gauss_weight;
0072
0073 volume=volume+Jdet*gauss_weight;
0074
0075 end
0076
0077
0078 deviatoricstress=deviatoricstress/volume;