Index: /issm/trunk/src/m/classes/public/mechanicalproperties.m
===================================================================
--- /issm/trunk/src/m/classes/public/mechanicalproperties.m	(revision 629)
+++ /issm/trunk/src/m/classes/public/mechanicalproperties.m	(revision 629)
@@ -0,0 +1,116 @@
+function md=mechanicalproperties(md,vx,vy)
+%MECHANICALPROPERTIES - compute stress and strain rate for a goven velocity
+%
+%   this routine computes the components of the stress tensor
+%   strain rate tensor and their respective principal directions.
+%   the results are in the model md: md.stress and md.strainrate
+%
+%   Usage:
+%      md=mechanicalproperties(md,vx,vy)
+%
+%   Example:
+%      md=mechanicalproperties(md,md.vx,md.vy);
+%      md=mechanicalproperties(md,md.vx_obs,md.vy_obs);
+
+%some checks
+if length(vx)~=md.numberofgrids | length(vy)~=md.numberofgrids,
+	error(['the input velocity should be of size ' num2str(md.numberofgrids) '!'])
+end
+if ~strcmpi(md.type,'2d')
+	error('only 2d model supported yet');
+end
+if any(md.elements_type(:,1)~=macayealenum),
+	disp('Warning: the model has some non macayeal elements. These will be treated like MacAyeal''s elements');
+end
+
+%initialization
+numberofelements=md.numberofelements;
+index=md.elements;
+summation=[1;1;1];
+alpha=zeros(numberofelements,3);
+beta=zeros(numberofelements,3);
+%gamma=zeros(numberofelements,3);
+directionsstress=zeros(numberofelements,4);
+directionsstrain=zeros(numberofelements,4);
+valuesstress=zeros(numberofelements,2);
+valuesstrain=zeros(numberofelements,2);
+
+%compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma
+x1=md.x(index(:,1)); x2=md.x(index(:,2)); x3=md.x(index(:,3));
+y1=md.y(index(:,1)); y2=md.y(index(:,2)); y3=md.y(index(:,3));
+invdet=1./(x1.*(y2-y3)-x2.*(y1-y3)+x3.*(y1-y2));
+alpha=[invdet.*(y2-y3) invdet.*(y3-y1) invdet.*(y1-y2)];
+beta =[invdet.*(x3-x2) invdet.*(x1-x3) invdet.*(x2-x1)];
+%gamma=[invdet.*(x2.*y3-x3.*y2) invdet.*(y1.*x3-y3.*x1) invdet.*(x1.*y2-x2.*y1)];
+clear invdet x1 x2 x3 y1 y2 y3
+
+%compute shear
+vxlist=vx(index)/md.yts;
+vylist=vy(index)/md.yts;
+ux=(vxlist.*alpha)*summation;
+uy=(vxlist.*beta)*summation;
+vx=(vylist.*alpha)*summation;
+vy=(vylist.*beta)*summation;						
+uyvx=(vx+uy)./2;
+clear vxlist vylist
+
+%compute viscosity
+nu=zeros(numberofelements,1);
+B_bar=md.B(index)*summation/3;
+power=(md.n-1)./(2*md.n);
+second_inv=(ux.^2+vy.^2+((uy+vx).^2)/4+ux.*vy);
+%some corrections
+location=find(second_inv~=0);
+nu(location)=B_bar(location)./(second_inv(location).^power(location));
+location=find(second_inv==0 & power~=0);
+nu(location)=10^18; 	%arbitrary maximum viscosity to apply where there is no effective shear
+location=find(second_inv==0 & power==0);
+nu(location)=B_bar(location);
+clear B_bar location second_inv power
+mean(nu)
+
+%compute stress
+tau_xx=nu.*ux;
+tau_yy=nu.*vy;
+tau_xy=nu.*uyvx;
+
+%compute principal properties of stress
+for i=1:numberofelements,
+
+	%compute stress and strainrate matrices
+	stress=[tau_xx(i) tau_xy(i)
+	tau_xy(i)  tau_yy(i)];
+	strain=[ux(i) uyvx(i)
+	uyvx(i)  vy(i)];
+
+	%eigen values and vectors
+	[directions,value]=eig(stress);
+	valuesstress(i,:)=[value(1,1) value(2,2)];
+	directionsstress(i,:)=directions(:)';
+	[directions,value]=eig(strain);
+	valuesstrain(i,:)=[value(1,1) value(2,2)];
+	directionsstrain(i,:)=directions(:)';
+end
+
+%plug onto the model
+%NB: Matlab sorts the eigen value in increasing order, we want the reverse
+stress=struct('xx',[],'yy',[],'xy',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'effectivevalue',[]);
+stress.xx=tau_xx;
+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.effectivevalue=1/sqrt(2)*sqrt(stress.xx.^2+stress.yy.^2+2*stress.xy.^2);
+strainrate=struct('xx',[],'yy',[],'xy',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'effectivevalue',[]);
+strainrate.xx=ux;
+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.effectivevalue=1/sqrt(2)*sqrt(strainrate.xx.^2+strainrate.yy.^2+2*strainrate.xy.^2);
+md.stress=stress;
+md.strainrate=strainrate;
Index: /issm/trunk/src/m/classes/public/process_solve_options.m
===================================================================
--- /issm/trunk/src/m/classes/public/process_solve_options.m	(revision 628)
+++ /issm/trunk/src/m/classes/public/process_solve_options.m	(revision 629)
@@ -51,18 +51,19 @@
 %check solution type is supported
 if ~(strcmpi(analysis_type,'control') |  ...
-	strcmpi(analysis_type,'diagnostic') |  ...
-	strcmpi(analysis_type,'prognostic') |  ...
-	strcmpi(analysis_type,'thermal') |  ...
-	strcmpi(analysis_type,'mesh') |  ...
-	strcmpi(analysis_type,'mesh2grid') |  ...
-	strcmpi(analysis_type,'transient') ),
+		strcmpi(analysis_type,'diagnostic') |  ...
+		strcmpi(analysis_type,'prognostic') |  ...
+		strcmpi(analysis_type,'thermal') |  ...
+		strcmpi(analysis_type,'parameters') |  ...
+		strcmpi(analysis_type,'mesh') |  ...
+		strcmpi(analysis_type,'mesh2grid') |  ...
+		strcmpi(analysis_type,'transient') ),
 	error(['process_solve_options error message: analysis_type ' analysis_type ' not supported yet!']);
 end
 if ~(strcmpi(sub_analysis_type,'none') |  ...
-	strcmpi(sub_analysis_type,'steady') |  ...
-	strcmpi(sub_analysis_type,'horiz') |  ...
-	strcmpi(sub_analysis_type,'vert') |  ...
-	strcmpi(sub_analysis_type,'') |  ...
-	strcmpi(sub_analysis_type,'transient') ),
+		strcmpi(sub_analysis_type,'steady') |  ...
+		strcmpi(sub_analysis_type,'horiz') |  ...
+		strcmpi(sub_analysis_type,'vert') |  ...
+		strcmpi(sub_analysis_type,'') |  ...
+		strcmpi(sub_analysis_type,'transient') ),
 	error(['process_solve_options error message: sub_analysis_type ' sub_analysis_type ' not supported yet!']);
 end
Index: /issm/trunk/src/m/classes/public/solve.m
===================================================================
--- /issm/trunk/src/m/classes/public/solve.m	(revision 628)
+++ /issm/trunk/src/m/classes/public/solve.m	(revision 629)
@@ -72,4 +72,8 @@
 elseif strcmpi(md.analysis_type,'thermal'),
 	md=thermal(md);
+
+elseif strcmpi(md.analysis_type,'parameters'),
+	md=parameters(md);
+
 else
 	error('solution type not supported for this model configuration.');