Index: /issm/trunk-jpl/src/m/mech/analyticaldamage.py
===================================================================
--- /issm/trunk-jpl/src/m/mech/analyticaldamage.py	(revision 17621)
+++ /issm/trunk-jpl/src/m/mech/analyticaldamage.py	(revision 17622)
@@ -1,8 +1,8 @@
 import numpy as npy
-from pairoptions import pairoptions
 from averaging import averaging
 from plotmodel import plotmodel
+from thomasparams import thomasparams
 
-def analyticaldamage(md,*args):
+def analyticaldamage(md,**kwargs):
 	'''
 	ANALYTICALDAMAGE - compute damage for an ice shelf 
@@ -14,11 +14,13 @@
 	
 	   Available options:
-			- 'eq'			: analytical equation to use in the calculation.  Must be one of:
+			-eq			: analytical equation to use in the calculation.  Must be one of:
 									'Weertman1D' for a confined ice shelf free to flow in one direction
 									'Weertman2D' for an unconfined ice shelf free to spread in any direction
 									'Thomas' for a 2D ice shelf, taking into account full strain rate tensor (default)
-			- 'smoothing'	: the amount of smoothing to be applied to the strain rate data.
+			-smoothing	: the amount of smoothing to be applied to the strain rate data.
 									Type 'help averaging' for more information on its usage.
-			- 'sigmab'		: a compressive backstress term to be subtracted from the driving stress 
+			-coordsys	: coordinate system for calculating the strain rate
+						components. Must be one of:
+			-sigmab		: a compressive backstress term to be subtracted from the driving stress 
 									in the damage calculation
 	
@@ -34,76 +36,33 @@
 	
 	   Usage:
-	      [damage,B,backstress]=analyticaldamage(md,options)
+	      damage,B,backstress=analyticaldamage(md,kwargs)
 	
 	   Example:
-	      [damage,B,backstress]=analyticaldamage(md,'eq','Weertman2D','smoothing',2,'backstress',10e3)
+	      damage,B,backstress=analyticaldamage(md,eq='Weertman2D',smoothing=2,sigmab=10e3)
 	'''
+
+	#unpack kwargs
+	eq=kwargs.pop('eq','Thomas')
+	if 'eq' in kwargs: del kwargs['eq']
+	smoothing=kwargs.pop('smoothing',0)
+	if 'smoothing' in kwargs: del kwargs['smoothing']
+	coordsys=kwargs.pop('coordsys','longitudinal')
+	if 'coordsys' in kwargs: del kwargs['coordsys']
+	sigmab=kwargs.pop('sigmab',0)
+	if 'sigmab' in kwargs: del kwargs['sigmab']
+	assert len(kwargs)==0, 'error, unexpected or misspelled kwargs'
+
+	if len(sigmab)==1:
+		sigmab=sigmab*npy.ones((md.mesh.numberofvertices,))
 
 	# check inputs
 	if 'strainrate' not in md.results.__dict__:
 		raise StandardError('md.results.strainrate not present.  Calculate using md=mechanicalproperties(md,vx,vy)')
