Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 3832)
+++ /issm/trunk/src/c/Makefile.am	(revision 3833)
@@ -50,6 +50,4 @@
 					./objects/Loads/Friction.h\
 					./objects/Loads/Friction.cpp\
-					./objects/Loads/Friction2.h\
-					./objects/Loads/Friction2.cpp\
 					./objects/DakotaPlugin.h\
 					./objects/DakotaPlugin.cpp\
@@ -499,6 +497,4 @@
 					./objects/Loads/Friction.h\
 					./objects/Loads/Friction.cpp\
-					./objects/Loads/Friction2.h\
-					./objects/Loads/Friction2.cpp\
 					./objects/DakotaPlugin.h\
 					./objects/DakotaPlugin.cpp\
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 3832)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 3833)
@@ -1180,14 +1180,6 @@
 	double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
 	double  viscosity;
-	double  alpha2_list[numgrids2d];
 	double  alpha2_gauss;
-
-	double  vx_list[numgrids];
-	double  vy_list[numgrids];
-	double  vz_list[numgrids];
-	double  thickness_list[numgrids];
-	double  bed_list[numgrids];
-	double  dragcoefficient_list[numgrids];
-	double  drag_p,drag_q;
+	Friction* friction=NULL;
 
 	/*parameters: */
@@ -1281,37 +1273,6 @@
 	if((onbed==1) && (shelf==0)){
 
-		/*Build alpha2_list used by drag stiffness matrix*/
-		Friction* friction=NewFriction();
-
-		/*Initialize all fields: */
-		friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
-		strcpy(friction->element_type,"2d");
-
-		inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],6,VxEnum);
-		inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],6,VyEnum);
-		inputs->GetParameterValues(&vz_list[0],&gaussgrids[0][0],6,VzEnum);
-		inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],6,DragCoefficientEnum);
-		inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],6,BedEnum);
-		inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],6,ThicknessEnum);
-		inputs->GetParameterValue(&drag_p,DragPEnum);
-		inputs->GetParameterValue(&drag_q,DragQEnum);
-
-		friction->gravity=matpar->GetG();
-		friction->rho_ice=matpar->GetRhoIce();
-		friction->rho_water=matpar->GetRhoWater();
-		friction->K=&dragcoefficient_list[0];
-		friction->bed=&bed_list[0];
-		friction->thickness=&thickness_list[0];
-		friction->vx=&vx_list[0];
-		friction->vy=&vy_list[0];
-		friction->vz=&vz_list[0];
-		friction->p=drag_p;
-		friction->q=drag_q;
-
-		/*Compute alpha2_list: */
-		FrictionGetAlpha2(&alpha2_list[0],friction);
-
-		/*Erase friction object: */
-		DeleteFriction(&friction);
+		/*build friction object, used later on: */
+		friction=new Friction("2d",inputs,matpar);
 
 		for(i=0;i<numgrids2d;i++){
@@ -1357,5 +1318,5 @@
 
 			/*Calculate DL on gauss point */
-			tria->GetParameterValue(&alpha2_gauss,&alpha2_list[0],gauss_coord_tria);
+			friction->GetAlpha2(&alpha2_gauss, gauss_coord_tria,VxEnum,VyEnum,VzEnum);
 
 			DLStokes[0][0]=alpha2_gauss*gauss_weight*Jdet2d;
@@ -1384,4 +1345,8 @@
 			}
 		}
+	
+		/*Free ressources:*/
+		delete friction;
+
 	} //if ( (onbed==1) && (shelf==0))
 
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 3832)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 3833)
@@ -1299,5 +1299,5 @@
 
 	/*friction: */
-	Friction2* friction2=NULL;
+	Friction* friction=NULL;
 	double alpha2;
 
@@ -1324,5 +1324,5 @@
 	/*build friction object, used later on: */
 	if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!");
-	friction2=new Friction2("2d",inputs,matpar);
+	friction=new Friction("2d",inputs,matpar);
 
 	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
@@ -1339,5 +1339,5 @@
 
 		/*Friction: */
-		friction2->GetAlpha2(&alpha2, gauss_l1l2l3,VxAverageEnum,VyAverageEnum,VzAverageEnum);
+		friction->GetAlpha2(&alpha2, gauss_l1l2l3,VxAverageEnum,VyAverageEnum,VzAverageEnum);
 
 		// If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
@@ -1380,5 +1380,5 @@
 	xfree((void**)&third_gauss_area_coord);
 	xfree((void**)&gauss_weights);
-	delete friction2;
+	delete friction;
 
 }	
@@ -2902,5 +2902,5 @@
 	int    drag_type;
 	double basalfriction;
