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