Changeset 13476

Timestamp:

09/28/12 13:33:15 (12 years ago)

Author:

cborstad

Message:

CHG: changed mechanicalproperties to sort eigenvalues in descending order by absolute value. Previously the first principal eigenalue was not guaranteed to have the greatest absolute value (as it should), a result of the way the eig() function in matlab sorts results.

File:

• issm/trunk-jpl/src/m/mech/mechanicalproperties.m (modified) (4 diffs)

issm/trunk-jpl/src/m/mech/mechanicalproperties.m

@@ -76 +76 @@
         %eigen values and vectors
         [directions,value]=eig(stress);
+        %sort by absolute value of eigenvalue in descending order
+        [val,idx]=sort(abs(diag(value)),'descend');
+        %re-order eigenvalues and associated vectors 
+        value=value(idx,idx);
+        directions=directions(:,idx);
         valuesstress(i,:)=[value(1,1) value(2,2)];
         directionsstress(i,:)=directions(:)';
         [directions,value]=eig(strain);
+        %same for strainrate
+        [val,idx]=sort(abs(diag(value)),'descend');
+        value=value(idx,idx);
+        directions=directions(:,idx);
         valuesstrain(i,:)=[value(1,1) value(2,2)];
         directionsstrain(i,:)=directions(:)';
@@ -89 +98 @@
 stress.yy=tau_yy;
 stress.xy=tau_xy;
-stress.principalvalue2=valuesstress(:,1);
-stress.principalaxis2=directionsstress(:,1:2);
-stress.principalvalue1=valuesstress(:,2);
-stress.principalaxis1=directionsstress(:,3:4);
+stress.principalvalue1=valuesstress(:,1);
+stress.principalaxis1=directionsstress(:,1:2);
+stress.principalvalue2=valuesstress(:,2);
+stress.principalaxis2=directionsstress(:,3:4);
 stress.effectivevalue=1/sqrt(2)*sqrt(stress.xx.^2+stress.yy.^2+2*stress.xy.^2);
 md.results.stress=stress;
@@ -100 +109 @@
 strainrate.yy=vy;
 strainrate.xy=uyvx;
-strainrate.principalvalue2=valuesstrain(:,1)*(365.25*24*3600); %strain rate in 1/a instead of 1/s
-strainrate.principalaxis2=directionsstrain(:,1:2);
-strainrate.principalvalue1=valuesstrain(:,2)*(365.25*24*3600); %strain rate in 1/a instead of 1/s
-strainrate.principalaxis1=directionsstrain(:,3:4);
+strainrate.principalvalue1=valuesstrain(:,1)*(365.25*24*3600); %strain rate in 1/a instead of 1/s
+strainrate.principalaxis1=directionsstrain(:,1:2);
+strainrate.principalvalue2=valuesstrain(:,2)*(365.25*24*3600); %strain rate in 1/a instead of 1/s
+strainrate.principalaxis2=directionsstrain(:,3:4);
 strainrate.effectivevalue=1/sqrt(2)*sqrt(strainrate.xx.^2+strainrate.yy.^2+2*strainrate.xy.^2);
 md.results.strainrate=strainrate;
@@ -111 +120 @@
 deviatoricstress.yy=tau_yy;
 deviatoricstress.xy=tau_xy;
-deviatoricstress.principalvalue2=valuesstress(:,1);
-deviatoricstress.principalaxis2=directionsstress(:,1:2);
-deviatoricstress.principalvalue1=valuesstress(:,2);
-deviatoricstress.principalaxis1=directionsstress(:,3:4);
+deviatoricstress.principalvalue1=valuesstress(:,1);
+deviatoricstress.principalaxis1=directionsstress(:,1:2);
+deviatoricstress.principalvalue2=valuesstress(:,2);
+deviatoricstress.principalaxis2=directionsstress(:,3:4);
 deviatoricstress.effectivevalue=1/sqrt(2)*sqrt(stress.xx.^2+stress.yy.^2+2*stress.xy.^2);
 md.results.deviatoricstress=deviatoricstress;