Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 17390)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 17391)
@@ -117,7 +117,163 @@
 ElementMatrix* DamageEvolutionAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
+	/* Check if ice in element */
+	if(!element->IsIceInElement()) return NULL;
+
+	/*Intermediaries*/
+	Element*    basalelement;
+	int         meshtype,dim;
+	int         stabilization;
+	IssmDouble  Jdet,dt,D_scalar,h;
+	IssmDouble  vel,vx,vy,dvxdx,dvydy,dvx[2],dvy[2];
+	IssmDouble *xyz_list  = NULL;
+
+	/*Get problem dimension and basal element*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			dim = 2;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnBed()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			dim = 2;
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = basalelement->GetNumberOfNodes();
+
+	/*Initialize Element vector*/
+	ElementMatrix* Ke     = basalelement->NewElementMatrix();
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
+
+	/*Retrieve all inputs and parameters*/
+	basalelement->GetVerticesCoordinates(&xyz_list);
+	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
+	basalelement->FindParam(&stabilization,DamageStabilizationEnum);
+	Input* vxaverage_input=NULL;
+	Input* vyaverage_input=NULL;
+	if(meshtype==Mesh2DhorizontalEnum){
+		vxaverage_input=basalelement->GetInput(VxEnum); _assert_(vxaverage_input);
+		vyaverage_input=basalelement->GetInput(VyEnum); _assert_(vyaverage_input);
+	}
+	else{
+		if(dim==1){
+			vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input);
+		}
+		if(dim==2){
+			vxaverage_input=element->GetInput(VxAverageEnum); _assert_(vxaverage_input);
+			vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input);
+		}
+	}
+	h=basalelement->CharacteristicLength();
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=basalelement->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		basalelement->NodalFunctions(basis,gauss);
+		
+		vxaverage_input->GetInputValue(&vx,gauss);
+		vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+		if(dim==2){
+			vyaverage_input->GetInputValue(&vy,gauss);
+			vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+		}
+
+		D_scalar=gauss->weight*Jdet;
+		TripleMultiply(basis,1,numnodes,1,
+					&D_scalar,1,1,0,
+					basis,1,numnodes,0,
+					&Ke->values[0],1);
+
+		GetB(B,basalelement,dim,xyz_list,gauss);
+		GetBprime(Bprime,basalelement,dim,xyz_list,gauss);
+
+		dvxdx=dvx[0];
+		if(dim==2) dvydy=dvy[1];
+		D_scalar=dt*gauss->weight*Jdet;
+
+		D[0*dim+0]=D_scalar*dvxdx;
+		if(dim==2) D[1*dim+1]=D_scalar*dvydy;
+
+		TripleMultiply(B,dim,numnodes,1,
+					D,dim,dim,0,
+					B,dim,numnodes,0,
+					&Ke->values[0],1);
+
+		D[0*dim+0]=D_scalar*vx;
+		if(dim==2) D[1*dim+1]=D_scalar*vy;
+
+		TripleMultiply(B,dim,numnodes,1,
+					D,dim,dim,0,
+					Bprime,dim,numnodes,0,
+					&Ke->values[0],1);
+
+		if(stabilization==2){
+			if(dim==1){
+				vel=fabs(vx)+1.e-8;
+				D[0]=h/(2*vel)*vx*vx;
+			}
+			else{
+				/*Streamline upwinding*/
+				vel=sqrt(vx*vx+vy*vy)+1.e-8;
+				D[0*dim+0]=h/(2*vel)*vx*vx;
+				D[1*dim+0]=h/(2*vel)*vy*vx;
+				D[0*dim+1]=h/(2*vel)*vx*vy;
+				D[1*dim+1]=h/(2*vel)*vy*vy;
+			}
+		}
+		else if(stabilization==1){
+			vxaverage_input->GetInputAverage(&vx);
+			if(dim==2) vyaverage_input->GetInputAverage(&vy);
+			D[0*dim+0]=h/2.0*fabs(vx);
+			if(dim==2) D[1*dim+1]=h/2.0*fabs(vy);
+		}
+		if(stabilization==1 || stabilization==2){
+			if(dim==1) D[0]=D_scalar*D[0];
+			else{
+				D[0*dim+0]=D_scalar*D[0*dim+0];
+				D[1*dim+0]=D_scalar*D[1*dim+0];
+				D[0*dim+1]=D_scalar*D[0*dim+1];
+				D[1*dim+1]=D_scalar*D[1*dim+1];
+			}
+
+			TripleMultiply(Bprime,dim,numnodes,1,
+						D,dim,dim,0,
+						Bprime,dim,numnodes,0,
+						&Ke->values[0],1);
+		}
+
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(D);
+	delete gauss;
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+	return Ke;
+}/*}}}*/
+ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	if(!element->IsIceInElement()) return NULL;
+
 	/*Intermediaries*/
 	int      meshtype;
 	Element* basalelement;