-	Friction2* friction2=NULL;
+	Friction* friction=NULL;
 	double alpha2,vx,vy;
 	double geothermalflux_value;
@@ -2936,5 +2936,5 @@
 	inputs->GetParameterValue(&drag_type,DragTypeEnum);
 	if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!");
-	friction2=new Friction2("3d",inputs,matpar);
+	friction=new Friction("3d",inputs,matpar);
 	
 	/* Ice/ocean heat exchange flux on ice shelf base */
@@ -2957,5 +2957,5 @@
 		inputs->GetParameterValue(&geothermalflux_value, &gauss_coord[0],GeothermalFluxEnum);
 	
-		friction2->GetAlpha2(&alpha2,&gauss_coord[0],VxAverageEnum,VyAverageEnum,VzAverageEnum);
+		friction->GetAlpha2(&alpha2,&gauss_coord[0],VxAverageEnum,VyAverageEnum,VzAverageEnum);
 		inputs->GetParameterValue(&vx, &gauss_coord[0],VxAverageEnum);
 		inputs->GetParameterValue(&vy, &gauss_coord[0],VyAverageEnum);
@@ -2981,5 +2981,5 @@
 	xfree((void**)&third_gauss_area_coord);
 	xfree((void**)&gauss_weights);
-	delete friction2;
+	delete friction;
 
 }
@@ -3943,16 +3943,161 @@
 
 	/* grid data: */
-	double vx_list[numgrids];
-	double vy_list[numgrids];
 	double adjx_list[numgrids];
 	double adjy_list[numgrids];
-	double thickness_list[numgrids];
-	double bed_list[numgrids];
-	double dragcoefficient_list[numgrids];
-	double drag_p;
-	double drag_q;
-	int    drag_type;
-
+
+	/* 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  dk[NDOF2]; 
+	double  vx,vy;
+	double  lambda,mu;
+	double  bed,thickness,Neff;
+	double  alpha_complement;
+	int     drag_type;
+	double  drag;
+	Friction* friction=NULL;
+
+	/*element vector at the gaussian points: */
+	double  grade_g[numgrids]={0.0};
+	double  grade_g_gaussian[numgrids];
+
+	/* Jacobian: */
+	double Jdet;
+
+	/*nodal functions: */
+	double l1l2l3[3];
+
+	/* strain rate: */
+	double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+
+	/*inputs: */
+	bool shelf;
+
+	/*parameters: */
+	double  cm_noisedmp;
+	double  cm_mindmp_slope;
+	double  cm_mindmp_value;
+	double  cm_maxdmp_value;
+	double  cm_maxdmp_slope;
+
+	/*retrieve inputs :*/
+	inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+
+	/*retrieve some parameters: */
+	this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
+	this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum);
+	this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum);
+	this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum);
+	this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum);
+
+
+	/*Get out if shelf*/
+	if(shelf)return;
+
+	/* Get node coordinates and dof list: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetDofList1(&doflist1[0]);
+
+	/*Build frictoin element, needed later: */
+	inputs->GetParameterValue(&drag_type,DragTypeEnum);
+	friction=new Friction("2d",inputs,matpar);
+
+	/* 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);
+
+	/* 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 (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum);
+		else alpha_complement=0;
+	
+		/*Recover alpha_complement and k: */
+		inputs->GetParameterValue(&drag, gauss_l1l2l3,DragCoefficientEnum);
+
+		/*recover lambda and mu: */
+		inputs->GetParameterValue(&lambda, gauss_l1l2l3,AdjointxEnum);
+		inputs->GetParameterValue(&mu, gauss_l1l2l3,AdjointyEnum);
+			
+		/*recover vx and vy: */
+		inputs->GetParameterValue(&vx, gauss_l1l2l3,VxEnum);
+		inputs->GetParameterValue(&vy, gauss_l1l2l3,VyEnum);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+		
+		/* Get nodal functions value at gaussian point:*/
+		GetNodalFunctions(l1l2l3, gauss_l1l2l3);
+
+		/*Get nodal functions derivatives*/
+		GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
+
+		/*Get k derivative: dk/dx */
+		inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);
+
+		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
+		for (i=0;i<numgrids;i++){
+
+			//standard term dJ/dki
+			grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i];
+
+			//noise dampening d/dki(1/2*(dk/dx)^2)
+			grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
+			
+			//min dampening
+			if(drag<cm_mindmp_value){ 
+				grade_g_gaussian[i]+=cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
+			}
+
+			//max dampening
+			if(drag>cm_maxdmp_value){ 
+				grade_g_gaussian[i]+= - cm_maxdmp_slope*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);
+	delete friction;
+
+}
+/*}}}*/
+/*FUNCTION Tria::GradjDragStokes {{{1*/
+void  Tria::GradjDragStokes(Vec grad_g,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];
+	double       dh1dh3[NDOF2][numgrids];
+
+	/* grid data: */
 	double drag;