-	if md.mesh.dimension!=2:
-		raise StandardError('only 2D model supported currently')
+	if not '2d' in md.mesh.__doc__:
+		raise StandardError('only 2d (planview) 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'
 
-	# process options
-	options = pairoptions(*args)
-	eq = options.getfieldvalue('eq','Thomas')
-	smoothing = options.getfieldvalue('smoothing',0)
-	sigmab = options.getfieldvalue('sigmab',0)
-	if len(sigmab==1):
-		sigmab=sigmab*npy.ones(md.mesh.numberofvertices,)
-	
-	# average element strain rates onto vertices
-	e1=averaging(md,md.results.strainrate.principalvalue1,smoothing)/md.constants.yts # convert to s^-1
-	e2=averaging(md,md.results.strainrate.principalvalue2,smoothing)/md.constants.yts
-	exx=averaging(md,md.results.strainrate.xx,smoothing)/md.constants.yts
-	eyy=averaging(md,md.results.strainrate.yy,smoothing)/md.constants.yts
-	exy=averaging(md,md.results.strainrate.xy,smoothing)/md.constants.yts
-	
-	# checks: any of e1 or e2 equal to zero?
-	pos=npy.nonzero(e1==0)
-	if npy.any(pos==1):
-		print 'WARNING: first principal strain rate equal to zero.  Value set to 1e-13 s^-1'
-		e1[pos]=1e-13
-	pos=npy.nonzero(e2==0)
-	if npy.any(pos==1):
-		disp('WARNING: second principal strain rate equal to zero.  Value set to 1e-13 s^-1');
-		e2[pos]=1e-13
-	
-	## old method using principal strain rates {{{
-	## ex=maximum principal tensile strain rate
-	#ex=e1;
-	#a=e2./e1;
-	#pos=find(e1<0 & e2>0); # longitudinal compression and lateral tension
-	#a(pos)=e1(pos)./e2(pos);
-	#ex(pos)=e2(pos);
-	#pos2=find(e1<0 & e2<0 & abs(e1)<abs(e2)); # lateral and longitudinal compression
-	#a(pos2)=e1(pos2)./e2(pos2);
-	#ex(pos2)=e2(pos2);
-	#pos3=find(e1>0 & e2>0 & abs(e1)<abs(e2)); # lateral and longitudinal tension 
-	#a(pos3)=e1(pos3)./e2(pos3);
-	#ex(pos3)=e2(pos3);
-	#id=find(e1<0 & e2<0);
-	#a(id)=-a(id); # where both strain rates are compressive, enforce negative alpha
-	#
-	## }}}
-	
-	# new method using longitudinal strain rates defined by observed velocity vector
-	velangle=npy.arctan(md.initialization.vy/md.initialization.vx)
-	ex=0.5*(exx+eyy)+0.5*(exx-eyy)*npy.cos(2.*velangle)+exy*npy.sin(2.*velangle)
-	ey=exx+eyy-ex # trace of strain rate tensor is invariant
-	exy=-0.5*(exx-eyy)*npy.sin(2.*velangle)+exy*npy.cos(2.*velangle)
-	a=ey/ex
-	b=exy/ex
-	pos=npy.nonzero(npy.logical_and(ex<0,ey<0))
-	#length(pos)
-	a[pos]=-a[pos]
-	
-	# a < -1 in areas of strong lateral compression or longitudinal compression
-	# and theta is undefined at a = -2
-	pos=npy.nonzero(abs((abs(a)-2))<1e-3)
-	a[pos]=-2+1e-3
+	a,b,theta,ex=thomasparams(md,eq=eq,smoothing=smoothing,coordsys=coordsys)
 	
 	# spreading stress
@@ -117,18 +76,5 @@
 	n=averaging(md,md.materials.rheology_n,0)
 	
-	if eq=='Weertman1D':
-		theta=1./8*npy.ones(md.mesh.numberofvertices,)
-		a=npy.zeros(md.mesh.numberofvertices,)
-	elif eq=='Weertman2D':
-		theta=1./9*npy.ones(md.mesh.numberofvertices,1)
-		a=npy.ones(md.mesh.numberofvertices,)
-	elif eq=='Thomas':
-		theta=((1+a+a**2+b**2)**((n-1)/2))/(abs(2+a)**n)
-	else:
-		raise StandardError('argument passed to "eq" not valid.  Type "help analyticaldamage" for usage')
-	
-	D=1-(1+a+a**2+b**2)**((n-1)/(2*n))/abs(ex)**(1./n)*(T-sigmab)/B/(2+a)/npy.sign(ex)
-	
-	backstress=npy.zeros(md.mesh.numberofvertices,)
+	D=1.-(1.+a+a**2+b**2)**((n-1.)/(2.*n))/npy.abs(ex)**(1./n)*(T-sigmab)/B/(2.+a)/npy.sign(ex)
 	
 	# D>1 where (2+a).*sign(ex)<0, compressive regions where high backstress needed
@@ -136,18 +82,23 @@
 	D[pos]=0
 	
-	# backstress to bring damage to zero
-	backstress[pos]=T[pos]-(1-D[pos])*B[pos]*npy.sign(ex[pos])*(2+a[pos])*abs(ex[pos])**(1./n[pos])/(1+a[pos]+a[pos]**2)**((n[pos]-1)/2./n[pos])
+	backstress=npy.zeros((md.mesh.numberofvertices,))
+	
+	# backstress to bring D down to one 
+	backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*npy.sign(ex[pos])*(2.+a[pos])*npy.abs(ex[pos])**(1./n[pos])/(1.+a[pos]+a[pos]**2)**((n[pos]-1.)/2./n[pos])
 	
 	pos=npy.nonzero(D<0)
+	#mask=ismember(1:md.mesh.numberofvertices,pos);
 	D[pos]=0
-
+	
 	# backstress to bring negative damage to zero
-	backstress[pos]=T[pos]-(1-D[pos])*B[pos]*npy.sign(ex[pos])*(2+a[pos])*abs(ex[pos])**(1./n[pos])/(1+a[pos]+a[pos]**2)**((n[pos]-1)/2./n[pos])
+	backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*npy.sign(ex[pos])*(2.+a[pos])*npy.abs(ex[pos])**(1./n[pos])/(1.+a[pos]+a[pos]**2)**((n[pos]-1.)/2./n[pos])
 	
 	pos=npy.nonzero(backstress<0)
 	backstress[pos]=0
 	
-	# increased rigidity to bring negative damage to zero
-	B[pos]=npy.sign(ex[pos])/(2+a[pos])*(1+a[pos]+a[pos]**2)**((n[pos]-1)/2./n[pos])*T[pos]/abs(ex[pos])**(1./n[pos]);
+	# rigidity from Thomas relation for D=0 and backstress=0
+	B=npy.sign(ex)/(2.+a)*(1.+a+a**2)**((n-1.)/2./n)*T/(npy.abs(ex)**(1./n))
+	pos=npy.nonzero(B<0)
+	B[pos]=md.materials.rheology_B[pos]
 	
 	damage=D
Index: /issm/trunk-jpl/src/m/mech/backstressfrominversion.py
===================================================================
--- /issm/trunk-jpl/src/m/mech/backstressfrominversion.py	(revision 17622)
+++ /issm/trunk-jpl/src/m/mech/backstressfrominversion.py	(revision 17622)
@@ -0,0 +1,75 @@
+import numpy as npy
+from averaging import averaging
+from thomasparams import thomasparams
+
+def backstressfrominversion(md,**kwargs):
+	'''
+	Compute ice shelf backstress from inversion results.
+
+	This routine computes backstress based on the analytical formalism of
+	Thomas (1973) and Borstad et al. (2013, The Cryoshpere).  The model
+	must contain inversion results for ice rigidity.  Strain rates must
+	also be included, either from observed or modeled velocities.  Ice
+	rigidity B is assumed to be parameterized by the ice temperature in
+	md.materials.rheology_B.
+
+   Available options:
+		- 'tempmask'	: mask the inverted rigidity to be no more than
+							appropriate for the temperature of the ice?  
+							Boolean, defaults to false.
+		- 'smoothing'	: the amount of smoothing to be applied to the strain rate data.
+								Type 'help averaging' for more information on its
+								usage. Defaults to 0.
+		- 'coordsys'	: coordinate system for calculating the strain rate
+							components. Must be one of: 
+				'longitudinal': x axis aligned along a flowline at every point (default)
+				'principal': x axis aligned along maximum principal strain rate
+					at every point
+				'xy': x and y axes same as in polar stereographic projection 
+
+   Return values:
+		'backstress' is the inferred backstress necessary to balance the
+		analytical solution (keeping damage within its appropriate limits, e.g. D
+		in [0,1]).
+
+   Usage:
+      backstress=backstressfrominversion(md,options)
+
+   Example:
+      backstress=backstressfrominversion(md,'smoothing',2,'coordsys','longitudinal','tempmask',true);
+	'''
+
+	# unpack kwargs
+	tempmask=kwargs.pop('tempmask',False)
+	if 'tempmask' in kwargs: del kwargs['maxiter']
+	smoothing=kwargs.pop('smoothing',0)
+	if 'smoothing' in kwargs: del kwargs['smoothing']
+	coordsys=kwargs.pop('coordsys','longitudinal')
+	if 'coordsys' in kwargs: del kwargs['coordsys']
+	assert len(kwargs)==0, 'error, unexpected or misspelled kwargs'
+
+	# some checks
+	if not hasattr(md.results,'strainrate'):
+		raise StandardError('md.results.strainrate not present.  Calculate using md=mechanicalproperties(md,vx,vy)')
+	if not '2d' in md.mesh.__doc__:
+		raise StandardError('only 2d (planview) model supported currently')
+	if any(md.flowequation.element_equation!=2):
+		raise StandardError('Warning: the model has some non-SSA elements.  These will be treated like SSA elements')
+
+	T=0.5*md.materials.rho_ice*md.constants.g*(1-md.materials.rho_ice/md.materials.rho_water)*md.geometry.thickness
+	n=averaging(md,md.materials.rheology_n,0)
+	B=md.materials.rheology_B
+	Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar.reshape(-1,)
+	
+	a0,b0,theta0,ex0=thomasparams(md,eq='Thomas',smoothing=smoothing,coordsys=coordsys)
+	
+	if tempmask:
+		Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar
+		pos=npy.nonzero(Bi>md.materials.rheology_B)
+		Bi[pos]=md.materials.rheology_B[pos]
+	
+	# analytical backstress solution
+	backstress=T-Bi*npy.sign(ex0)*(2+a0)*npy.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))
+	backstress[npy.nonzero(backstress<0)]=0
+
+	return backstress
Index: /issm/trunk-jpl/src/m/mech/damagefrominversion.py
===================================================================
--- /issm/trunk-jpl/src/m/mech/damagefrominversion.py	(revision 17622)
+++ /issm/trunk-jpl/src/m/mech/damagefrominversion.py	(revision 17622)
@@ -0,0 +1,32 @@
+import numpy as npy
+
+def damagefrominversion(md):
+	'''
+	compute ice shelf damage from inversion results
+
+	This routine computes damage based on the analytical formalism of Borstad et
+	al. (2013, The Cryosphere).  The model must contain inversion results for
+	ice rigidity.  Ice rigidity B is assumed to be parameterized by the ice
+	temperature in md.materials.rheology_B. 
+	
+	Usage:
+		damage=damagefrominversion(md)
+	
+	Example:
+		damage=damagefrominversion(md)
+	'''
+
+	# check inputs
+	if not hasattr(md.results,'strainrate'):
+		raise StandardError('md.results.strainrate is not present.  Calculate using md=mechanicalproperties(md,vx,vy)')
+	if not '2d' in md.mesh.__doc__:
+		raise StandardError('only 2d (planview) model supported currently')
+	if any(md.flowequation.element_equation!=2):
+		raise StandardError('Warning: the model has some non-SSA elements.  These will be treated like SSA elements')
+
+	damage=npy.zeros((md.mesh.numberofvertices,))
+	# Damage where Bi softer than B(T)
+	pos=npy.nonzero(md.results.StressbalanceSolution.MaterialsRheologyBbar<md.materials.rheology_B)
+	damage[pos]=1.-md.results.StressbalanceSolution.MaterialsRheologyBbar[pos]/md.materials.rheology_B[pos]
+
+	return damage
Index: /issm/trunk-jpl/src/m/mech/mechanicalproperties.py
===================================================================
--- /issm/trunk-jpl/src/m/mech/mechanicalproperties.py	(revision 17621)
+++ /issm/trunk-jpl/src/m/mech/mechanicalproperties.py	(revision 17622)
@@ -23,7 +23,7 @@
 	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 md.mesh.dimension!=2:
+	#	raise StandardError('only 2D model supported currently')
 
 	if npy.any(md.flowequation.element_equation!=2):
Index: /issm/trunk-jpl/src/m/mech/thomasparams.py
===================================================================
--- /issm/trunk-jpl/src/m/mech/thomasparams.py	(revision 17622)
+++ /issm/trunk-jpl/src/m/mech/thomasparams.py	(revision 17622)
@@ -0,0 +1,145 @@
+import numpy as npy
+from averaging import averaging
+
+def thomasparams(md,**kwargs):
+	'''
+	compute Thomas' geometric parameters for an ice shelf 
+
+	This routine computes geometric parameters representing ratios between
+	components of the horizontal strain rate tensor for an ice shelf, as
+	originally developed in Thomas (1973).  The model must contain computed
+	strain rates, either from observed or modeled ice velocities.
+
+   Available options:
+	 -eq			: analytical equation to use in the calculation.  Must be one of:
+				'Thomas' for a 2D ice shelf, taking into account full strain rate
+					tensor (default)
+				'Weertman1D' for a confined ice shelf free to flow in one direction
+				'Weertman2D' for an unconfined ice shelf free to spread in any direction
+
+	 -smoothing	: an integer smoothing parameter for the averaging function
+						(default 0) Type 'help averaging' for more information on its usage.
+
+	 -coordsys	: coordinate system for calculating the strain rate
+						components. Must be one of:
+				'longitudinal': x axis aligned along a flowline at every point (default)
+				'principal': x axis aligned along maximum principal strain rate
+					at every point
+				'xy': x and y axes same as in polar stereographic projection 
+
+   Return values: 
+
+		'alpha' which is the ratio e_yy/e_xx between components of the strain
+		rate tensor
+
+		'beta' which is the ratio e_xy/e_xx between components of the strain rate
+		tensor
+
+		'theta' which is a combination of alpha and beta arising from the form of
+		the equivalent stress
+
+		'exx' is the strain rate along a coordinate system defined by 'coordsys' 
+
+		'sigxx' is the deviatoric stress along a coordinate system defined by 'coordsys' 
+
+   Usage: 
+		alpha,beta,theta,exx,sigxx=thomasparams(md)
+
+   Example: 
+		alpha,beta,theta,exx,sigxx=thomasparams(md,eq='Thomas',smoothing=2,coordsys='longitudinal')
+	'''
+
+	#unpack kwargs
+	eq=kwargs.pop('eq','Thomas')
+	if 'eq' in kwargs: del kwargs['eq']
+	smoothing=kwargs.pop('smoothing',0)
+	if 'smoothing' in kwargs: del kwargs['smoothing']
+	coordsys=kwargs.pop('coordsys','longitudinal')
+	if 'coordsys' in kwargs: del kwargs['coordsys']
+	assert len(kwargs)==0, 'error, unexpected or misspelled kwargs'
+
+	# some checks
+	if not hasattr(md.results,'strainrate'):
+		raise StandardError('md.results.strainrate not present.  Calculate using md=mechanicalproperties(md,vx,vy)')
+	if not '2d' in md.mesh.__doc__:
+		raise StandardError('only 2d (planview) model supported currently')
+	if any(md.flowequation.element_equation!=2):
+		raise StandardError('Warning: the model has some non-SSA elements.  These will be treated like SSA elements')
+
+	# average element strain rates onto vertices
+	e1=averaging(md,md.results.strainrate.principalvalue1,smoothing)/md.constants.yts # convert to s^-1
+	e2=averaging(md,md.results.strainrate.principalvalue2,smoothing)/md.constants.yts
+	exx=averaging(md,md.results.strainrate.xx,smoothing)/md.constants.yts
+	eyy=averaging(md,md.results.strainrate.yy,smoothing)/md.constants.yts
+	exy=averaging(md,md.results.strainrate.xy,smoothing)/md.constants.yts
+	
+	# checks: any of e1 or e2 equal to zero?
+	pos=npy.nonzero(e1==0)
+	if npy.any(pos==1):
+		print 'WARNING: first principal strain rate equal to zero.  Value set to 1e-13 s^-1'
+		e1[pos]=1.e-13
+	pos=npy.nonzero(e2==0)
+	if npy.any(pos==1):
+		print 'WARNING: second principal strain rate equal to zero.  Value set to 1e-13 s^-1'
+		e2[pos]=1.e-13
+	
+	# rheology
+	n=averaging(md,md.materials.rheology_n,0)
+	B=md.materials.rheology_B
+	
+	if coordsys=='principal':
+		b=npy.zeros((md.mesh.numberofvertices,))
+		ex=e1
+		a=e2/e1
+		pos=npy.nonzero(npy.logical_and(e1<0,e2>0)) # longitudinal compression and lateral tension
+		a[pos]=e1[pos]/e2[pos]
+		ex[pos]=e2[pos]
+		pos2=npy.nonzero(e1<0 & e2<0 & npy.abs(e1)<npy.abs(e2)) # lateral and longitudinal compression
+		a[pos2]=e1[pos2]/e2[pos2]
+		ex[pos2]=e2[pos2]
+		pos3=npy.nonzero(e1>0 & e2>0 & npy.abs(e1)<npy.abs(e2)) # lateral and longitudinal tension
+		a[pos3]=e1[pos3]/e2[pos3]
+		ex[pos3]=e2[pos3]
+		ind=npy.nonzero(e1<0 & e2<0)
+		a[ind]=-a[ind] # where both strain rates are compressive, enforce negative alpha
+		sigxx=(npy.abs(ex)/((1.+a+a**2)**((n-1.)/2.)))**(1./n)*B
+	elif coordsys=='xy':
+		ex=exx
+		a=eyy/exx
+		b=exy/exx
+	elif coordsys=='longitudinal':
+		# using longitudinal strain rates defined by observed velocity vector
+		velangle=npy.arctan(md.initialization.vy/md.initialization.vx)
+		pos=npy.nonzero(md.initialization.vx==0)
+		velangle[pos]=npy.pi/2
+		ex=0.5*(exx+eyy)+0.5*(exx-eyy)*npy.cos(2.*velangle)+exy*npy.sin(2.*velangle)
+		ey=exx+eyy-ex # trace of strain rate tensor is invariant
+		exy=-0.5*(exx-eyy)*npy.sin(2.*velangle)+exy*npy.cos(2.*velangle)
+		a=ey/ex
+		b=exy/ex
+		sigxx=abs(ex)**(1./n-1.)*ex/((1.+a+a**2+b**2)**((n-1.)/(2.*n)))*B
+	else:
+		raise ValueError('argument passed to "coordsys" not valid')
+	
+	# a < -1 in areas of strong lateral compression or longitudinal compression and
+	# theta flips sign at a = -2
+	pos=npy.nonzero(npy.abs((npy.abs(a)-2.))<1.e-3)
+	if len(pos)>0:
+		print 'Warning: ', len(pos), ' vertices have alpha within 1e-3 of -2'
+	a[pos]=-2+1e-3
+
+	if eq=='Weertman1D':
+		theta=1./8
+		a=npy.zeros((md.mesh.numberofvertices,))
+	elif eq=='Weertman2D':
+		theta=1./9
+		a=npy.ones((md.mesh.numberofvertices,1))
+	elif eq=='Thomas':
+		theta=((1.+a+a**2+b**2)**((n-1.)/2.))/(npy.abs(2.+a)**n)
+	else:
+		raise ValueError('argument passed to "eq" not valid')
+
+	alpha=a
+	beta=b
+
+	return alpha,beta,theta,ex