+	IssmDouble  Jdet,dt;
+	IssmDouble  f,damage;
+	IssmDouble* xyz_list = NULL;
 
 	/*Get basal element*/
@@ -134,35 +290,18 @@
 	}
 
-	/*Intermediaries */
-	int         stabilization;
-	IssmDouble  Jdet,dt,D_scalar;
-	IssmDouble  vel,vx,vy,dvxdx,dvydy,dvx[2],dvy[2];
-	IssmDouble *xyz_list  = NULL;
-
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = basalelement->GetNumberOfNodes();
 
-	/*Initialize Element vector*/
-	ElementMatrix* Ke     = basalelement->NewElementMatrix(NoneApproximationEnum);
-	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
-	IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
-	IssmDouble     D[2][2]={0.};
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = basalelement->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
 	basalelement->GetVerticesCoordinates(&xyz_list);
 	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
-	basalelement->FindParam(&stabilization,DamageStabilizationEnum);
-	Input* vxaverage_input=NULL;
-	Input* vyaverage_input=NULL;
-	if(meshtype==Mesh2DhorizontalEnum){
-		vxaverage_input=basalelement->GetInput(VxEnum); _assert_(vxaverage_input);
-		vyaverage_input=basalelement->GetInput(VyEnum); _assert_(vyaverage_input);
-	}
-	else{
-		vxaverage_input=basalelement->GetInput(VxAverageEnum); _assert_(vxaverage_input);
-		vyaverage_input=basalelement->GetInput(VyAverageEnum); _assert_(vyaverage_input);
-	}
-	IssmDouble h=basalelement->CharacteristicLength();
+	this->CreateDamageFInput(basalelement);
+	Input* damaged_input = basalelement->GetInput(DamageDEnum); _assert_(damaged_input);
+	Input* damagef_input = basalelement->GetInput(DamageFEnum); _assert_(damagef_input);
+
 
 	/* Start  looping on the number of gaussian points: */
