StressBedCompute

PURPOSE ^

STRESSBEDCOMPUTE - compute the stress on the basal boundary condition

SYNOPSIS ^

function stress_bed=StressBedCompute(m,inputs,type);

DESCRIPTION ^

STRESSBEDCOMPUTE - compute the stress on the basal boundary condition

   return a vector of size (numberofelements,1), holding the stress for 
   every element on bed

   Usage:
      stress_bed=StressBedCompute(m,inputs,type)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function stress_bed=StressBedCompute(m,inputs,type);
0002 %STRESSBEDCOMPUTE - compute the stress on the basal boundary condition
0003 %
0004 %   return a vector of size (numberofelements,1), holding the stress for
0005 %   every element on bed
0006 %
0007 %   Usage:
0008 %      stress_bed=StressBedCompute(m,inputs,type)
0009 
0010 %global variables
0011 global cluster gridset
0012 
0013 %recover fem model fields
0014 elements=m.elements;
0015 grids=m.grids;
0016 materials=m.materials;
0017 loads=m.loads;
0018 gridset=m.gridset;
0019 
0020 %figure out active elements that will take part in the stiffness and load generation
0021 [n1,n2]=GetNumberOfActiveElements(elements);
0022 
0023 %initialization
0024 stress_bed=zeros((n2-n1)+1,1);
0025 
0026 if strcmpi(type,'2d')
0027     disp('stress_bed not computed for 2d meshes')
0028     return
0029 end
0030 
0031 %initialize vectors
0032 stress_bed=struct('xx',[],'yy',[],'zz',[],'xy',[],'xz',[],'yz',[],'stress_n','stress_nn','normal_x',[],'normal_y',[],'normal_z',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'principalvalue3',[],'principalaxis3',[],'effectivevalue',[]);
0033 stress_bed1=zeros((n2-n1)+1,6);
0034 normal1=zeros((n2-n1)+1,3);
0035 stress_n1=zeros((n2-n1)+1,3);
0036 stress_nn1=zeros((n2-n1)+1,1);
0037 A1=zeros((n2-n1)+1,1); Vx1=zeros((n2-n1)+1,1); Vy1=zeros((n2-n1)+1,1); Vz1=zeros((n2-n1)+1,1);
0038 A2=zeros((n2-n1)+1,1); Vx2=zeros((n2-n1)+1,1); Vy2=zeros((n2-n1)+1,1); Vz2=zeros((n2-n1)+1,1);
0039 A3=zeros((n2-n1)+1,1); Vx3=zeros((n2-n1)+1,1); Vy3=zeros((n2-n1)+1,1); Vz3=zeros((n2-n1)+1,1);
0040 
0041 %Go through all elements and call the stress_bed routine, then compute eigen values and vector
0042 for n=n1:n2,
0043     if ~isempty(elements(n).element),
0044         [stress_vector,normal]=StressBed(elements(n).element,grids,materials,inputs);
0045 
0046         if ~isnan(stress_vector)
0047             stress_matrix=[stress_vector(1) stress_vector(4) stress_vector(5)
0048                       stress_vector(4)  stress_vector(2)  stress_vector(6)
0049                       stress_vector(5)  stress_vector(6)  stress_vector(3)];
0050             stress_n=stress_matrix*normal;
0051             stress_nn=normal'*stress_matrix*normal;
0052 
0053             %eigen values and vectors
0054             [directions,value]=eig(stress_matrix);
0055 
0056             %Plug into global vectors
0057             stress_bed1(n,:)=stress_vector;
0058             stress_n1(n,:)=stress_n;
0059             stress_nn1(n)=stress_nn;
0060             normal1(n,:)=normal';
0061             A1(n,1)=value(1,1); A2(n,1)=value(2,2); A3(n,1)=value(3,3);
0062             Vx1(n,1)=directions(1,1); Vx2(n,1)=directions(1,2); Vx3(n,1)=directions(1,3);
0063             Vy1(n,1)=directions(2,1); Vy2(n,1)=directions(2,2); Vy3(n,1)=directions(2,3);
0064             Vz1(n,1)=directions(3,1); Vz2(n,1)=directions(3,2); Vz3(n,1)=directions(3,3);
0065         else
0066             stress_bed1(n,:)=NaN*ones(1,6);
0067             stress_n1(n,:)=NaN*ones(1,3);
0068             stress_nn1(n)=NaN;
0069             normal1(n,:)=NaN*ones(1,3);
0070             A1(n,1)=NaN; A2(n,1)=NaN; A3(n,1)=NaN;
0071             Vx1(n,1)=NaN; Vx2(n,1)=NaN; Vx3(n,1)=NaN;
0072             Vy1(n,1)=NaN; Vy2(n,1)=NaN; Vy3(n,1)=NaN;
0073             Vz1(n,1)=NaN; Vz2(n,1)=NaN; Vz3(n,1)=NaN;
0074         end
0075     end
0076 end
0077 
0078 %plug results into outputs
0079 %NB: Matlab sorts the eigen value in increasing order, we want the reverse
0080 %components
0081 stress_bed.xx=stress_bed1(:,1);
0082 stress_bed.yy=stress_bed1(:,2);
0083 stress_bed.zz=stress_bed1(:,3);
0084 stress_bed.xy=stress_bed1(:,4);
0085 stress_bed.xz=stress_bed1(:,5);
0086 stress_bed.yz=stress_bed1(:,6);
0087 %principal axis
0088 stress_bed.principalvalue3=A1;
0089 stress_bed.principalaxis3=[Vx1 Vy1 Vz1];
0090 stress_bed.principalvalue2=A2;
0091 stress_bed.principalaxis2=[Vx2 Vy2 Vz2];
0092 stress_bed.principalvalue1=A3;
0093 stress_bed.principalaxis1=[Vx3 Vy3 Vz3];
0094 %norm or effective value
0095 stress_bed.effectivevalue=1/sqrt(2)*sqrt(stress_bed.xx.^2+stress_bed.yy.^2+stress_bed.zz.^2+2*stress_bed.xy.^2+2*stress_bed.xz.^2+2*stress_bed.yz.^2);
0096 %force
0097 stress_bed.normal_x=normal1(:,1);
0098 stress_bed.normal_y=normal1(:,2);
0099 stress_bed.normal_z=normal1(:,3);
0100 stress_bed.stress_n=stress_n1;
0101 stress_bed.stress_nn=stress_nn1;

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