Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 18469)
+++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 18470)
@@ -74,10 +74,7 @@
 	int finiteelement;
 
 	iomodel->Constant(&finiteelement,DamageElementinterpEnum);
-
-	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(1,MeshVertexonbaseEnum);
 	::CreateNodes(nodes,iomodel,DamageEvolutionAnalysisEnum,finiteelement);
-	iomodel->DeleteData(1,MeshVertexonbaseEnum);
 }/*}}}*/
 void DamageEvolutionAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
 
@@ -108,74 +105,60 @@
 
 	/* Check if ice in element */
 	if(!element->IsIceInElement()) return NULL;
-
 	/*Intermediaries*/
-	Element*    basalelement;
 	int         domaintype,dim;
 	int         stabilization;
-	IssmDouble  Jdet,dt,D_scalar,h;
-	IssmDouble  vel,vx,vy,dvxdx,dvydy,dvx[2],dvy[2];
+	IssmDouble  Jdet,dt,D_scalar,h,hx,hy,hz;
+	IssmDouble  vel,vx,vy,vz,dvxdx,dvydy,dvzdz,dvx[3],dvy[3],dvz[3];
 	IssmDouble *xyz_list  = NULL;
 
-	/*Get problem dimension and basal element*/
+	/*Get problem dimension*/
 	element->FindParam(&domaintype,DomainTypeEnum);
 	switch(domaintype){
-		case Domain2DhorizontalEnum:
-			basalelement = element;
-			dim = 2;
-			break;
-		case Domain3DEnum:
-			if(!element->IsOnBase()) return NULL;
-			basalelement = element->SpawnBasalElement();
-			dim = 2;
-			break;
-		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+		case Domain2DhorizontalEnum: dim = 2; break;
+		case Domain3DEnum:           dim = 3; break;
+		default: _error_("Not implemented yet");
 	}
 
 	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = basalelement->GetNumberOfNodes();
+	int numnodes = element->GetNumberOfNodes();
 
 	/*Initialize Element vector*/
