Stress

PURPOSE ^

STRESS - compute the stress on each element

SYNOPSIS ^

function stress=Stress(triaelem,grids,materials,inputs)

DESCRIPTION ^

STRESS - compute the stress on each element

   Usage:
      stress=Stress(triaelem,grids,materials,inputs)

   See also GETSTRAINRATE, GETVISCOSITY2D

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function stress=Stress(triaelem,grids,materials,inputs)
0002 %STRESS - compute the stress on each element
0003 %
0004 %   Usage:
0005 %      stress=Stress(triaelem,grids,materials,inputs)
0006 %
0007 %   See also GETSTRAINRATE, GETVISCOSITY2D
0008 
0009     %initialize
0010     stress=zeros(3,1);
0011     volume=0;
0012 
0013     %some variables
0014     NDOF1=1;
0015     numgrids=3;
0016 
0017     %recover material parameters
0018     matice=materials(triaelem.matid).material;
0019     B=matice.B;
0020 
0021     %Get all element grid data:
0022     xyz_list=getgriddata(triaelem,grids);
0023 
0024     %recover extra inputs
0025     [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0026     [flow_law_param flow_law_is_present]=recover_input(inputs,'B');
0027 
0028     %we need velocities to compute thermal profiles (even if it is a zero
0029     %vector).
0030     if ~velocity_is_present,
0031         error('Stress error message: input velocity not present!');
0032     end
0033 
0034     %initialize vxvyvz_list
0035     vxvy_list=zeros(numgrids,2);
0036     pressure=zeros(numgrids,1);
0037 
0038     %Build row indices for elementary vector.
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     %Update material parameter that deals with ice rigidity:
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     % Get gaussian points and weights
0057     [num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights]=GaussTria(2);
0058 
0059     %Start  looping on the number of gaussian points:
0060     for igarea=1:num_area_gauss,
0061 
0062         %Pick up the gaussian point and its weight:
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         %Build Stress
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;         %stress=deviatoricstress-pressure
0074 
0075         %Get Jacobian determinant:
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 %for ig=1:num_area_gauss,
0083 
0084 %Divide stress, integrated over volume, by volume.
0085 stress=stress/volume;

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003