+	double alpha_complement;
+	Friction* friction=NULL;
 
 	/* gaussian points: */
@@ -3967,12 +4112,10 @@
 
 	/* parameters: */
+	double  vx,vy,vz;
+	double  lambda,mu,xi;
+	double  bed,thickness,Neff;
+	double  surface_normal[3];
+	double  bed_normal[3];
 	double  dk[NDOF2]; 
-	double  vx,vy;
-	double  lambda,mu;
-	double  bed,thickness,Neff;
-	
-	/*drag: */
-	double alpha_complement_list[numgrids];
-	double alpha_complement;
 
 	/*element vector at the gaussian points: */
@@ -3991,4 +4134,5 @@
 	/*inputs: */
 	bool shelf;
+	int  drag_type;
 
 	/*parameters: */
@@ -4001,4 +4145,5 @@
 	/*retrieve inputs :*/
 	inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+	inputs->GetParameterValue(&drag_type,DragTypeEnum);
 
 	/*retrieve some parameters: */
@@ -4009,5 +4154,4 @@
 	this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum);
 
-
 	/*Get out if shelf*/
 	if(shelf)return;
@@ -4017,13 +4161,7 @@
 	GetDofList1(&doflist1[0]);
 
-	/*Recover inputs: */
-	inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum);
-	inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum);
-	inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],3,ThicknessEnum);
-	inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],3,BedEnum);
-	inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],3,DragCoefficientEnum);
-	inputs->GetParameterValue(&drag_p,DragPEnum);
-	inputs->GetParameterValue(&drag_q,DragQEnum);
+	/*Build frictoin element, needed later: */
 	inputs->GetParameterValue(&drag_type,DragTypeEnum);
+	friction=new Friction("2d",inputs,matpar);
 
 	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
@@ -4038,243 +4176,7 @@
 		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
 
-		/*Build alpha_complement_list: */
-		if (drag_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=&dragcoefficient_list[0];
-			friction->bed=&bed_list[0];
-			friction->thickness=&thickness_list[0];
-			friction->vx=&vx_list[0];
-			friction->vy=&vy_list[0];
-			friction->p=drag_p;
-			friction->q=drag_q;
-
-			
-			if(friction->p!=1) ISSMERROR("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);
-		inputs->GetParameterValue(&drag, gauss_l1l2l3,DragCoefficientEnum);
-
-		/*recover lambda and mu: */
-		inputs->GetParameterValue(&lambda, gauss_l1l2l3,AdjointxEnum);
-		inputs->GetParameterValue(&mu, gauss_l1l2l3,AdjointyEnum);
-			
-		/*recover vx and vy: */
-		inputs->GetParameterValue(&vx, gauss_l1l2l3,VxEnum);
-		inputs->GetParameterValue(&vy, gauss_l1l2l3,VyEnum);
-
-		/* Get Jacobian determinant: */
-		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-		
-		/* Get nodal functions value at gaussian point:*/
-		GetNodalFunctions(l1l2l3, gauss_l1l2l3);
-
-		/*Get nodal functions derivatives*/
-		GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
-
-		/*Get k derivative: dk/dx */
-		inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);
-
-		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
-		for (i=0;i<numgrids;i++){
-
-			//standard term dJ/dki
-			grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i];
-
-			//noise dampening d/dki(1/2*(dk/dx)^2)
-			grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
-			
-			//min dampening
-			if(drag<cm_mindmp_value){ 
-				grade_g_gaussian[i]+=cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
-			}
-
-			//max dampening
-			if(drag>cm_maxdmp_value){ 
-				grade_g_gaussian[i]+= - cm_maxdmp_slope*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);
-
-}
-/*}}}*/
-/*FUNCTION Tria::GradjDragStokes {{{1*/
-void  Tria::GradjDragStokes(Vec grad_g,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];
-	double       dh1dh3[NDOF2][numgrids];
-
-	/* grid data: */
-	double vx_list[numgrids];
-	double vy_list[numgrids];
-	double vz_list[numgrids];
-	double drag;
-	double  thickness_list[numgrids];
-	double  bed_list[numgrids];
-	double  dragcoefficient_list[numgrids];
-	double  drag_p,drag_q;
-	double alpha_complement_list[numgrids];
-	double alpha_complement;
-
-	/* 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];
-	double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
-
-	/* parameters: */
-	double  vx,vy,vz;
-	double  lambda,mu,xi;
-	double  bed,thickness,Neff;
-	double  surface_normal[3];
-	double  bed_normal[3];
-	double  dk[NDOF2]; 
-
-	/*element vector at the gaussian points: */
-	double  grade_g[numgrids]={0.0};
-	double  grade_g_gaussian[numgrids];
-
-	/* Jacobian: */
-	double Jdet;
-
-	/*nodal functions: */
-	double l1l2l3[3];
-
-	/* strain rate: */
-	double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
-
-	/*inputs: */
-	bool shelf;
-	int  drag_type;
-
-	/*parameters: */
-	double  cm_noisedmp;
-	double  cm_mindmp_slope;
-	double  cm_mindmp_value;
-	double  cm_maxdmp_value;
-	double  cm_maxdmp_slope;
-
-	/*retrieve inputs :*/
-	inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
-	inputs->GetParameterValue(&drag_type,DragTypeEnum);
-
-	/*retrieve some parameters: */
-	this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
-	this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum);
-	this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum);
-	this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum);
-	this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum);
-
-	/*Get out if shelf*/
-	if(shelf)return;
-
-	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
-	GetDofList1(&doflist1[0]);
-
-	/*Build alpha_complement_list: */
-	if (drag_type==2){
-
-		/*Allocate friction object: */
-		Friction* friction=NewFriction();
-
-		inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum);
-		inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum);
-		inputs->GetParameterValues(&vz_list[0],&gaussgrids[0][0],3,VzAverageEnum);
-		inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],3,DragCoefficientEnum);
-		inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],3,BedEnum);
-		inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],3,ThicknessEnum);
-		inputs->GetParameterValue(&drag_p,DragPEnum);
-		inputs->GetParameterValue(&drag_q,DragQEnum);
-
-
-		/*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=&dragcoefficient_list[0];
-		friction->bed=&bed_list[0];
-		friction->thickness=&thickness_list[0];
-		friction->vx=&vx_list[0];
-		friction->vy=&vy_list[0];
-		friction->p=drag_p;
-		friction->q=drag_q;
-
-
-		if(friction->p!=1) ISSMERROR("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;
-	}
-
-	/* 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);
-
-	/* 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);
-
-		/*Recover alpha_complement and k: */
-		GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3);
+		/*Recover alpha_complement and drag: */
+		if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum);
+		else alpha_complement=0;
 		inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum);
 
