Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 16898)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 16899)
@@ -138,7 +138,4 @@
 /*Finite Element Analysis*/
 ElementMatrix* HydrologyDCInefficientAnalysis::CreateKMatrix(Element* element){/*{{{*/
-	_error_("not implemented yet");
-}/*}}}*/
-ElementVector* HydrologyDCInefficientAnalysis::CreatePVector(Element* element){/*{{{*/
 
 	/*Intermediaries*/
@@ -160,4 +157,79 @@
 
 	/*Intermediaries */
+	IssmDouble  D_scalar,Jdet,dt;
+	IssmDouble *xyz_list  = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = basalelement->GetNumberOfNodes();
+
+	/*Initialize Element vector*/
+	ElementMatrix* Ke     = basalelement->NewElementMatrix();
+	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble     D[2][2]={0.};
+
+	/*Retrieve all inputs and parameters*/
+	basalelement->GetVerticesCoordinates(&xyz_list);
+	IssmDouble sediment_storing       = SedimentStoring(basalelement);
+	IssmDouble sediment_transmitivity = basalelement->GetMaterialParameter(HydrologydcSedimentTransmitivityEnum);
+	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
+
+	/* 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);
+
+		/*Diffusivity*/
+		D_scalar=sediment_transmitivity*gauss->weight*Jdet;
+		if(dt!=0.) D_scalar=D_scalar*dt;
+		D[0][0]=D_scalar;
+		D[1][1]=D_scalar;
+		GetB(B,element,xyz_list,gauss); 
+		TripleMultiply(B,2,numnodes,1,
+					&D[0][0],2,2,0,
+					B,2,numnodes,0,
+					&Ke->values[0],1);
+
+		/*Transient*/
+		if(dt!=0.){
+			basalelement->NodalFunctions(&basis[0],gauss);
+			D_scalar=sediment_storing*gauss->weight*Jdet;
+			TripleMultiply(basis,numnodes,1,0,
+						&D_scalar,1,1,0,
+						basis,1,numnodes,0,
+						&Ke->values[0],1);
+		}
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+	return Ke;
+}/*}}}*/
+ElementVector* HydrologyDCInefficientAnalysis::CreatePVector(Element* element){/*{{{*/
+
+	/*Intermediaries*/
+	int      meshtype;
+	Element* basalelement;
+
+	/*Get basal element*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnBed()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Intermediaries */
 	IssmDouble dt,scalar,water_head;
 	IssmDouble water_load,transfer;
@@ -210,4 +282,31 @@
 	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return pe;
+}/*}}}*/
+void HydrologyDCInefficientAnalysis::GetB(IssmDouble* B,Element* element,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=[ 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++){
+		B[numnodes*0+i] = dbasis[0*numnodes+i];
+		B[numnodes*1+i] = dbasis[1*numnodes+i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
 void HydrologyDCInefficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h	(revision 16898)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h	(revision 16899)
@@ -27,4 +27,5 @@
 
 		/*Intermediaries*/
+		void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 		IssmDouble SedimentStoring(Element* element);
 		void ElementizeEplMask(FemModel* femmodel);
Index: /issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp	(revision 16898)
+++ /issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp	(revision 16899)
@@ -71,5 +71,61 @@
 /*Finite Element Analysis*/
 ElementMatrix* L2ProjectionEPLAnalysis::CreateKMatrix(Element* element){/*{{{*/
-	_error_("not implemented yet");
+
+	/*Intermediaries*/
+	int      meshtype;
+	Element* basalelement;
+
+	/*Get basal element*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			break;
+		case Mesh2DverticalEnum:
+			if(!element->IsOnBed()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnBed()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Intermediaries */
+	IssmDouble  D,Jdet;
+	IssmDouble *xyz_list  = NULL;
+
+	/*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);
+
+	/*Retrieve all inputs and parameters*/
+	basalelement->GetVerticesCoordinates(&xyz_list);
+
+	/* 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);
+		D=gauss->weight*Jdet;
+
+		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);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+	return Ke;
 }/*}}}*/
 ElementVector* L2ProjectionEPLAnalysis::CreatePVector(Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16898)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16899)
@@ -4790,5 +4790,5 @@
 
 	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussPenta(2,3);
+	gauss=new GaussPenta(3,3);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 16898)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 16899)
@@ -260,4 +260,5 @@
 		case HydrologydcEplCompressibilityEnum:     return this->epl_compressibility;
 		case HydrologydcWaterCompressibilityEnum:   return this->water_compressibility;
+		case HydrologydcSedimentTransmitivityEnum:  return this->sediment_transmitivity;
 		case ConstantsGEnum:                        return this->g;
 		default: _error_("Enum "<<EnumToStringx(enum_in)<<" not supported yet");