DeviatoricStress

PURPOSE ^

DEVIATORICSTRESS - computes the deviatoric stress of a pentaelem

SYNOPSIS ^

function deviatoricstress=DeviatoricStress(pentaelem,grids,materials,inputs)

DESCRIPTION ^

DEVIATORICSTRESS - computes the deviatoric stress of a pentaelem
 
   Usage:
      deviatoricstress=DeviatoricStress(pentaelem,grids,materials,inputs)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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