@@ -4338,9 +4240,10 @@
 	VecSetValues(grad_g,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
 
-cleanup_and_return: 
+	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);
+	delete friction;
 
 }
Index: /issm/trunk/src/c/objects/Loads/Friction.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Friction.cpp	(revision 3832)
+++ /issm/trunk/src/c/objects/Loads/Friction.cpp	(revision 3833)
@@ -1,160 +1,116 @@
-/*!\file: Friction.cpp
- * \brief: wrapper for all friction parameters
- */ 
+/*!\file Friction.c
+ * \brief: implementation of the Friction object
+ */
 
+/*Headers:*/
+/*{{{1*/
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "stdio.h"
+#include <string.h>
 #include "../objects.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
 #include "../../shared/shared.h"
 #include "../../include/include.h"
+/*}}}*/	
 
-/*--------------------------------------------------
-	NewFriction
-  --------------------------------------------------*/
-
-/*FUNCTION NewFriction{{{1*/
-Friction* NewFriction(void)
-{
-	/* create a new Friction object */
-
-	return (Friction*)xmalloc(sizeof(Friction));
+/*methods: */
+/*FUNCTION Friction::Friction() {{{1*/
+Friction::Friction(){
+	this->element_type=NULL;
+	this->inputs=NULL;
+	this->matpar=NULL;
 }
 /*}}}*/
-/*FUNCTION FrictionEcho{{{1*/
-void  FrictionEcho(Friction* friction){
+/*FUNCTION Friction::Friction(char* element_type, Inputs* inputs,Matpar* matpar){{{1*/
+Friction::Friction(char* element_type_in,Inputs* inputs_in,Matpar* matpar_in){
 
-	int i;
-	
-	printf("Friction echo: \n");
-	printf("   element_type: %s\n",friction->element_type);
-	printf("   gravity: %g\n",friction->gravity);
-	printf("   rho_ice: %g\n",friction->rho_ice);
-	printf("   rho_water: %g\n",friction->rho_water);
-	printf("   p: %g\n",friction->p);
-	printf("   q: %g\n",friction->q);
-	printf("   analysis_type: %i\n",friction->analysis_type);
-
-	printf("K: ");
-	for(i=0;i<3;i++)printf("%g ",friction->K[i]);
-	printf("\n");
-
-	printf("bed: ");
-	for(i=0;i<3;i++)printf("%g ",friction->bed[i]);
-	printf("\n");
-
-	printf("thickness: ");
-	for(i=0;i<3;i++)printf("%g ",friction->thickness[i]);
-	printf("\n");
-
-	printf("vx: ");
-	for(i=0;i<3;i++)printf("%g ",friction->vx[i]);
-	printf("\n            ");
-	printf("vy: ");
-	for(i=0;i<3;i++)printf("%g ",friction->vy[i]);
-	printf("\n");
-	if(friction->vz){
-		printf("vz: ");
-		for(i=0;i<3;i++)printf("%g ",friction->vz[i]);
-		printf("\n");
-	}
+	this->inputs=inputs_in;
+	this->element_type=(char*)xmalloc((strlen(element_type_in)+1)*sizeof(char));
+	strcpy(this->element_type,element_type);
+	this->matpar=matpar_in;
 
 }
 /*}}}*/
-/*FUNCTION FrictionInit {{{1*/
-/*--------------------------------------------------
-	FrictionInit
-  --------------------------------------------------*/
-
-int FrictionInit(Friction* friction)
-{
-	friction->element_type=NULL;
-	friction->gravity=UNDEF;
-	friction->rho_ice=UNDEF;
-	friction->rho_water=UNDEF;
-	friction->K=NULL; 
-	friction->bed=NULL;
-	friction->thickness=NULL;
-	friction->vx=NULL;
-	friction->vy=NULL;
-	friction->vz=NULL;
-	friction->p=UNDEF;
-	friction->q=UNDEF;
-	friction->analysis_type=UNDEF;
-	
-	return 1;
+/*FUNCTION Friction::~Friction() {{{1*/
+Friction::~Friction(){
+	xfree((void**)&element_type);
 }
 /*}}}*/
-/*FUNCTION DeleteFriction {{{1*/
-/*--------------------------------------------------
-	DeleteFriction
-  --------------------------------------------------*/
-
-void DeleteFriction(Friction** pfriction)
-{
-	Friction* friction = *pfriction;
-
-	/*Just  erase element_type: */
-	xfree((void**)&friction->element_type);
-
-	/*Erase entire structure: */
-	xfree((void**)pfriction);
+/*FUNCTION Friction::Echo {{{1*/
+void Friction::Echo(void){
+	printf("Friction:\n");
+	printf("   element_type: %s\n",this->element_type);
+	inputs->Echo();
+	matpar->Echo();
 }
 /*}}}*/
-/*FUNCTION FrictionGetAlpha2 {{{1*/
-/*--------------------------------------------------
-	FrictionGetAlpha2
-  --------------------------------------------------*/
-void  FrictionGetAlpha2(double* alpha2, Friction* friction){
+/*FUNCTION Friction::GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum){{{1*/
+void Friction::GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum){
 
 	/*This routine calculates the basal friction coefficient 
-	alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p*
-	alpha2 is assumed to be double[3]: */
+	alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
 
 	/*diverse: */
-	int     i;
-	const int     numgrids=3;
-	double  Neff[numgrids];
 	double  r,s;
-	double  velocity_x[numgrids];
-	double  velocity_y[numgrids];
-	double  velocity_z[numgrids];
-	double  velocity_mag[numgrids];
+	double  drag_p, drag_q;
+	double  gravity,rho_ice,rho_water;
+	double  Neff;
+	double  thickness,bed;
+	double  vx,vy,vz,vmag;
+	double  drag_coefficient;
+	double  alpha2;
+
+
+	/*Recover parameters: */
+	inputs->GetParameterValue(&drag_p,DragPEnum);
+	inputs->GetParameterValue(&drag_q,DragQEnum);
+	inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);
+	inputs->GetParameterValue(&bed, gauss,BedEnum);
+	inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);
+
+	/*Get material parameters: */
+	gravity=matpar->GetG();
+	rho_ice=matpar->GetRhoIce();
+	rho_water=matpar->GetRhoWater();
 
 	//compute r and q coefficients: */
-	r=friction->q/friction->p;
-	s=1./friction->p;
+	r=drag_q/drag_p;
+	s=1./drag_p;
 		
 	//From bed and thickness, compute effective pressure when drag is viscous:
-	for(i=0;i<numgrids;i++){
+	Neff=gravity*(rho_ice*thickness+rho_water*bed);
+	
+	/*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
+	the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
+	for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
+	replace it by Neff=0 (ie, equival it to an ice shelf)*/
+	if (Neff<0)Neff=0;
 
-		Neff[i]=friction->gravity*(friction->rho_ice*friction->thickness[i]+friction->rho_water*friction->bed[i]);
+	if(strcmp(element_type,"2d")){
+		inputs->GetParameterValue(&vx, gauss,vxenum);
+		inputs->GetParameterValue(&vy, gauss,vyenum);
+		vmag=sqrt(pow(vx,2)+pow(vy,2));
+	}
+	else{
+		inputs->GetParameterValue(&vx, gauss,vxenum);
+		inputs->GetParameterValue(&vy, gauss,vyenum);
+		inputs->GetParameterValue(&vz, gauss,vzenum);
+		vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
+	}
+	
+	alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
 
-		/*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
-		the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
-		for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
-		replace it by Neff=0 (ie, equival it to an ice shelf)*/
-		if (Neff[i]<0)Neff[i]=0;
-
-		//We need the velocity magnitude to evaluate the basal stress:
-		if(strcmp(friction->element_type,"2d")){
-			velocity_x[i]=friction->vx[i];//velocities of size numgridsx2
-			velocity_y[i]=friction->vy[i];
-			velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2));
-		}
-		else{
-			velocity_x[i]=friction->vx[i]; //velocities of size numgridsx3
-			velocity_y[i]=friction->vy[i];
-			velocity_z[i]=friction->vz[i];
-			velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2)+pow(velocity_z[i],2));
-		}
-	
-		alpha2[i]=pow(friction->K[i],2)*pow(Neff[i],r)*pow(velocity_mag[i],(s-1));
-	}
+	/*Assign output pointers:*/
+	*palpha2=alpha2;
 }
 /*}}}*/
