Stress

PURPOSE ^

STRESS - compute the stress on each element

SYNOPSIS ^

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

DESCRIPTION ^

STRESS - compute the stress on each element

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

   See also STRESSBED, STRESSSURFACE, GETSTRAINRATE, GETSTRAINRATESTOKES, GETVISCOSITY3D

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function stress=Stress(pentaelem,grids,materials,inputs)
0002 %STRESS - compute the stress on each element
0003 %
0004 %   Usage:
0005 %      stress=Stress(pentaelem,grids,materials,inputs)
0006 %
0007 %   See also STRESSBED, STRESSSURFACE, GETSTRAINRATE, GETSTRAINRATESTOKES, GETVISCOSITY3D
0008 
0009     %initialize
0010     stress=zeros(6,1);
0011     volume=0;
0012 
0013     %some variables
0014     NDOF1=1;
0015     numgrids=6;
0016 
0017     %recover material parameters
0018     matice=materials(pentaelem.matid).material;
0019     matpar=materials(end).constants;
0020 
0021     gravity=matpar.g;
0022     rho_ice=matpar.rho_ice;
0023     B=matice.B;
0024 
0025     %Get all element grid data:
0026     xyz_list=getgriddata(pentaelem,grids);
0027 
0028     %recover extra inputs
0029     [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0030     [surface_param surface_is_present]=recover_input(inputs,'surface');
0031     [flow_law_param flow_law_is_present]=recover_input(inputs,'B');
0032     [elementonstokes_param elementonstokes_is_present]=recover_input(inputs,'elementonstokes');
0033 
0034     %we need velocities to compute thermal profiles (even if it is a zero
0035     %vector).
0036     if ~velocity_is_present &~elementonstokes_is_present
0037         error('Stress error message: missing inputs!');
0038     end
0039 
0040     %initialize vxvyvz_list
0041     vxvyvz_list=zeros(numgrids,3);
0042     pressure=zeros(numgrids,1);
0043 
0044     %Build row indices for elementary vector.
0045     for i=1:numgrids,
0046         doflist=grids(pentaelem.g(i)).grid.doflist;
0047         for j=1:3,
0048             dof=doflist(j);
0049             vxvyvz_list(i,j)=velocity_param(dof);
0050         end
0051         pressure(i)=velocity_param(doflist(4));    
0052         dof=doflist(1);
0053         if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
0054         if(surface_is_present) surface_list(i)=surface_param(dof);end;
0055     end
0056     
0057     %Update material parameter that deals with ice rigidity:
0058     if flow_law_is_present,
0059         B_param=GetParameterValue(pentaelem,B_list,gauss_coord);
0060         matice.B=B_param; clear B_param.
0061     end
0062     
0063     % Get gaussian points and weights
0064     area_order=2;
0065     num_vert_gauss=2;
0066     [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);
0067 
0068     %Start  looping on the number of gaussian points:
0069     for igarea=1:num_area_gauss,
0070         for igvert=1:num_vert_gauss,
0071 
0072             %Pick up the gaussian point and its weight:
0073             gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
0074             gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
0075             
0076             %Build deviatoric stress
0077             strainrate_g=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord);
0078             viscosity=GetViscosity3d(matice,strainrate_g);
0079             deviatoricstress=viscosity*strainrate_g;
0080 
0081             %Build Pressure
0082             if elementonstokes_param(pentaelem.id)==stokesenum();
0083                 pressure_g=GetParameterValue(pentaelem,pressure,gauss_coord);
0084             else %use pattyn's and MacAyeal's assumption: P=sigma_zz'+rho_ice*g*(s-z)
0085                 z_av=mean(xyz_list(:,3));
0086                 if surface_is_present,
0087                     s_av=mean(surface_list);
0088                 else
0089                     s_av=mean(pentaelem.s);
0090                 end
0091                 pressure_g=deviatoricstress(3)+rho_ice*gravity*(s_av-z_av);
0092             end
0093 
0094             %Build stress
0095             stress_g=deviatoricstress;
0096             stress_g([1;2;3],1)=stress_g(1:3,1)-pressure_g;         %stress=deviatoricstress-pressure
0097 
0098             %Get Jacobian determinant:
0099             Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
0100 
0101             stress=stress+stress_g*Jdet*gauss_weight;
0102             volume=volume+Jdet*gauss_weight;
0103 
0104         end %for ig=1:vert_gauss,
0105     end %for ig=1:num_area_gauss,
0106 
0107 %Divide stress, integrated over volume, by volume.
0108 stress=stress/volume;

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