Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5809)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5810)
@@ -809,5 +809,5 @@
 			CreatePVectorDiagnosticPattyn( pg);
 			CreatePVectorDiagnosticStokes( pg);
-			//CreatePVectorDiagnosticPattynStokes( pg);
+			CreatePVectorCouplingPattynStokes( pg);
 		}
 		else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation));
@@ -3160,4 +3160,146 @@
 }
 /*}}}*/
+/*FUNCTION Penta::CreatePVectorCouplingPattynStokes {{{1*/
+void Penta::CreatePVectorCouplingPattynStokes( Vec pg){
+
+	/*indexing: */
+	int i,j;
+	const int numdof=NUMVERTICES*NDOF4;
+	int*      doflist=NULL;
+
+	/*parameters: */
+	double		   xyz_list_tria[NUMVERTICES2D][3];
+	double         xyz_list[NUMVERTICES][3];
+	double		   bed_normal[3];
+
+	/* gaussian points: */
+	int     ig;
+	GaussPenta *gauss=NULL;
+
+	double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
+	double  viscosity, w, alpha2_gauss;
+	double  dw[3];
+
+	/*matrices: */
+	double     Pe_gaussian[24]={0.0}; //for the six nodes
+	double     l1l6[6]; //for the six nodes of the penta
+	double     dh1dh6[3][6]; //for the six nodes of the penta
+	double     Jdet;
+	double     Jdet2d;
+
+	Tria*     tria=NULL;
+	Friction* friction=NULL;
+
+	/*parameters: */
+	double stokesreconditioning;
+	int    approximation;
+	int    analysis_type;
+
+	/*retrieve inputs :*/
+	inputs->GetParameterValue(&approximation,ApproximationEnum);
+
+	/*retrieve some parameters: */
+	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+	this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
+
+	/*If on water or not Stokes, skip load: */
+	if(IsOnWater() || approximation!=PattynStokesApproximationEnum) return;
+
+	/* Get node coordinates and dof list: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+	GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
+
+	/*Retrieve all inputs we will be needing: */
+	Input* vx_input=inputs->GetInput(VxEnum);               ISSMASSERT(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);               ISSMASSERT(vy_input);
+	Input* vz_input=inputs->GetInput(VzEnum);               ISSMASSERT(vz_input);
+	Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);   ISSMASSERT(vzpattyn_input);
+
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussPenta(5,5);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		/*Compute strain rate and viscosity: */
+		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
+		matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+
+		/* Get nodal functions */
+		GetNodalFunctionsP1(&l1l6[0], gauss);
+		GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
+
+		/*Compute the derivative of VzPattyn*/
+		vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
+
+		/* Build gaussian vector */
+		for(i=0;i<NUMVERTICES;i++){
+			Pe_gaussian[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i]/2;
+			Pe_gaussian[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i]/2;
+			Pe_gaussian[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dh1dh6[0][i]+dw[1]*dh1dh6[1][i]+dw[2]*dh1dh6[2][i])/2;
+			Pe_gaussian[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i];
+		}
+	}
+	delete gauss;
+
+	/*Deal with 2d friction at the bedrock interface: */
+	if(IsOnBed() && !IsOnShelf()){
+
+		for(i=0;i<NUMVERTICES2D;i++){
+			for(j=0;j<3;j++){
+				xyz_list_tria[i][j]=xyz_list[i][j];
+			}
+		}
+
+		/*build friction object, used later on: */
+		friction=new Friction("3d",inputs,matpar,analysis_type);
+
+		/* Start looping on the number of gauss 2d (nodes on the bedrock) */
+		gauss=new GaussPenta(0,1,2,2);
+		for(ig=gauss->begin();ig<gauss->end();ig++){
+
+			gauss->GaussPoint(ig);
+
+			/*Get the Jacobian determinant */
+			GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
+
+			/*Get L matrix : */
+			GetNodalFunctionsP1(l1l6, gauss);
+
+			/*Get normal vecyor to the bed */
+			BedNormal(&bed_normal[0],xyz_list_tria);
+
+			/*Get Viscosity*/
+			this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
+			matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+
+			/*Get friction*/
+			friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
+
+			/*Get w and its derivatives*/
+			vzpattyn_input->GetParameterValue(&w, gauss);
+			vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
+
+			for(i=0;i<NUMVERTICES2D;i++){
+				Pe_gaussian[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+viscosity*dw[2]*bed_normal[0])*l1l6[i];
+				Pe_gaussian[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+viscosity*dw[2]*bed_normal[1])*l1l6[i];
+				Pe_gaussian[i*NDOF4+2]+=Jdet2d*gauss->weight*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*l1l6[i];
+			}
+		}
+		/*Free ressources:*/
+		delete gauss;
+	} //if ( (IsOnBed()==1) && (IsOnShelf()==1))
+
+	/*Add Pe_reduced to global vector pg: */
+	VecSetValues(pg,numdof,doflist,(const double*)Pe_gaussian,ADD_VALUES);
+
+	/*Free ressources:*/
+	xfree((void**)&doflist);
+
+}
+/*}}}*/
 /*FUNCTION Penta::CreatePVectorDiagnosticHutter{{{1*/
 void  Penta::CreatePVectorDiagnosticHutter(Vec pg){
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5809)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5810)
@@ -135,4 +135,5 @@
 		void	  CreatePVectorAdjointMacAyeal( Vec pg);
 		void	  CreatePVectorAdjointPattyn( Vec pg);
+		void	  CreatePVectorCouplingPattynStokes( Vec pg);
 		void	  CreatePVectorDiagnosticHutter( Vec pg);
 		void	  CreatePVectorDiagnosticMacAyeal( Vec pg);