-/*FUNCTION FrictionGetAlphaComplement {{{1*/
-/*--------------------------------------------------
-	FrictionGetAlphaComplement
-  --------------------------------------------------*/
-void  FrictionGetAlphaComplement(double* alpha_complement, Friction* friction){
-
+/*FUNCTION Friction::GetAlphaComplement {{{1*/
+void Friction::GetAlphaComplement(double* palpha_complement, double* gauss,int vxenum,int vyenum){
+	
 	/* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. 
 	 * FrictionGetAlphaComplement is used in control methods on drag, and it computes: 
@@ -163,33 +119,51 @@
 	/*diverse: */
 	int     i;
-	const int     numgrids=3;
-	double  Neff[numgrids];
+	double  Neff;
 	double  r,s;
-	double  velocity_x[numgrids];
-	double  velocity_y[numgrids];
-	double  velocity_mag[numgrids];
+	double  vx;
+	double  vy;
+	double  vmag;
+	double  drag_p,drag_q;
+	double  drag_coefficient;
+	double  bed,thickness;
+	double  gravity,rho_ice,rho_water;
+	double  alpha_complement;
+
+	/*Recover parameters: */
+	inputs->GetParameterValue(&drag_p,DragPEnum);
+	inputs->GetParameterValue(&drag_q,DragQEnum);
+	inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);
+	inputs->GetParameterValue(&bed, gauss,BedEnum);
+	inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);
+
+	/*Get material parameters: */
+	gravity=matpar->GetG();
+	rho_ice=matpar->GetRhoIce();
+	rho_water=matpar->GetRhoWater();
+
 
 	//compute r and q coefficients: */