@@ -173,126 +312,11 @@
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		basalelement->NodalFunctions(basis,gauss);
-		D_scalar=gauss->weight*Jdet;
-
-		vxaverage_input->GetInputValue(&vx,gauss);
-		vyaverage_input->GetInputValue(&vy,gauss);
-		vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-		vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
-
-		TripleMultiply(basis,1,numnodes,1,
-					&D_scalar,1,1,0,
-					basis,1,numnodes,0,
-					&Ke->values[0],1);
-
-		GetB(B,element,xyz_list,gauss);
-		GetBprime(Bprime,element,xyz_list,gauss);
-
-		dvxdx=dvx[0];
-		dvydy=dvy[1];
-		D_scalar=dt*gauss->weight*Jdet;
-
-		D[0][0]=D_scalar*dvxdx;
-		D[1][1]=D_scalar*dvydy;
-		TripleMultiply(B,2,numnodes,1,
-					&D[0][0],2,2,0,
-					B,2,numnodes,0,
-					&Ke->values[0],1);
-
-		D[0][0]=D_scalar*vx;
-		D[1][1]=D_scalar*vy;
-		TripleMultiply(B,2,numnodes,1,
-					&D[0][0],2,2,0,
-					Bprime,2,numnodes,0,
-					&Ke->values[0],1);
-
-		if(stabilization==2){
-			/*Streamline upwinding*/
-			vel=sqrt(vx*vx+vy*vy)+1.e-8;
-			D[0][0]=h/(2*vel)*vx*vx;
-			D[1][0]=h/(2*vel)*vy*vx;
-			D[0][1]=h/(2*vel)*vx*vy;
-			D[1][1]=h/(2*vel)*vy*vy;
-		}
-		else if(stabilization==1){
-			/*SSA*/
-			vxaverage_input->GetInputAverage(&vx);
-			vyaverage_input->GetInputAverage(&vy);
-			D[0][0]=h/2.0*fabs(vx);
-			D[0][1]=0.;
-			D[1][0]=0.;
-			D[1][1]=h/2.0*fabs(vy);
-		}
-		if(stabilization==1 || stabilization==2){
-			D[0][0]=D_scalar*D[0][0];
-			D[1][0]=D_scalar*D[1][0];
-			D[0][1]=D_scalar*D[0][1];
-			D[1][1]=D_scalar*D[1][1];
-			TripleMultiply(Bprime,2,numnodes,1,
-						&D[0][0],2,2,0,
-						Bprime,2,numnodes,0,
-						&Ke->values[0],1);
-		}
-
-	}
-
-	/*Clean up and return*/
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(basis);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(Bprime);
-	delete gauss;
-	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
-	return Ke;
-}/*}}}*/
-ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/
-
-	/*Intermediaries*/
-	int      meshtype;
-	Element* basalelement;
-
-	/*Get basal element*/
-	element->FindParam(&meshtype,MeshTypeEnum);
-	switch(meshtype){
-		case Mesh2DhorizontalEnum:
-			basalelement = element;
-			break;
-		case Mesh3DEnum:
-			if(!element->IsOnBed()) return NULL;
-			basalelement = element->SpawnBasalElement();
-			break;
-		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
-	}
-
-	/*Intermediaries */
-	IssmDouble  Jdet,dt;
-	IssmDouble  f,damage;
-	IssmDouble* xyz_list = NULL;
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Initialize Element vector and other vectors*/
-	ElementVector* pe    = element->NewElementVector();
-	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
-
-	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinates(&xyz_list);
-	element->FindParam(&dt,TimesteppingTimeStepEnum);
-	this->CreateDamageFInput(element);
-	Input* damaged_input = element->GetInput(DamageDbarEnum); _assert_(damaged_input);
-	Input* damagef_input = element->GetInput(DamageFEnum);    _assert_(damagef_input);
-
-	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-
-		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		element->NodalFunctions(basis,gauss);
 
 		damaged_input->GetInputValue(&damage,gauss);
 		damagef_input->GetInputValue(&f,gauss);
 
-		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i];
+		for(int i=0;i<numnodes;i++){
+			pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i];
+		}
 	}
 
@@ -304,5 +328,5 @@
 	return pe;
 }/*}}}*/
-void DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
 	 * For node i, Bi can be expressed in the actual coordinate system
@@ -324,6 +348,7 @@
 	/*Build B: */
 	for(int i=0;i<numnodes;i++){
-		B[numnodes*0+i] = basis[i];
-		B[numnodes*1+i] = basis[i];
+		for(int j=0;j<dim;j++){
+			B[numnodes*j+i] = basis[i];
+		}
 	}
 
@@ -331,5 +356,5 @@
 	xDelete<IssmDouble>(basis);
 }/*}}}*/
-void DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
 	 * For node i, Bi' can be expressed in the actual coordinate system
@@ -346,11 +371,12 @@
 
 	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
 	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 	/*Build B': */
 	for(int i=0;i<numnodes;i++){
-		Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
-		Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
+		for(int j=0;j<dim;j++){
+			Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
+		}
 	}
 
@@ -364,33 +390,44 @@
 void DamageEvolutionAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
+	int meshtype;
 	IssmDouble  max_damage;
 	int			*doflist = NULL;
