Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp	(revision 17178)
+++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp	(revision 17179)
@@ -48,9 +48,163 @@
 }/*}}}*/
 ElementMatrix* ExtrudeFromTopAnalysis::CreateKMatrix(Element* element){/*{{{*/
-	_error_("not implemented yet");
-}/*}}}*/
+
+	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=CreateKMatrixVolume(element);
+	ElementMatrix* Ke2=CreateKMatrixSurface(element);
+	ElementMatrix* Ke3=CreateKMatrixBed(element);
+	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
+
+	/*clean-up and return*/
+	delete Ke1;
+	delete Ke2;
+	delete Ke3;
+	return Ke;
+}/*}}}*/
+ElementMatrix* ExtrudeFromTopAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	IssmDouble  Jdet,D;
+	IssmDouble *xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementMatrix* Ke     = element->NewElementMatrix();
+	IssmDouble*    B      = xNew<IssmDouble>(numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+
+	/* 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);
+		element->NodalFunctions(Bprime,gauss);
+		GetB(B,element,xyz_list,gauss);
+		D=gauss->weight*Jdet;
+
+		TripleMultiply(B,1,numnodes,1,
+					&D,1,1,0,
+					Bprime,1,numnodes,0,
+					&Ke->values[0],1);
+	} 
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(Bprime);
+	delete gauss;
+	return Ke;
+}
+/*}}}*/
+ElementMatrix* ExtrudeFromTopAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/
+
+	if(!element->IsOnSurface()) return NULL;
+
+	/*Intermediaries */
+	IssmDouble  Jdet,D,normal[2];
+	IssmDouble *xyz_list_top = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementMatrix* Ke    = element->NewElementMatrix();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesTop(&xyz_list_top);
+	element->NormalTop(&normal[0],xyz_list_top);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGaussTop(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
+		element->NodalFunctions(basis,gauss);
+		D = - gauss->weight*Jdet*normal[1]; 
+
+		TripleMultiply(basis,1,numnodes,1,
+					&D,1,1,0,
+					basis,1,numnodes,0,
+					&Ke->values[0],1);
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list_top);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return Ke;
+}
+/*}}}*/
+ElementMatrix* ExtrudeFromTopAnalysis::CreateKMatrixBed(Element* element){/*{{{*/
+
+	if(!element->IsOnBed()) return NULL;
+
+	/*Intermediaries */
+	IssmDouble  Jdet,D,normal[2];
+	IssmDouble *xyz_list_base = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementMatrix* Ke    = element->NewElementMatrix();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	element->NormalBase(&normal[0],xyz_list_base);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGaussBase(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+		element->NodalFunctions(basis,gauss);
+		D = - gauss->weight*Jdet*normal[1]; 
+
+		TripleMultiply(basis,1,numnodes,1,
+					&D,1,1,0,
+					basis,1,numnodes,0,
+					&Ke->values[0],1);
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list_base);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return Ke;
+}
+/*}}}*/
 ElementVector* ExtrudeFromTopAnalysis::CreatePVector(Element* element){/*{{{*/
-_error_("not implemented yet");
-}/*}}}*/
+	return NULL;
+}/*}}}*/
+void ExtrudeFromTopAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*	Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
+		where hi is the interpolation function for node i.*/
+
+	/*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++){
+		B[i] = dbasis[1*numnodes+i];  
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
+}
+/*}}}*/
 void ExtrudeFromTopAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	   _error_("not implemented yet");
Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h	(revision 17178)
+++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h	(revision 17179)
@@ -25,5 +25,9 @@
 		ElementMatrix* CreateJacobianMatrix(Element* element);
 		ElementMatrix* CreateKMatrix(Element* element);
+		ElementMatrix* CreateKMatrixVolume(Element* element);
+		ElementMatrix* CreateKMatrixSurface(Element* element);
+		ElementMatrix* CreateKMatrixBed(Element* element);
 		ElementVector* CreatePVector(Element* element);
+		void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
 		void InputUpdateFromSolution(IssmDouble* solution,Element* element);