Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18174)
+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18175)
@@ -186,18 +186,23 @@
 	/*Compute element divergence*/
 
 	/*Intermediaries*/
+	int        dim;
 	IssmDouble Jdet;
 	IssmDouble divergence=0.;
 	IssmDouble dvx[3],dvy[3],dvz[3];
 	IssmDouble *xyz_list = NULL;
 
 	/*Get inputs and parameters*/
-	Input* vx_input        = this->GetInput(VxEnum); _assert_(vx_input);
-	Input* vy_input        = this->GetInput(VyEnum); _assert_(vy_input);
-	Input* vz_input        = this->GetInput(VzEnum); _assert_(vz_input);
+	this->FindParam(&dim,DomainDimensionEnum);
+	Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input);
+	Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input);
+	Input* vz_input = NULL;
+	if(dim==3){
+		vz_input = this->GetInput(VzEnum); _assert_(vz_input);
+	}
 	this->GetVerticesCoordinates(&xyz_list);
 
-	Gauss* gauss=this->NewGaussBase(5);
+	Gauss* gauss=this->NewGauss(5);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 		gauss->GaussPoint(ig);
 		this->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -205,9 +210,14 @@
 		/*Get strain rate assuming that epsilon has been allocated*/
 		vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
 		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
-		vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
+		if(dim==2){
+			divergence += (dvx[0]+dvy[1])*gauss->weight*Jdet;
+		}
+		else{
+			vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
+			divergence += (dvx[0]+dvy[1]+dvz[2])*gauss->weight*Jdet;
+		}
 
-		divergence += (dvx[0]+dvy[1]+dvz[2])*gauss->weight*Jdet;
 	}
 
 	/*Clean up and return*/