StressBed

PURPOSE ^

STRESSBED - compute the stress on the bed

SYNOPSIS ^

function [stress_bed,normal]=StressBed(pentaelem,grids,materials,inputs)

DESCRIPTION ^

STRESSBED - compute the stress on the bed

   This routine compute the stress_bed of each element, return NaN if not on bed

   Usage:
      [stress_bed,normal]=StressBed(pentaelem,grids,materials,inputs)

   See also STRESS, STRESSSURFACE, GETSTRAINRATE, GETVISCOSITY3D

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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