Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 1210)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 1211)
@@ -1225,9 +1225,9 @@
 	
 	/*Bail out if this element does not touch the bedrock: */
-	if (!onbed){
-		return;
-	}
-	else{
-		
+	if (!onbed) return;
+
+	if (sub_analysis_type==HorizAnalysisEnum()){
+
+		/*MacAyeal or Pattyn*/
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 		tria->GradjDrag( grad_g,inputs,analysis_type,sub_analysis_type);
@@ -1235,4 +1235,13 @@
 		return;
 	}
+	else if (sub_analysis_type==StokesAnalysisEnum()){
+
+		/*Stokes*/
+		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+		tria->GradjDragStokes( grad_g,inputs,analysis_type,sub_analysis_type);
+		delete tria;
+		return;
+	}
+	else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
 }
 
Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 1210)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 1211)
@@ -2266,5 +2266,4 @@
 
 	/* parameters: */
-	double  viscosity2;
 	double  dvx[NDOF2];
 	double  dvy[NDOF2]; 
@@ -2430,4 +2429,243 @@
 
 #undef __FUNCT__ 
+#define __FUNCT__ "Tria::GradjDragStokes"
+void  Tria::GradjDragStokes(Vec grad_g,void* vinputs,int analysis_type,int sub_analysis_type){
+
+	int i;
+
+	/* node data: */
+	const int    numgrids=3;
+	const int    NDOF2=2;
+	double       xyz_list[numgrids][3];
+	int          doflist1[numgrids];
+	int          numberofdofspernode;
+
+	/* grid data: */
+	double vx_list[numgrids];
+	double vy_list[numgrids];
+	double vz_list[numgrids];
+	double vxvyvz_list[numgrids][3];
+	double adjx_list[numgrids];
+	double adjy_list[numgrids];
+	double adjz_list[numgrids];
+	double adjxyz_list[numgrids][3];
+
+	double drag;
+	int    dofs1[1]={0};
+	int    dofs3[3]={0,1,2};
+
+	/* gaussian points: */
+	int     num_gauss,ig;
+	double* first_gauss_area_coord  =  NULL;
+	double* second_gauss_area_coord =  NULL;
+	double* third_gauss_area_coord  =  NULL;
+	double* gauss_weights           =  NULL;
+	double  gauss_weight;
+	double  gauss_l1l2l3[3];
+
+	/* parameters: */
+	double  vx,vy,vz;
+	double  lambda,mu,xi;
+	double  bed,thickness,Neff;
+	double  surface_normal[3];
+	double  bed_normal[3];
+
+	/*drag: */
+	double  pcoeff,qcoeff;
+	double alpha_complement_list[numgrids];
+	double alpha_complement;
+
+	/*element vector at the gaussian points: */
+	double  grade_g[numgrids];
+	double  grade_g_gaussian[numgrids];
+
+	/* Jacobian: */
+	double Jdet;
+
+	/*nodal functions: */
+	double l1l2l3[3];
+
+	/* strain rate: */
+	double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+
+	ParameterInputs* inputs=NULL;
+
+	/*recover pointers: */
+	inputs=(ParameterInputs*)vinputs;
+
+	/* Get node coordinates and dof list: */
+	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetDofList1(&doflist1[0]);
+
+	/* Set grade_g to 0: */
+	for(i=0;i<numgrids;i++) grade_g[i]=0.0;
+
+	/* recover input parameters: */
+	inputs->Recover("drag",&k[0],1,dofs1,numgrids,(void**)nodes);
+	inputs->Recover("bed",&b[0],1,dofs1,numgrids,(void**)nodes);
+	inputs->Recover("thickness",&h[0],1,dofs1,numgrids,(void**)nodes);
+	if(!inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs3,numgrids,(void**)nodes)){
+		throw ErrorException(__FUNCT__,"missing velocity input parameter");
+	}
+	if(!inputs->Recover("adjoint",&adjxyz_list[0][0],3,dofs3,numgrids,(void**)nodes)){
+		throw ErrorException(__FUNCT__,"missing adjoint input parameter");
+	}
+
+	/*Initialize parameter lists: */
+	for(i=0;i<numgrids;i++){
+		vx_list[i]=vxvyvz_list[i][0];
+		vy_list[i]=vxvyvz_list[i][1];
+		vz_list[i]=vxvyvz_list[i][2];
+		adjx_list[i]=adjxyz_list[i][0];
+		adjy_list[i]=adjxyz_list[i][1];
+		adjz_list[i]=adjxyz_list[i][2];
+	}
+
+	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
+
+#ifdef _DEBUGELEMENTS_
+	if(my_rank==RANK && id==ELID){ 
+		printf("   gaussian points: \n");
+		for(i=0;i<num_gauss;i++){
+			printf("    %g %g %g : %g\n",first_gauss_area_coord[i],second_gauss_area_coord[i],third_gauss_area_coord[i],gauss_weights[i]);
+		}
+	}
+#endif
+
+	/* Start  looping on the number of gaussian points: */
+	for (ig=0; ig<num_gauss; ig++){
+		/*Pick up the gaussian point: */
+		gauss_weight=*(gauss_weights+ig);
+		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+		/*Build alpha_complement_list: */
+		if (!shelf && (friction_type==2)){
+
+			/*Allocate friction object: */
+			Friction* friction=NewFriction();
+
+			/*Initialize all fields: */
+			friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
+			strcpy(friction->element_type,"2d");
+
+			friction->gravity=matpar->GetG();
+			friction->rho_ice=matpar->GetRhoIce();
+			friction->rho_water=matpar->GetRhoWater();
+			friction->K=&k[0];
+			friction->bed=&b[0];
+			friction->thickness=&h[0];
+			friction->velocities=&vxvyvz_list[0][0];
+			friction->p=p;
+			friction->q=q;
+
+			if(friction->p!=1) throw ErrorException(__FUNCT__,"non-linear friction not supported yet in control methods!");
+
+			/*Compute alpha complement: */
+			FrictionGetAlphaComplement(&alpha_complement_list[0],friction);
+
+			/*Erase friction object: */
+			DeleteFriction(&friction);
+		}
+		else{
+			alpha_complement_list[0]=0;
+			alpha_complement_list[1]=0;
+			alpha_complement_list[2]=0;
+		}
+
+		/*Recover alpha_complement and k: */
+		GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3);
+		GetParameterValue(&drag, &k[0],gauss_l1l2l3);
+#ifdef _DEBUG_ 
+		printf("Drag complement: %20.20lf Drag: %20.20lf\n",alpha_complement,drag);
+#endif
+
+		/*recover lambda mu and xi: */
+		GetParameterValue(&lambda, &adjx_list[0],gauss_l1l2l3);
+		GetParameterValue(&mu, &adjy_list[0],gauss_l1l2l3);
+		GetParameterValue(&xi, &adjz_list[0],gauss_l1l2l3);
+#ifdef _DEBUG_ 
+		printf("Adjoint vector %20.20lf %20.20lf\n",lambda,mu);
+#endif
+
+		/*recover vx vy and vz: */
+		GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
+		GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
+		GetParameterValue(&vz, &vz_list[0],gauss_l1l2l3);
+#ifdef _DEBUG_ 
+		printf("Velocity vector %20.20lf %20.20lf\n",vx,vy);
+
+		/*Get normal vecyor to the bed */
+		SurfaceNormal(&surface_normal[0],xyz_list);
+
+		bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
+		bed_normal[1]=-surface_normal[1];
+		bed_normal[2]=-surface_normal[2];
+#endif
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+#ifdef _DEBUG_ 
+		printf("Element id %i Jacobian determinant: %lf\n",TriaElementGetID(this),Jdet);
+#endif
+
+		/* Get nodal functions value at gaussian point:*/
+		GetNodalFunctions(l1l2l3, gauss_l1l2l3);
+
+		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
+		for (i=0;i<numgrids;i++){
+			grade_g_gaussian[i]=(
+						-lambda*(2*drag*alpha_complement*(vx - vz*bed_normal[0]*bed_normal[2]))
+						-mu    *(2*drag*alpha_complement*(vy - vz*bed_normal[1]*bed_normal[2]))
+						-xi    *(2*drag*alpha_complement*(-vx*bed_normal[0]*bed_normal[2]-vy*bed_normal[1]*bed_normal[2]))
+						)*Jdet*gauss_weight*l1l2l3[i]; 
+		}
+
+		/*Add gradje_g_gaussian vector to gradje_g: */
+		for( i=0; i<numgrids; i++)grade_g[i]+=grade_g_gaussian[i];
+	}
+
+	/*Add grade_g to global vector grad_g: */
+	VecSetValues(grad_g,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
+
+cleanup_and_return: 
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&gauss_weights);
+
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::SurfaceNormal"
+
+void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
+
+	int i;
+	double v13[3];
+	double v23[3];
+	double normal[3];
+	double normal_norm;
+
+	for (i=0;i<3;i++){
+		v13[i]=xyz_list[0][i]-xyz_list[2][i];
+		v23[i]=xyz_list[1][i]-xyz_list[2][i];
+	}
+
+	normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
+	normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
+	normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
+
+	normal_norm=sqrt( pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2) );
+
+	*(surface_normal)=normal[0]/normal_norm;
+	*(surface_normal+1)=normal[1]/normal_norm;
+	*(surface_normal+2)=normal[2]/normal_norm;
+
+}
+
+#undef __FUNCT__ 
 #define __FUNCT__ "Tria::GradjB"
 void  Tria::GradjB(Vec grad_g,void* vinputs,int analysis_type,int sub_analysis_type){
Index: /issm/trunk/src/c/objects/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Tria.h	(revision 1210)
+++ /issm/trunk/src/c/objects/Tria.h	(revision 1211)
@@ -96,4 +96,6 @@
 		void  Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type);
 		void  GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type);
+		void  GradjDragStokes(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type);
+		void  SurfaceNormal(double* surface_normal, double xyz_list[3][3]);
 		void  GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type);
 		double Misfit(void* inputs,int analysis_type,int sub_analysis_type);