Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp	(revision 17171)
+++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp	(revision 17172)
@@ -101,8 +101,266 @@
 }/*}}}*/
 ElementMatrix* FreeSurfaceBaseAnalysis::CreateKMatrix(Element* element){/*{{{*/
-	_error_("not implemented yet");
+
+	/*Intermediaries*/
+	int         meshtype,dim,stabilization;
+	Element*    basalelement = NULL;
+	IssmDouble *xyz_list  = NULL;
+	IssmDouble  Jdet,D_scalar,dt,h;
+	IssmDouble  vel,vx,vy;
+
+	/*Get basal element*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			dim = 2;
+			break;
+		case Mesh2DverticalEnum:
+			if(!element->IsOnBed()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			dim = 1;
+			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(NoneApproximationEnum);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    D      = xNew<IssmDouble>(dim*dim);
+
+	/*Retrieve all inputs and parameters*/
+	basalelement->GetVerticesCoordinates(&xyz_list);
+	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
+	basalelement->FindParam(&stabilization,MasstransportStabilizationEnum);
+	Input* vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
+	Input* vy_input=NULL;
+	if(dim>1) basalelement->GetInput(VyEnum); _assert_(vy_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);
+
+		vx_input->GetInputValue(&vx,gauss);
+		if(dim==2) vy_input->GetInputValue(&vy,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,element,dim,xyz_list,gauss);
+		GetBprime(Bprime,element,dim,xyz_list,gauss);
+
+		D_scalar=dt*gauss->weight*Jdet;
+
+		for(int i=0;i<dim*dim;i++) D[i]=0.;
+		D[0] = D_scalar*vx;
+		if(dim==2) D[1*dim+1] = D_scalar*vy;
+		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){
+			/*Streamline upwinding*/
+			if(dim==1){
+			 vel=fabs(vx)+1.e-8;
+			 D[0] = h/(2.*vel)*vx;
+			}
+			else{
+			 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*/
+			if(dim==1){
+				vx_input->GetInputAverage(&vx);
+				D[0]=h/2.*fabs(vx);
+			}
+			else{
+				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);
+			}
+		}
+		if(stabilization==1 || stabilization==2){
+			for(int i=0;i<dim*dim;i++) D[i]=D_scalar*D[i];
+			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* FreeSurfaceBaseAnalysis::CreatePVector(Element* element){/*{{{*/
-_error_("not implemented yet");
+	/*Intermediaries*/
+	int         meshtype,dim;
+	IssmDouble  Jdet,dt;
+	IssmDouble  mb,mb_correction,bed,vz;
+	Element*    basalelement = NULL;
+	IssmDouble *xyz_list  = NULL;
+
+	/*Get basal element*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			dim = 2;
+			break;
+		case Mesh2DverticalEnum:
+			if(!element->IsOnBed()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			dim = 1;
+			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 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);
+	Input* mb_input            = basalelement->GetInput(BasalforcingsMeltingRateEnum);   _assert_(mb_input);
+	Input* mb_correction_input = basalelement->GetInput(BasalforcingsMeltingRateCorrectionEnum);
+	Input* bed_input           = basalelement->GetInput(BedEnum);                        _assert_(bed_input);
+	Input* vz_input      = NULL;
+	switch(dim){
+		case 1: vz_input = basalelement->GetInput(VyEnum); _assert_(vz_input); break;
+		case 2: vz_input = basalelement->GetInput(VzEnum); _assert_(vz_input); break;
+		default: _error_("not implemented");
+	}
+
+	/*Initialize mb_correction to 0, do not forget!:*/
+	/* 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);
+
+		vz_input->GetInputValue(&vz,gauss);
+		mb_input->GetInputValue(&mb,gauss);
+		bed_input->GetInputValue(&bed,gauss);
+		if(mb_correction_input)
+		 mb_correction_input->GetInputValue(&mb_correction,gauss);
+		else
+		 mb_correction=0.;
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(bed+dt*(mb-mb_correction) + dt*vz)*basis[i];
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+	return pe;
+
+}/*}}}*/
+void FreeSurfaceBaseAnalysis::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
+	 * by: 
+	 *       Bi=[ N ]
+	 *          [ N ]
+	 * where N is the finiteelement function for node i.
+	 *
+	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	element->NodalFunctions(basis,gauss);
+
+	/*Build B: */
+	for(int i=0;i<numnodes;i++){
+		for(int j=0;j<dim;i++){
+			B[numnodes*j+i] = basis[i];
+		}
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
+}/*}}}*/
+void FreeSurfaceBaseAnalysis::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
+	 * by: 
+	 *       Bi_prime=[ dN/dx ]
+	 *                [ dN/dy ]
+	 * where N is the finiteelement function for node i.
+	 *
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+	/*Build B': */
+	for(int i=0;i<numnodes;i++){
+		for(int j=0;j<dim;j++){
+			Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
+		}
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
+
 }/*}}}*/
 void FreeSurfaceBaseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h	(revision 17171)
+++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h	(revision 17172)
@@ -26,4 +26,6 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
+		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/analyses/FreeSurfaceTopAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp	(revision 17171)
+++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp	(revision 17172)
@@ -110,8 +110,261 @@
 }/*}}}*/
 ElementMatrix* FreeSurfaceTopAnalysis::CreateKMatrix(Element* element){/*{{{*/
-	_error_("not implemented yet");
+
+	/*Intermediaries*/
+	int         meshtype,dim,stabilization;
+	Element*    topelement = NULL;
+	IssmDouble *xyz_list  = NULL;
+	IssmDouble  Jdet,D_scalar,dt,h;
+	IssmDouble  vel,vx,vy;
+
+	/*Get top element*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			topelement = element;
+			dim = 2;
+			break;
+		case Mesh2DverticalEnum:
+			if(!element->IsOnSurface()) return NULL;
+			topelement = element->SpawnTopElement();
+			dim = 1;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnSurface()) return NULL;
+			topelement = element->SpawnTopElement();
+			dim = 2;
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = topelement->GetNumberOfNodes();
+
+	/*Initialize Element vector*/
+	ElementMatrix* Ke     = topelement->NewElementMatrix(NoneApproximationEnum);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    D      = xNew<IssmDouble>(dim*dim);
+
+	/*Retrieve all inputs and parameters*/
+	topelement->GetVerticesCoordinates(&xyz_list);
+	topelement->FindParam(&dt,TimesteppingTimeStepEnum);
+	topelement->FindParam(&stabilization,MasstransportStabilizationEnum);
+	Input* vx_input=topelement->GetInput(VxEnum); _assert_(vx_input);
+	Input* vy_input=NULL;
+	if(dim>1) topelement->GetInput(VyEnum); _assert_(vy_input);
+	h = topelement->CharacteristicLength();
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=topelement->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		topelement->NodalFunctions(basis,gauss);
+
+		vx_input->GetInputValue(&vx,gauss);
+		if(dim==2) vy_input->GetInputValue(&vy,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,element,dim,xyz_list,gauss);
+		GetBprime(Bprime,element,dim,xyz_list,gauss);
+
+		D_scalar=dt*gauss->weight*Jdet;
+
+		for(int i=0;i<dim*dim;i++) D[i]=0.;
+		D[0] = D_scalar*vx;
+		if(dim==2) D[1*dim+1] = D_scalar*vy;
+		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){
+			/*Streamline upwinding*/
+			if(dim==1){
+				vel=fabs(vx)+1.e-8;
+				D[0] = h/(2.*vel)*vx;
+			}
+			else{
+				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*/
+			if(dim==1){
+				vx_input->GetInputAverage(&vx);
+				D[0]=h/2.*fabs(vx);
+			}
+			else{
+				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);
+			}
+		}
+		if(stabilization==1 || stabilization==2){
+			for(int i=0;i<dim*dim;i++) D[i]=D_scalar*D[i];
+			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){topelement->DeleteMaterials(); delete topelement;};
+	return Ke;
 }/*}}}*/
 ElementVector* FreeSurfaceTopAnalysis::CreatePVector(Element* element){/*{{{*/
-_error_("not implemented yet");
+	/*Intermediaries*/
+	int         meshtype,dim;
+	IssmDouble  Jdet,dt;
+	IssmDouble  ms,surface,vz;
+	Element*    topelement = NULL;
+	IssmDouble *xyz_list  = NULL;
+
+	/*Get top element*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			topelement = element;
+			dim = 2;
+			break;
+		case Mesh2DverticalEnum:
+			if(!element->IsOnSurface()) return NULL;
+			topelement = element->SpawnTopElement();
+			dim = 1;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnSurface()) return NULL;
+			topelement = element->SpawnTopElement();
+			dim = 2;
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = topelement->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = topelement->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	topelement->GetVerticesCoordinates(&xyz_list);
+	topelement->FindParam(&dt,TimesteppingTimeStepEnum);
+	Input* ms_input      = topelement->GetInput(SurfaceforcingsMassBalanceEnum);  _assert_(ms_input);
+	Input* surface_input = topelement->GetInput(SurfaceEnum);                     _assert_(surface_input);
+	Input* vz_input      = NULL;
+	switch(dim){
+		case 1: vz_input = topelement->GetInput(VyEnum); _assert_(vz_input); break;
+		case 2: vz_input = topelement->GetInput(VzEnum); _assert_(vz_input); break;
+		default: _error_("not implemented");
+	}
+
+	/*Initialize mb_correction to 0, do not forget!:*/
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=topelement->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		topelement->NodalFunctions(basis,gauss);
+
+		ms_input->GetInputValue(&ms,gauss);
+		vz_input->GetInputValue(&vz,gauss);
+		surface_input->GetInputValue(&surface,gauss);
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(surface + dt*ms + dt*vz)*basis[i];
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	if(meshtype!=Mesh2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;};
+	return pe;
+
+}/*}}}*/
+void FreeSurfaceTopAnalysis::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
+	 * by: 
+	 *       Bi=[ N ]
+	 *          [ N ]
+	 * where N is the finiteelement function for node i.
+	 *
+	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	element->NodalFunctions(basis,gauss);
+
+	/*Build B: */
+	for(int i=0;i<numnodes;i++){
+		for(int j=0;j<dim;i++){
+			B[numnodes*j+i] = basis[i];
+		}
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
+}/*}}}*/
+void FreeSurfaceTopAnalysis::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
+	 * by: 
+	 *       Bi_prime=[ dN/dx ]
+	 *                [ dN/dy ]
+	 * where N is the finiteelement function for node i.
+	 *
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+	/*Build B': */
+	for(int i=0;i<numnodes;i++){
+		for(int j=0;j<dim;j++){
+			Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
+		}
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
+
 }/*}}}*/
 void FreeSurfaceTopAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h	(revision 17171)
