Index: /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp	(revision 16888)
+++ /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp	(revision 16889)
@@ -60,5 +60,80 @@
 /*Finite Element Analysis*/
 ElementMatrix* BalancevelocityAnalysis::CreateKMatrix(Element* element){/*{{{*/
-	_error_("not implemented yet");
+
+	/*Intermediaries */
+	IssmDouble  dhdt,mb,ms,Jdet;
+	IssmDouble  h,gamma,thickness;
+	IssmDouble  hnx,hny,dhnx[2],dhny[2];
+	IssmDouble *xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke     = element->NewElementMatrix();
+	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    HNx    = xNew<IssmDouble>(numnodes);
+	IssmDouble*    HNy    = xNew<IssmDouble>(numnodes);
+	IssmDouble*    H      = xNew<IssmDouble>(numnodes);
+	IssmDouble*    Nx     = xNew<IssmDouble>(numnodes);
+	IssmDouble*    Ny     = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all Inputs and parameters: */
+	element->GetVerticesCoordinates(&xyz_list);
+	Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input);
+	h = element->CharacteristicLength();
+
+	/*Get vector N for all nodes and build HNx and HNy*/
+	element->GetInputListOnNodes(Nx,SurfaceSlopeXEnum);
+	element->GetInputListOnNodes(Ny,SurfaceSlopeYEnum);
+	element->GetInputListOnNodes(H,ThicknessEnum);
+	for(int i=0;i<numnodes;i++){
+		IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]+1.e-10);
+		HNx[i] = -H[i]*Nx[i]/norm;
+		HNy[i] = -H[i]*Ny[i]/norm;
+	}
+
+	/*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);
+
+		H_input->GetInputValue(&thickness,gauss);
+		if(thickness<50.) thickness=50.;
+		element->ValueP1DerivativesOnGauss(&dhnx[0],HNx,xyz_list,gauss);
+		element->ValueP1DerivativesOnGauss(&dhny[0],HNy,xyz_list,gauss);
+		element->ValueP1OnGauss(&hnx,HNx,gauss);
+		element->ValueP1OnGauss(&hny,HNy,gauss);
+
+		gamma=h/(2.*thickness+1.e-10);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
+							(basis[i]+gamma*(basis[i]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+i]*hnx + dbasis[1*numnodes+i]*hny))*
+							(basis[j]*(dhnx[0]+dhny[1])  + dbasis[0*numnodes+j]*hnx + dbasis[1*numnodes+j]*hny)
+							);
+			}
+		}
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
+	xDelete<IssmDouble>(H);
+	xDelete<IssmDouble>(Nx);
+	xDelete<IssmDouble>(Ny);
+	xDelete<IssmDouble>(HNx);
+	xDelete<IssmDouble>(HNy);
+	xDelete<IssmDouble>(B);
+	delete gauss;
+	return Ke;
 }/*}}}*/
 ElementVector* BalancevelocityAnalysis::CreatePVector(Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.cpp	(revision 16888)
+++ /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.cpp	(revision 16889)
@@ -46,5 +46,49 @@
 /*Finite Element Analysis*/
 ElementMatrix* SmoothedSurfaceSlopeXAnalysis::CreateKMatrix(Element* element){/*{{{*/
-	_error_("not implemented yet");
+
+	/* Intermediaries */
+	IssmDouble  Jdet,thickness,l=8.;
+	IssmDouble *xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke     = element->NewElementMatrix();
+	IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_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);
+		thickness_input->GetInputValue(&thickness,gauss);
+		if(thickness<50.) thickness=50.;
+
+		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
+							basis[i]*basis[j]
+							+(l*thickness)*(l*thickness)*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
+							);
+			}
+		}
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(dbasis);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list);
+	return Ke;
 }/*}}}*/
 ElementVector* SmoothedSurfaceSlopeXAnalysis::CreatePVector(Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.cpp	(revision 16888)
+++ /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.cpp	(revision 16889)
@@ -46,5 +46,49 @@
 /*Finite Element Analysis*/
 ElementMatrix* SmoothedSurfaceSlopeYAnalysis::CreateKMatrix(Element* element){/*{{{*/
-	_error_("not implemented yet");
+
+	/* Intermediaries */
+	IssmDouble  Jdet,thickness,l=8.;
+	IssmDouble *xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke     = element->NewElementMatrix();
+	IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_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);
+		thickness_input->GetInputValue(&thickness,gauss);
+		if(thickness<50.) thickness=50.;
+
+		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
+							basis[i]*basis[j]
+							+(l*thickness)*(l*thickness)*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
+							);
+			}
+		}
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(dbasis);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list);
+	return Ke;
 }/*}}}*/
 ElementVector* SmoothedSurfaceSlopeYAnalysis::CreatePVector(Element* element){/*{{{*/