Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 17330)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 17331)
@@ -261,5 +261,5 @@
 	/*Intermediaries */
 	int        stabilization;
-	int        meshtype;
+	int        meshtype,dim;
 	IssmDouble Jdet,D_scalar,dt,h;
 	IssmDouble vel,vx,vy,dvxdx,dvydy;
@@ -267,4 +267,13 @@
 	IssmDouble* xyz_list = NULL;
 
+	/*Get problem dimension*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum:   dim = 1; break;
+		case Mesh2DhorizontalEnum: dim = 2; break;
+		case Mesh3DEnum:           dim = 2; break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = element->GetNumberOfNodes();
@@ -273,7 +282,7 @@
 	ElementMatrix* Ke     = element->NewElementMatrix();
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
-	IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
-	IssmDouble     D[2][2]={0.};
+	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
 
 	/*Retrieve all inputs and parameters*/
@@ -290,5 +299,6 @@
 	else{
 		vxaverage_input=element->GetInput(VxAverageEnum); _assert_(vxaverage_input);
-		vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input);
+		if(dim==2) vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input);
+
 	}
 	h = element->CharacteristicLength();
@@ -303,7 +313,9 @@
 
 		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);
+		if(dim==2){
+			vyaverage_input->GetInputValue(&vy,gauss);
+			vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+		}
 
 		D_scalar=gauss->weight*Jdet;
@@ -313,48 +325,60 @@
 					&Ke->values[0],1);
 
-		GetB(B,element,xyz_list,gauss);
-		GetBprime(Bprime,element,xyz_list,gauss);
+		GetB(B,element,dim,xyz_list,gauss);
+		GetBprime(Bprime,element,dim,xyz_list,gauss);
 
 		dvxdx=dvx[0];
-		dvydy=dvy[1];
+		if(dim==2) 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,
+		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][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,
+		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){
-			/*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;
+			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){
 			/*SSA*/
 			vxaverage_input->GetInputAverage(&vx);
-			vyaverage_input->GetInputAverage(&vy);
-			D[0][0]=h/2.0*fabs(vx);
-			D[1][1]=h/2.0*fabs(vy);
+			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){
-			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,
+			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);
 		}
@@ -366,4 +390,5 @@
 	xDelete<IssmDouble>(B);
 	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(D);
 	delete gauss;
 	return Ke;
@@ -422,6 +447,6 @@
 
 		/*WARNING: B and Bprime are inverted compared to CG*/
-		GetB(Bprime,element,xyz_list,gauss);
-		GetBprime(B,element,xyz_list,gauss);
+		GetB(Bprime,element,2,xyz_list,gauss);
+		GetBprime(B,element,2,xyz_list,gauss);
 
 		D_scalar = - dt*gauss->weight*Jdet;
@@ -568,5 +593,5 @@
 	return pe;
 }/*}}}*/
-void MasstransportAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void MasstransportAnalysis::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
@@ -588,6 +613,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];
+		}
 	}
 
@@ -595,5 +621,5 @@
 	xDelete<IssmDouble>(basis);
 }/*}}}*/
-void MasstransportAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void MasstransportAnalysis::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
@@ -610,11 +636,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];
+		}
 	}
 
Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h	(revision 17330)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h	(revision 17331)
@@ -30,6 +30,6 @@
 		ElementVector* CreatePVectorCG(Element* element);
 		ElementVector* CreatePVectorDG(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);