+++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h	(revision 17172)
@@ -26,4 +26,6 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
+		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/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17171)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17172)
@@ -208,4 +208,5 @@
 		virtual void   SmbGradients(void)=0;
 	   virtual Element*   SpawnBasalElement(void)=0;
+		virtual Element*   SpawnTopElement(void)=0;
 		virtual void   ReduceMatrices(ElementMatrix* Ke,ElementVector* pe)=0;
 		virtual void   ResetCoordinateSystem()=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17171)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17172)
@@ -2605,4 +2605,14 @@
 	this->inputs->DeleteInput(VxAverageEnum);
 	this->inputs->DeleteInput(VyAverageEnum);
+
+	return tria;
+}
+/*}}}*/
+/*FUNCTION Penta::SpawnTopElement{{{*/
+Element*  Penta::SpawnTopElement(void){
+
+	_assert_(this->IsOnSurface());
+
+	Tria* tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1.
 
 	return tria;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17171)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17172)
@@ -121,4 +121,5 @@
 		void   SetTemporaryElementType(int element_type_in);
 	   Element* SpawnBasalElement(void);
+		Element* SpawnTopElement(void);
 		IssmDouble SurfaceArea(void);
 		void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17171)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17172)
@@ -121,4 +121,5 @@
 		int         NumberofNodesPressure(void){_error_("not implemented yet");};
 	   Element*    SpawnBasalElement(void){_error_("not implemented yet");};