-	r=friction->q/friction->p;
-	s=1./friction->p;
+	r=drag_q/drag_p;
+	s=1./drag_p;
 		
 	//From bed and thickness, compute effective pressure when drag is viscous:
-	for(i=0;i<numgrids;i++){
+	Neff=gravity*(rho_ice*thickness+rho_water*bed);
 
-		Neff[i]=friction->gravity*(friction->rho_ice*friction->thickness[i]+friction->rho_water*friction->bed[i]);
+	/*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
+	  the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
+	  for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
+	  replace it by Neff=0 (ie, equival it to an ice shelf)*/
+	if (Neff<0)Neff=0;
 
-		/*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
-		the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
-		for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
-		replace it by Neff=0 (ie, equival it to an ice shelf)*/
-		if (Neff[i]<0)Neff[i]=0;
+	//We need the velocity magnitude to evaluate the basal stress:
+	inputs->GetParameterValue(&vx, gauss,vxenum);
+	inputs->GetParameterValue(&vy, gauss,vyenum);
+	vmag=sqrt(pow(vx,2)+pow(vy,2));
+	
+	alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
 
-		//We need the velocity magnitude to evaluate the basal stress:
-		velocity_x[i]=friction->vx[i]; //velocities of size numgridsx2
-		velocity_y[i]=friction->vy[i];
-		velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2));
-	
-		alpha_complement[i]=pow(Neff[i],r)*pow(velocity_mag[i],(s-1));
-	}
+	/*Assign output pointers:*/
+	*palpha_complement=alpha_complement;
+
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Loads/Friction.h
===================================================================
--- /issm/trunk/src/c/objects/Loads/Friction.h	(revision 3832)
+++ /issm/trunk/src/c/objects/Loads/Friction.h	(revision 3833)
@@ -1,38 +1,33 @@
-/*! \file Friction.h
- *  \brief Friction object, interface, declaration
+/*!\file Friction.h
+ * \brief: header file for friction object
  */
 
