Index: /issm/trunk-jpl/src/m/mech/mechanicalproperties.py
===================================================================
--- /issm/trunk-jpl/src/m/mech/mechanicalproperties.py	(revision 15844)
+++ /issm/trunk-jpl/src/m/mech/mechanicalproperties.py	(revision 15844)
@@ -0,0 +1,141 @@
+import numpy as npy
+from GetNodalFunctionsCoeff import GetNodalFunctionsCoeff
+from results import results
+
+def 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.results
+	
+   Usage:
+      md=mechanicalproperties(md,vx,vy)
+	
+   Example:
+      md=mechanicalproperties(md,md.initialization.vx,md.initialization.vy)
+      md=mechanicalproperties(md,md.inversion.vx_obs,md.inversion.vy_obs)
+	"""
+
+	#some checks
+	if len(vx)!=md.mesh.numberofvertices or len(vy)!=md.mesh.numberofvertices:
+		raise ValueError('the input velocity should be of size ' + md.mesh.numberofvertices)
+
+	if md.mesh.dimension!=2:
+		raise StandardError('only 2D model supported currently')
+
+	if npy.any(md.flowequation.element_equation!=2):
+		print 'Warning: the model has some non SSA elements. These will be treated like SSA elements'
+	
+	#initialization
+	numberofelements=md.mesh.numberofelements
+	index=md.mesh.elements
+	summation=npy.array([[1],[1],[1]])
+	directionsstress=npy.zeros((numberofelements,4))
+	directionsstrain=npy.zeros((numberofelements,4))
+	valuesstress=npy.zeros((numberofelements,2))
+	valuesstrain=npy.zeros((numberofelements,2))
+	
+	#compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma
+	alpha,beta=GetNodalFunctionsCoeff(index,md.mesh.x,md.mesh.y)[0:2]
+	
+	#compute shear
+	vxlist=vx[index-1]/md.constants.yts
+	vylist=vy[index-1]/md.constants.yts
+	ux=npy.dot((vxlist*alpha),summation).reshape(-1,)
+	uy=npy.dot((vxlist*beta),summation).reshape(-1,)
+	vx=npy.dot((vylist*alpha),summation).reshape(-1,)
+	vy=npy.dot((vylist*beta),summation).reshape(-1,)
+	uyvx=(vx+uy)/2.
+	#clear vxlist vylist
+	
+	#compute viscosity
+	nu=npy.zeros((numberofelements,))
+	B_bar=npy.dot(md.materials.rheology_B[index-1],summation/3.).reshape(-1,)
+	power=(md.materials.rheology_n-1.)/(2.*md.materials.rheology_n)
+	second_inv=(ux**2.+vy**2.+((uy+vx)**2.)/4.+ux*vy).reshape(-1,)
+	
+	#some corrections
+	location=npy.nonzero(npy.logical_and(second_inv==0,power!=0))
+	nu[location]=10^18 	#arbitrary maximum viscosity to apply where there is no effective shear
+	
+	if 'matice' in md.materials.__module__:
+		location=npy.nonzero(second_inv)
+		nu[location]=B_bar[location]/(second_inv[location]**power[location])
+		location=npy.nonzero(npy.logical_and(second_inv==0,power==0))
+		nu[location]=B_bar[location]
+	elif 'matdamageice' in md.materials.__module__:
+		Zinv=npy.dot(md.materials.rheology_Z[index-1],summation/3.).reshape(-1,)
+		location=npy.nonzero(second_inv)
+		nu[location]=Zinv[location]*B_bar[location]/npy.power(second_inv[location],power[location])
+		location=npy.nonzero(npy.logical_and(second_inv==0,power==0))
+		nu[location]=Zinv[location]*B_bar[location]
+		#clear Zinv
+	else:
+		raise StandardError('class of md.materials (' + md.materials.__module__ + ') not recognized or not supported')
+	
+	#compute stress
+	tau_xx=nu*ux
+	tau_yy=nu*vy
+	tau_xy=nu*uyvx
+	
+	#compute principal properties of stress
+	for i in npy.arange(numberofelements):
+	
+		#compute stress and strainrate matrices
+		stress=npy.array([ [tau_xx[i], tau_xy[i]], [tau_xy[i], tau_yy[i]] ])
+		strain=npy.array([ [ux[i], uyvx[i]], [uyvx[i], vy[i]] ])
+	
+		#eigenvalues and vectors for stress
+		value,directions=npy.linalg.eig(stress);
+		idx=abs(value).argsort()[::-1] # sort in descending order
+		value=value[idx]
+		directions=directions[:,idx]
+		valuesstress[i,:]=[value[0],value[1]]
+		directionsstress[i,:]=directions.transpose().flatten()
+
+		#eigenvalues and vectors for strain
+		value,directions=npy.linalg.eig(strain);
+		idx=abs(value).argsort()[::-1] # sort in descending order
+		value=value[idx]
+		directions=directions[:,idx]
+		valuesstrain[i,:]=[value[0],value[1]]
+		directionsstrain[i,:]=directions.transpose().flatten()
+
+	##plug onto the model
+	##NB: Matlab sorts the eigen value in increasing order, we want the reverse
+	stress=results()
+	stress.xx=tau_xx
+	stress.yy=tau_yy
+	stress.xy=tau_xy
+	stress.principalvalue1=valuesstress[:,0]
+	stress.principalaxis1=directionsstress[:,0:2]
+	stress.principalvalue2=valuesstress[:,1]
+	stress.principalaxis2=directionsstress[:,2:4]
+	stress.effectivevalue=1./npy.sqrt(2.)*npy.sqrt(stress.xx**2+stress.yy**2+2.*stress.xy**2)
+	md.results.stress=stress
+	
+	strainrate=results()
+	strainrate.xx=ux*md.constants.yts #strain rate in 1/a instead of 1/s
+	strainrate.yy=vy*md.constants.yts 
+	strainrate.xy=uyvx*md.constants.yts 
+	strainrate.principalvalue1=valuesstrain[:,0]*md.constants.yts 
+	strainrate.principalaxis1=directionsstrain[:,0:2]
+	strainrate.principalvalue2=valuesstrain[:,1]*md.constants.yts 
+	strainrate.principalaxis2=directionsstrain[:,2:4]
+	strainrate.effectivevalue=1./npy.sqrt(2.)*npy.sqrt(strainrate.xx**2+strainrate.yy**2+2.*strainrate.xy**2)
+	md.results.strainrate=strainrate
+	
+	deviatoricstress=results()
+	deviatoricstress.xx=tau_xx
+	deviatoricstress.yy=tau_yy
+	deviatoricstress.xy=tau_xy
+	deviatoricstress.principalvalue1=valuesstress[:,0]
+	deviatoricstress.principalaxis1=directionsstress[:,1:2]
+	deviatoricstress.principalvalue2=valuesstress[:,1]
+	deviatoricstress.principalaxis2=directionsstress[:,2:4]
+	deviatoricstress.effectivevalue=1./npy.sqrt(2.)*npy.sqrt(stress.xx**2+stress.yy**2+2.*stress.xy**2)
+	md.results.deviatoricstress=deviatoricstress
+
+	return md