-
+	Element*   basalelement=NULL;
+
+	element->FindParam(&meshtype,MeshTypeEnum);
+	if(meshtype!=Mesh2DhorizontalEnum){
+		if(!element->IsOnBed()) return;
+		basalelement=element->SpawnBasalElement();
+	}
+	else{
+		basalelement = element;
+	}
 	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
+	int numnodes = basalelement->GetNumberOfNodes();
 
 	/*Fetch dof list and allocate solution vector*/
-	element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-	IssmDouble* values    = xNew<IssmDouble>(numnodes);
+	basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	IssmDouble* newdamage = xNew<IssmDouble>(numnodes);
 
 	/*Get user-supplied max_damage: */
-	element->FindParam(&max_damage,DamageMaxDamageEnum);
+	basalelement->FindParam(&max_damage,DamageMaxDamageEnum);
 
 	/*Use the dof list to index into the solution vector: */
 	for(int i=0;i<numnodes;i++){
-		values[i]=solution[doflist[i]];
+		newdamage[i]=solution[doflist[i]];
 		/*Check solution*/
-		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(newdamage[i])) _error_("NaN found in solution vector");
 		/*Enforce D < max_damage and D > 0 */
-		if(values[i]>max_damage) values[i]=max_damage;
-		else if(values[i]<0.)    values[i]=0.;
+		if(newdamage[i]>max_damage) newdamage[i]=max_damage;
+		else if(newdamage[i]<0.)    newdamage[i]=0.;
 	}
 
 	/*Get all inputs and parameters*/
-	element->AddInput(DamageDbarEnum,values,P1Enum);
+	basalelement->AddBasalInput(DamageDEnum,newdamage,P1Enum);
 
 	/*Free ressources:*/
-	xDelete<IssmDouble>(values);
+	xDelete<IssmDouble>(newdamage);
 	xDelete<int>(doflist);
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 }/*}}}*/
 void DamageEvolutionAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
@@ -427,5 +464,5 @@
 	Input* sigma_xy_input  = element->GetInput(StressTensorxyEnum);     _assert_(sigma_xy_input);
 	Input* sigma_yy_input  = element->GetInput(StressTensoryyEnum);     _assert_(sigma_yy_input);
-	Input* damage_input    = element->GetInput(DamageDbarEnum); _assert_(damage_input);
+	Input* damage_input    = element->GetInput(DamageDEnum); _assert_(damage_input);
 
 	/*retrieve the desired type of equivalent stress*/
@@ -446,11 +483,10 @@
 		s_yy=sigma_yy/(1.-damage);
 
-		if(equivstress==1){ /* von Mises */
-			J2s=1./sqrt(2.)*sqrt(s_xx*s_xx + s_yy*s_yy + s_xy*s_xy);
-			Chi=sqrt(3.0)*J2s;
+		if(equivstress==0){ /* von Mises */
+			Chi=sqrt(s_xx*s_xx - s_xx*s_yy + s_yy*s_yy + 3*s_xy*s_xy);
 		}
 		Psi=Chi-stress_threshold;
 		PosPsi=max(Psi,0.);
-		NegPsi=max(-Psi,0.);
+		NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
 
 		f[iv]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1.-damage),-c3);
Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h	(revision 17390)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h	(revision 17391)
@@ -26,6 +26,6 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
 		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
 		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 17390)
+++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 17391)
@@ -713,5 +713,5 @@
 
 	//First recover damage  using the element: */
-	element->GetInputValue(&damage,node,DamageDbarEnum);
+	element->GetInputValue(&damage,node,DamageDEnum);
 
 	//Recover our data:
Index: /issm/trunk-jpl/src/c/cores/damage_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/damage_core.cpp	(revision 17390)
+++ /issm/trunk-jpl/src/c/cores/damage_core.cpp	(revision 17391)
@@ -19,5 +19,4 @@
 	char   **requested_outputs = NULL;
 
-	if(VerboseSolution()) _printf0_("   computing damage\n");
 	
 	//first recover parameters common to all solutions
@@ -33,4 +32,5 @@
 	}
 
+	if(VerboseSolution()) _printf0_("   computing damage\n");
 	femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
 	solutionsequence_linear(femmodel);