+		Element*    SpawnTopElement(void){_error_("not implemented yet");};
 		IssmDouble  StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
 		IssmDouble  PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17171)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17172)
@@ -2261,4 +2261,23 @@
 }
 /*}}}*/
+/*FUNCTION Tria::SpawnTopElement{{{*/
+Element*  Tria::SpawnTopElement(void){
+
+	int index1,index2;
+	int meshtype;
+
+	this->parameters->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			return this;
+		case Mesh2DverticalEnum:
+			_assert_(HasEdgeOnSurface());
+			this->EdgeOnSurfaceIndices(&index1,&index2);
+			return SpawnSeg(index1,index2);
+		default:
+			_error_("not implemented yet");
+	}
+}
+/*}}}*/
 /*FUNCTION Tria::SurfaceArea {{{*/
 IssmDouble Tria::SurfaceArea(void){
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17171)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17172)
@@ -119,4 +119,5 @@
 		void	      SmbGradients();
 	   Element*    SpawnBasalElement(void);
+		Element*    SpawnTopElement(void);
 		int         VelocityInterpolation();
 		IssmDouble  PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/shared/Exceptions/Exceptions.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Exceptions/Exceptions.cpp	(revision 17171)
+++ /issm/trunk-jpl/src/c/shared/Exceptions/Exceptions.cpp	(revision 17172)
@@ -45,4 +45,5 @@
 
 	file_line= what_line;
+	this->Report();
 
 }/*}}}*/