-#ifndef _FRICTION_H
-#define _FRICTION_H
+#ifndef _FRICTION_H_
+#define _FRICTION_H_
 
-/*!Friction declaration: */
-struct Friction {
+/*Headers:*/
+/*{{{1*/
+class Inputs;
+class Matpar;
+/*}}}*/
 
-	char*   element_type;
-	double  gravity;
-	double  rho_ice;
-	double  rho_water;
-	double* K;
-	double* bed;
-	double* thickness;
-	double* vx;
-	double* vy;
-	double* vz;
-	double  p;
-	double  q;
-	int     analysis_type;
+class Friction{
+
+	public:
+
+		char* element_type;
+		Inputs* inputs;
+		Matpar* matpar;
+
+		/*methods: */
+		Friction();
+		Friction(char* element_type, Inputs* inputs,Matpar* matpar);
+		~Friction();
+	
+		void  Echo(void);
+		void  GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum);
+		void  GetAlphaComplement(double* alpha_complement, double* gauss,int vxenum,int vyenum);
 
 };
 
-
-/* creation, initialisation: */
-
-	Friction*	NewFriction(void);
-	int	FrictionInit(Friction* friction);
-	void  DeleteFriction(Friction** pfriction);
-	void  FrictionGetAlpha2(double* alpha2, Friction* friction);
-	void  FrictionGetAlphaComplement(double* alpha_complement, Friction* friction);
-	void  FrictionEcho(Friction* friction);
-
-#endif  /* _FRICTION_H */
-
+#endif  /* _FRICTION_H_ */
Index: sm/trunk/src/c/objects/Loads/Friction2.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Friction2.cpp	(revision 3832)
+++ 	(revision )
@@ -1,169 +1,0 @@
-/*!\file Friction2.c
- * \brief: implementation of the Friction2 object
- */
-
-/*Headers:*/
-/*{{{1*/
-#ifdef HAVE_CONFIG_H
-	#include "config.h"
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#include "stdio.h"
-#include <string.h>
-#include "../objects.h"
-#include "../../EnumDefinitions/EnumDefinitions.h"
-#include "../../shared/shared.h"
-#include "../../include/include.h"
-/*}}}*/	
-
-/*methods: */
-/*FUNCTION Friction2::Friction2() {{{1*/
-Friction2::Friction2(){
-	this->element_type=NULL;
-	this->inputs=NULL;
-	this->matpar=NULL;
-}
-/*}}}*/
-/*FUNCTION Friction2::Friction2(char* element_type, Inputs* inputs,Matpar* matpar){{{1*/
-Friction2::Friction2(char* element_type_in,Inputs* inputs_in,Matpar* matpar_in){
-
-	this->inputs=inputs_in;
-	this->element_type=(char*)xmalloc((strlen(element_type_in)+1)*sizeof(char));
-	strcpy(this->element_type,element_type);
-	this->matpar=matpar_in;
-
-}
-/*}}}*/
-/*FUNCTION Friction2::~Friction2() {{{1*/
-Friction2::~Friction2(){
-	xfree((void**)&element_type);
-}
-/*}}}*/
-/*FUNCTION Friction2::Echo {{{1*/
-void Friction2::Echo(void){
-	printf("Friction2:\n");
-	printf("   element_type: %s\n",this->element_type);
-	inputs->Echo();
-	matpar->Echo();
-}
-/*}}}*/
-/*FUNCTION Friction2::GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum){{{1*/
-void Friction2::GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum){
-
-	/*This routine calculates the basal friction coefficient 
-	alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
-
-	/*diverse: */
-	double  r,s;
-	double  drag_p, drag_q;
-	double  gravity,rho_ice,rho_water;
-	double  Neff;
-	double  thickness,bed;
-	double  vx,vy,vz,vmag;
-	double  drag_coefficient;
-	double  alpha2;
-
-
-	/*Recover parameters: */
-	inputs->GetParameterValue(&drag_p,DragPEnum);
-	inputs->GetParameterValue(&drag_q,DragQEnum);
-	inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);
-	inputs->GetParameterValue(&bed, gauss,BedEnum);
-	inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);
-
-	/*Get material parameters: */
-	gravity=matpar->GetG();
-	rho_ice=matpar->GetRhoIce();
-	rho_water=matpar->GetRhoWater();
-
-	//compute r and q coefficients: */
-	r=drag_q/drag_p;
-	s=1./drag_p;
-		
-	//From bed and thickness, compute effective pressure when drag is viscous:
-	Neff=gravity*(rho_ice*thickness+rho_water*bed);
-	
-	/*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
-	the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
-	for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
-	replace it by Neff=0 (ie, equival it to an ice shelf)*/
-	if (Neff<0)Neff=0;
-
-	if(strcmp(element_type,"2d")){
-		inputs->GetParameterValue(&vx, gauss,vxenum);
-		inputs->GetParameterValue(&vy, gauss,vyenum);
-		vmag=sqrt(pow(vx,2)+pow(vy,2));
-	}
-	else{
-		inputs->GetParameterValue(&vx, gauss,vxenum);
-		inputs->GetParameterValue(&vy, gauss,vyenum);
-		inputs->GetParameterValue(&vz, gauss,vzenum);
-		vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
-	}
-	
-	alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
-
-	/*Assign output pointers:*/
-	*palpha2=alpha2;
-}
-/*}}}*/
-/*FUNCTION Friction2::GetAlphaComplement {{{1*/
-void Friction2::GetAlphaComplement(double* palpha_complement, double* gauss,int vxenum,int vyenum){
-	
-	/* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. 
-	 * FrictionGetAlphaComplement is used in control methods on drag, and it computes: 
-	 * alpha_complement= Neff ^r * vel ^s*/
-
-	/*diverse: */
-	int     i;
-	double  Neff;
-	double  r,s;
-	double  vx;
-	double  vy;
-	double  vmag;
-	double  drag_p,drag_q;
-	double  drag_coefficient;
-	double  bed,thickness;
-	double  gravity,rho_ice,rho_water;
-	double  alpha_complement;
-
-	/*Recover parameters: */
-	inputs->GetParameterValue(&drag_p,DragPEnum);
-	inputs->GetParameterValue(&drag_q,DragQEnum);
-	inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);
-	inputs->GetParameterValue(&bed, gauss,BedEnum);
-	inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);
-
-	/*Get material parameters: */
-	gravity=matpar->GetG();
-	rho_ice=matpar->GetRhoIce();
-	rho_water=matpar->GetRhoWater();
-
-
-	//compute r and q coefficients: */
-	r=drag_q/drag_p;
-	s=1./drag_p;
-		
-	//From bed and thickness, compute effective pressure when drag is viscous:
-	Neff=gravity*(rho_ice*thickness+rho_water*bed);
-
-	/*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
-	  the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
-	  for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
-	  replace it by Neff=0 (ie, equival it to an ice shelf)*/
-	if (Neff<0)Neff=0;
-
-	//We need the velocity magnitude to evaluate the basal stress:
-	inputs->GetParameterValue(&vx, gauss,vxenum);
-	inputs->GetParameterValue(&vy, gauss,vyenum);
-	vmag=sqrt(pow(vx,2)+pow(vy,2));
-	
-	alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
-
-	/*Assign output pointers:*/
-	*palpha_complement=alpha_complement;
-
-}
-/*}}}*/
Index: sm/trunk/src/c/objects/Loads/Friction2.h
===================================================================
--- /issm/trunk/src/c/objects/Loads/Friction2.h	(revision 3832)
+++ 	(revision )
@@ -1,33 +1,0 @@
-/*!\file Friction2.h
- * \brief: header file for friction2 object
- */
-
-#ifndef _FRICTION2_H_
-#define _FRICTION2_H_
-
-/*Headers:*/
-/*{{{1*/
-class Inputs;
-class Matpar;
-/*}}}*/
-
-class Friction2{
-
-	public:
-
-		char* element_type;
-		Inputs* inputs;
-		Matpar* matpar;
-
-		/*methods: */
-		Friction2();
-		Friction2(char* element_type, Inputs* inputs,Matpar* matpar);
-		~Friction2();
-	
-		void  Echo(void);
-		void  GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum);
-		void  GetAlphaComplement(double* alpha_complement, double* gauss,int vxenum,int vyenum);
-
-};
-
-#endif  /* _FRICTION2_H_ */
Index: /issm/trunk/src/c/objects/objects.h
===================================================================
--- /issm/trunk/src/c/objects/objects.h	(revision 3832)
+++ /issm/trunk/src/c/objects/objects.h	(revision 3833)
@@ -26,5 +26,4 @@
 /*Loads: */
 #include "./Loads/Friction.h"
-#include "./Loads/Friction2.h"
 #include "./Loads/Icefront.h"
 #include "./Loads/Numericalflux.h"