-	ElementMatrix* Ke     = basalelement->NewElementMatrix();
+	ElementMatrix* Ke     = element->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(domaintype==Domain2DhorizontalEnum){
-		vxaverage_input=basalelement->GetInput(VxEnum); _assert_(vxaverage_input);
-		vyaverage_input=basalelement->GetInput(VyEnum); _assert_(vyaverage_input);
+	element->GetVerticesCoordinates(&xyz_list);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	element->FindParam(&stabilization,DamageStabilizationEnum);
+	Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
+	Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
+	Input* vz_input = NULL;
+	if(dim==3){
+		vz_input=element->GetInput(VzEnum); _assert_(vz_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();
 
+	if(dim==2) h=element->CharacteristicLength();
+
 	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=basalelement->NewGauss(2);
+	Gauss* gauss=element->NewGauss(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 		gauss->GaussPoint(ig);
 
-		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		basalelement->NodalFunctions(basis,gauss);
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->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);
+		vx_input->GetInputValue(&vx,gauss);
+		vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+
+		if(dim==3){
+			vz_input->GetInputValue(&vz,gauss);
+			vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
 		}
 
 		D_scalar=gauss->weight*Jdet;
@@ -184,15 +167,17 @@
 					basis,1,numnodes,0,
 					&Ke->values[0],1);
 
-		GetB(B,basalelement,dim,xyz_list,gauss);
-		GetBprime(Bprime,basalelement,dim,xyz_list,gauss);
+		GetB(B,element,dim,xyz_list,gauss);
+		GetBprime(Bprime,element,dim,xyz_list,gauss);
 
 		dvxdx=dvx[0];
-		if(dim==2) dvydy=dvy[1];
+		dvydy=dvy[1];
+		if(dim==3) dvzdz=dvz[2];
 		D_scalar=dt*gauss->weight*Jdet;
 
 		D[0*dim+0]=D_scalar*dvxdx;
-		if(dim==2) D[1*dim+1]=D_scalar*dvydy;
+		D[1*dim+1]=D_scalar*dvydy;
+		if(dim==3) D[2*dim+2]=D_scalar*dvzdz;
 
 		TripleMultiply(B,dim,numnodes,1,
 					D,dim,dim,0,
@@ -200,7 +185,8 @@
 					&Ke->values[0],1);
 
 		D[0*dim+0]=D_scalar*vx;
-		if(dim==2) D[1*dim+1]=D_scalar*vy;
+		D[1*dim+1]=D_scalar*vy;
+		if(dim==3) D[2*dim+2]=D_scalar*vz;
 
 		TripleMultiply(B,dim,numnodes,1,
 					D,dim,dim,0,
@@ -208,9 +194,17 @@
 					&Ke->values[0],1);
 
 		if(stabilization==2){
-			if(dim==1){
-				vel=fabs(vx)+1.e-8;
-				D[0]=h/(2.0*vel)*vx*vx;
+			if(dim==3){
+				vel=sqrt(vx*vx+vy*vy+vz*vz)+1.e-8;
+				D[0*dim+0]=h/(2.0*vel)*vx*vx;
+				D[1*dim+0]=h/(2.0*vel)*vy*vx;
+				D[2*dim+0]=h/(2.0*vel)*vz*vx;
+				D[0*dim+1]=h/(2.0*vel)*vx*vy;
+				D[1*dim+1]=h/(2.0*vel)*vy*vy;
+				D[2*dim+1]=h/(2.0*vel)*vy*vz;
+				D[0*dim+2]=h/(2.0*vel)*vx*vz;
+				D[1*dim+2]=h/(2.0*vel)*vy*vz;
+				D[2*dim+2]=h/(2.0*vel)*vz*vz;
 			}
 			else{
 				/*Streamline upwinding*/
@@ -222,20 +216,39 @@
 			}
 		}
 		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(dim==2){
+				vx_input->GetInputAverage(&vx);
+				vy_input->GetInputAverage(&vy);
+				D[0*dim+0]=h/2.0*fabs(vx);
+				D[1*dim+1]=h/2.0*fabs(vy);
+			}
+			else if(dim==3){ 
+				element->ElementSizes(&hx,&hy,&hz);
+				vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
+				h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
+				D[0*dim+0]=h/(2.*vel)*fabs(vx*vx);  D[0*dim+1]=h/(2.*vel)*fabs(vx*vy); D[0*dim+2]=h/(2.*vel)*fabs(vx*vz);
+				D[1*dim+0]=h/(2.*vel)*fabs(vy*vx);  D[1*dim+1]=h/(2.*vel)*fabs(vy*vy); D[1*dim+2]=h/(2.*vel)*fabs(vy*vz);
+				D[2*dim+0]=h/(2.*vel)*fabs(vz*vx);  D[2*dim+1]=h/(2.*vel)*fabs(vz*vy); D[2*dim+2]=h/(2.*vel)*fabs(vz*vz);
+			}
 		}
 		if(stabilization==1 || stabilization==2){
-			if(dim==1) D[0]=D_scalar*D[0];
-			else{
+			if(dim==2){
 				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];
 			}
-
+			else if(dim==3){
+				D[0*dim+0]=D_scalar*D[0*dim+0];
+				D[1*dim+0]=D_scalar*D[1*dim+0];
+				D[2*dim+0]=D_scalar*D[2*dim+0];
+				D[0*dim+1]=D_scalar*D[0*dim+1];
+				D[1*dim+1]=D_scalar*D[1*dim+1];
+				D[2*dim+1]=D_scalar*D[2*dim+1];
+				D[0*dim+2]=D_scalar*D[0*dim+2];
+				D[1*dim+2]=D_scalar*D[1*dim+2];
+				D[2*dim+2]=D_scalar*D[2*dim+2];
+			}
 			TripleMultiply(Bprime,dim,numnodes,1,
 						D,dim,dim,0,
 						Bprime,dim,numnodes,0,
@@ -251,7 +264,6 @@
 	xDelete<IssmDouble>(Bprime);
 	xDelete<IssmDouble>(D);
 	delete gauss;
-	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return Ke;
 }/*}}}*/
 ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/
@@ -261,53 +273,54 @@
 
 	/*Intermediaries*/
 	int      domaintype,damagelaw;
-	Element* basalelement;
 	IssmDouble  Jdet,dt;
 	IssmDouble  f,damage;
 	IssmDouble* xyz_list = NULL;
-
-	/*Get basal element*/
+	/*Get element*/
 	element->FindParam(&domaintype,DomainTypeEnum);
-	switch(domaintype){
-		case Domain2DhorizontalEnum:
-			basalelement = element;
-			break;
-		case Domain3DEnum:
-			if(!element->IsOnBase()) return NULL;
-			basalelement = element->SpawnBasalElement();
-			break;
-		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-	}
 
-
 	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = basalelement->GetNumberOfNodes();
+	int numnodes = element->GetNumberOfNodes();
 
 	/*Initialize Element vector and other vectors*/
-	ElementVector* pe    = basalelement->NewElementVector();
+	ElementVector* pe    = element->NewElementVector();
 	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
-	basalelement->GetVerticesCoordinates(&xyz_list);
-	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
-	basalelement->FindParam(&damagelaw,DamageLawEnum);
-	if(damagelaw==1 | damagelaw==2){
-		this->CreateDamageFInputPralong(basalelement);
+	element->GetVerticesCoordinates(&xyz_list);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	element->FindParam(&damagelaw,DamageLawEnum);
+	switch(damagelaw){
+		case 1:
+			this->CreateDamageFInputPralong(element);
+			break;
+		case 2:
+			this->CreateDamageFInputPralong(element);
+			break;
+		case 3:
+			this->CreateDamageFInputExp(element);
+			break;
+		default:
+			_error_("not implemented yet");
 	}
-	else if(damagelaw==3){
-		this->CreateDamageFInputExp(basalelement);
+
+	Input* damaged_input = NULL;
+	Input* damagef_input = element->GetInput(DamageFEnum); _assert_(damagef_input);
+	if(domaintype==Domain2DhorizontalEnum){
+		damaged_input = element->GetInput(DamageDbarEnum); _assert_(damaged_input);
 	}
-	Input* damaged_input = basalelement->GetInput(DamageDEnum); _assert_(damaged_input);
-	Input* damagef_input = basalelement->GetInput(DamageFEnum); _assert_(damagef_input);
+	else{
+		damaged_input = element->GetInput(DamageDEnum); _assert_(damaged_input);
+	}
 
 
 	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=basalelement->NewGauss(2);
+	Gauss* gauss=element->NewGauss(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 		gauss->GaussPoint(ig);
 
-		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		basalelement->NodalFunctions(basis,gauss);
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
 
 		damaged_input->GetInputValue(&damage,gauss);
 		damagef_input->GetInputValue(&f,gauss);
@@ -316,11 +329,9 @@
 			pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i];
 		}
 	}
-
 	/*Clean up and return*/
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(basis);
-	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	delete gauss;
 	return pe;
 }/*}}}*/
@@ -392,25 +403,18 @@
 	int domaintype;
 	IssmDouble  max_damage;
 	int			*doflist = NULL;
-	Element*   basalelement=NULL;
 
 	element->FindParam(&domaintype,DomainTypeEnum);
-	if(domaintype!=Domain2DhorizontalEnum){
-		if(!element->IsOnBase()) return;
-		basalelement=element->SpawnBasalElement();
-	}
-	else{
-		basalelement = element;
-	}
+
 	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = basalelement->GetNumberOfNodes();
+	int numnodes = element->GetNumberOfNodes();
 
 	/*Fetch dof list and allocate solution vector*/
-	basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
 	IssmDouble* newdamage = xNew<IssmDouble>(numnodes);
 
 	/*Get user-supplied max_damage: */
-	basalelement->FindParam(&max_damage,DamageMaxDamageEnum);
+	element->FindParam(&max_damage,DamageMaxDamageEnum);
 
 	/*Use the dof list to index into the solution vector: */
 	for(int i=0;i<numnodes;i++){
@@ -424,16 +428,15 @@
 
 	/*Get all inputs and parameters*/
 	if(domaintype==Domain2DhorizontalEnum){
-		basalelement->AddInput(DamageDbarEnum,newdamage,element->GetElementType());
+		element->AddInput(DamageDbarEnum,newdamage,element->GetElementType());
 	}
 	else{
-		basalelement->AddBasalInput(DamageDEnum,newdamage,element->GetElementType());
+		element->AddInput(DamageDEnum,newdamage,element->GetElementType());
 	}
 
 	/*Free ressources:*/
 	xDelete<IssmDouble>(newdamage);
 	xDelete<int>(doflist);
-	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 }/*}}}*/
 void DamageEvolutionAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	/*Default, do nothing*/
@@ -445,10 +448,10 @@
 
 	/*Intermediaries */
 	IssmDouble c1,c2,c3,healing,stress_threshold;
-	IssmDouble s_xx,s_xy,s_yy,s1,s2,stmp;
+	IssmDouble s_xx,s_xy,s_xz,s_yy,s_yz,s_zz,s1,s2,s3,stmp;
 	IssmDouble J2s,Chi,Psi,PosPsi,NegPsi;
-	IssmDouble damage,tau_xx,tau_xy,tau_yy;
-	int equivstress,domaintype,damagelaw;
+	IssmDouble damage,tau_xx,tau_xy,tau_xz,tau_yy,tau_yz,tau_zz,stressMaxPrincipal;
+	int equivstress,domaintype,damagelaw,dim;
 
 	/*Fetch number of vertices and allocate output*/
 	int numnodes = element->GetNumberOfNodes();
@@ -463,13 +466,32 @@
 	element->FindParam(&domaintype,DomainTypeEnum);
 	element->FindParam(&damagelaw,DamageLawEnum);
 
-	/*Compute stress tensor: */
+	/*Get problem dimension*/
+	switch(domaintype){
+		case Domain2DhorizontalEnum: dim = 2; break;
+		case Domain3DEnum:           dim = 3; break;
+		default: _error_("not implemented");
+	}
+	/*Compute stress tensor and Stress Max Principal: */
 	element->ComputeDeviatoricStressTensor();
-
+	if(dim==3){
+		/*Only works in 3d because the pressure is defined*/
+		element->StressMaxPrincipalCreateInput();
+	}
 	/*retrieve what we need: */
 	Input* tau_xx_input  = element->GetInput(DeviatoricStressxxEnum);     _assert_(tau_xx_input);
 	Input* tau_xy_input  = element->GetInput(DeviatoricStressxyEnum);     _assert_(tau_xy_input);
 	Input* tau_yy_input  = element->GetInput(DeviatoricStressyyEnum);     _assert_(tau_yy_input);
+	Input* tau_xz_input  = NULL;
+	Input* tau_yz_input  = NULL;
+	Input* tau_zz_input  = NULL;
+	Input* stressMaxPrincipal_input = NULL;
+	if(dim==3){
+		tau_xz_input  = element->GetInput(DeviatoricStressxzEnum);     _assert_(tau_xz_input);
+		tau_yz_input  = element->GetInput(DeviatoricStressyzEnum);     _assert_(tau_yz_input);
+		tau_zz_input  = element->GetInput(DeviatoricStresszzEnum);     _assert_(tau_zz_input);
+		stressMaxPrincipal_input = element->GetInput(StressMaxPrincipalEnum); _assert_(stressMaxPrincipal_input);
+	}
 	Input* damage_input = NULL;
 	if(domaintype==Domain2DhorizontalEnum){
 		damage_input = element->GetInput(DamageDbarEnum); 	_assert_(damage_input);
@@ -490,37 +512,66 @@
 		tau_xx_input->GetInputValue(&tau_xx,gauss);
 		tau_xy_input->GetInputValue(&tau_xy,gauss);
 		tau_yy_input->GetInputValue(&tau_yy,gauss);
-	
+		if(dim==3){
+			tau_xz_input->GetInputValue(&tau_xz,gauss);
+			tau_yz_input->GetInputValue(&tau_yz,gauss);
+			tau_zz_input->GetInputValue(&tau_zz,gauss);
+		}
 		/*Calculate effective stress components*/
 		s_xx=tau_xx/(1.-damage);
 		s_xy=tau_xy/(1.-damage);
 		s_yy=tau_yy/(1.-damage);
-
+		if(dim==3){
+			s_xz=tau_xz/(1.-damage);
+			s_yz=tau_yz/(1.-damage);
+			s_zz=tau_zz/(1.-damage);
+		}
 		/*Calculate principal effective stresses*/
-		s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
-		s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
-		if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;}
+		if(dim==2){
+			s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
+			s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
+			if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;}
 
-		if(equivstress==0){ /* von Mises */
-			Chi=sqrt(s1*s1-s1*s2+s2*s2);
-		}
-		else if(equivstress==1){ /* max principal stress */
-			Chi=s1;
-		}
-		Psi=Chi-stress_threshold;
-		NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
+			if(equivstress==0){ /* von Mises */
+				Chi=sqrt(s1*s1-s1*s2+s2*s2);
+			}
+			else if(equivstress==1){ /* max principal stress */
+				Chi=s1;
+			}
+			Psi=Chi-stress_threshold;
+			NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
 
-		if(damagelaw==1){
-			PosPsi=max(Psi,0.);
-			f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
+			if(damagelaw==1){
+				PosPsi=max(Psi,0.);
+				f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
+			}
+			else if(damagelaw==2){
+				PosPsi=max(Psi,1.);
+				f[i]= c1*(pow(log10(PosPsi),c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
+			}
+			else _error_("damage law not supported");
 		}
-		else if(damagelaw==2){
-			PosPsi=max(Psi,1.);
-			f[i]= c1*(pow(log10(PosPsi),c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
+		else{
+			if(equivstress==1){/* max principal stress */
+				stressMaxPrincipal_input->GetInputValue(&stressMaxPrincipal,gauss);
+				Chi=stressMaxPrincipal/(1.-damage);
+			}
+			else if(equivstress==0){/* von Mises */
+				Chi=sqrt(((s_xx-s_yy)*(s_xx-s_yy)+(s_yy-s_zz)*(s_yy-s_zz)+(s_zz-s_xx)*(s_zz-s_xx)+6.*(s_xy*s_xy+s_yz*s_yz+s_xz*s_xz))/2.);
+			}
+			Psi=Chi-stress_threshold;
+			NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
+			if(damagelaw==1){
+				PosPsi=max(Psi,0.);
+				f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
+			}
+			else if(damagelaw==2){
+				PosPsi=max(Psi,1.);
+				f[i]= c1*(pow(log10(PosPsi),c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
+			}
+			else _error_("damage law not supported");
 		}
-		else _error_("damage law not supported");
 	}
-
 	/*Add input*/
 	element->AddInput(DamageFEnum,f,element->GetElementType());
 	
Index: ../trunk-jpl/src/c/cores/damage_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/damage_core.cpp	(revision 18469)
+++ ../trunk-jpl/src/c/cores/damage_core.cpp	(revision 18470)
@@ -33,7 +33,7 @@
 		if(VerboseSolution()) _printf0_("   saving results\n");
 		femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
 	}
-	
+
 	/*Free resources:*/	
 	if(numoutputs){
 		for(int i=0;i<numoutputs;i++){
Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18469)
+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18470)
@@ -1407,7 +1407,7 @@
 
 		/*Initialize maximum eigne value*/
 		if(numroots>0){
-			max = x[0];
+			max = fabs(x[0]);
 		}
 		else{
 			_error_("No eigen value found");
@@ -1415,7 +1415,7 @@
 
 		/*Get max*/
 		for(int i=1;i<numroots;i++){
-			if(x[i]>max) max = x[i];
+			if(fabs(x[i])>max) max = fabs(x[i]);
 		}
 
 		maxprincipal[iv]=max;
Index: ../trunk-jpl/src/c/classes/Node.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Node.cpp	(revision 18469)
+++ ../trunk-jpl/src/c/classes/Node.cpp	(revision 18470)
@@ -93,7 +93,6 @@
 				analysis_enum==L2ProjectionBaseAnalysisEnum || 
 				analysis_enum==BalancethicknessAnalysisEnum ||
 				analysis_enum==HydrologyDCInefficientAnalysisEnum ||
-				analysis_enum==DamageEvolutionAnalysisEnum || 
 				analysis_enum==HydrologyDCEfficientAnalysisEnum ||
 				analysis_enum==LevelsetAnalysisEnum ||
 				analysis_enum==ExtrapolationAnalysisEnum
