Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 97)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 98)
@@ -307,21 +307,13 @@
 	double B[3][numdofs];
 	double Bprime[3][numdofs];
-	double L[2][numdofs];
 	double D[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}};              // material matrix, simple scalar matrix.
 	double D_scalar;
-	double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
-	double DL_scalar;
 
 	/* local element matrices: */
 	double Ke_gg[numdofs][numdofs]; //local element stiffness matrix 
 	double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix evaluated at the gaussian point.
-	double Ke_gg_drag_gaussian[numdofs][numdofs]; //stiffness matrix contribution from drag
 	
 	double Jdet;
 	
-	/*slope: */
-	double  slope[2]={0.0,0.0};
-	double  slope_magnitude;
-
 	/*input parameters for structural analysis (diagnostic): */
 	double* velocity_param=NULL;
@@ -330,12 +322,4 @@
 	double  oldvxvy_list[numgrids][2];
 	double  thickness;
-
-	/*friction: */
-	double alpha2_list[numgrids]={0.0,0.0,0.0};
-	double alpha2;
-
-	double MAXSLOPE=.06; // 6 %
-	double MOUNTAINKEXPONENT=12;
-
 
 	/*recover extra inputs from users, at current convergence iteration: */
@@ -390,37 +374,4 @@
 
 
-	/*Build alpha2_list used by drag stiffness matrix*/
-	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=&vxvy_list[0][0];
-		friction->p=p;
-		friction->q=q;
-
-		/*Compute alpha2_list: */
-		FrictionGetAlpha2(&alpha2_list[0],friction);
-
-		/*Erase friction object: */
-		DeleteFriction(&friction);
-	}
-
-	#ifdef _DEBUGELEMENTS_
-	if(my_rank==RANK && id==ELID){ 
-		printf("   alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]);
-	}
-	#endif
-
 	/* 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, 2);
@@ -446,17 +397,4 @@
 		/*Compute thickness at gaussian point: */
 		GetParameterValue(&thickness, &h[0],gauss_l1l2l3);
-
-		// If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
-		//velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
-		if(!shelf){
-			GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3);
-			slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
-
-			if (slope_magnitude>MAXSLOPE){
-				alpha2_list[0]=pow(10,MOUNTAINKEXPONENT);
-				alpha2_list[1]=pow(10,MOUNTAINKEXPONENT);
-				alpha2_list[2]=pow(10,MOUNTAINKEXPONENT);
-			}
-		}
 
 		/*Get strain rate from velocity: */
@@ -501,9 +439,4 @@
 		GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
 
-		/*Get L matrix if viscous basal drag present: */
-		if((friction_type==2) && (!shelf)){
-			GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
-		}	
-
 		/*  Do the triple product tB*D*Bprime: */
 		TripleMultiply( &B[0][0],3,numdofs,1,
@@ -515,33 +448,6 @@
 		for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
 		
-		/*Now, take care of the basal friction if there is any: */
-		if((!shelf) && (friction_type==2)){
-
-			GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3);
-
-			if (velocity_param){
-				DL_scalar=alpha2*gauss_weight*Jdet;
-			}
-			else{
-				DL_scalar=0;
-			}
-				
-			for (i=0;i<2;i++){
-				DL[i][i]=DL_scalar;
-			}
-			
-			/*  Do the triple producte tL*D*L: */
-			TripleMultiply( &L[0][0],2,numdofs,1,
-						&DL[0][0],2,2,0,
-						&L[0][0],2,numdofs,0,
-						&Ke_gg_drag_gaussian[0][0],0);
-
-			for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_drag_gaussian[i][j];
-		}
-
-
 		#ifdef _DEBUGELEMENTS_
 		if(my_rank==RANK && id==ELID){ 
-			printf("      alpha2 %g\n",alpha2);
 			printf("      B:\n");
 			for(i=0;i<3;i++){
@@ -558,12 +464,4 @@
 				printf("\n");
 			}
-			printf("      L:\n");
-			for(i=0;i<2;i++){
-				for(j=0;j<numdofs;j++){
-					printf("%g ",L[i][j]);
-				}
-				printf("\n");
-			}
-
 		}
 		#endif
@@ -572,4 +470,11 @@
 	/*Add Ke_gg to global matrix Kgg: */
 	MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
+
+
+	/*Do not forget to include friction: */
+	if(!shelf){
+		CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type);
+	}
+
 
 	#ifdef _DEBUGELEMENTS_
@@ -597,5 +502,185 @@
 
 }
-	
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::CreateKMatrixDiagnosticHorizFriction"
+void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,ParameterInputs* inputs,int analysis_type){
+
+
+	/* local declarations */
+	int             i,j;
+
+	/* node data: */
+	const int    numgrids=3;
+	const int    numdofs=2*numgrids;
+	double       xyz_list[numgrids][3];
+	int          doflist[numdofs];
+	int          numberofdofspernode;
+	
+	/* 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];
+
+	/* matrices: */
+	double L[2][numdofs];
+	double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
+	double DL_scalar;
+
+	/* local element matrices: */
+	double Ke_gg[numdofs][numdofs]; //local element stiffness matrix 
+	double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix contribution from drag
+	
+	double Jdet;
+	
+	/*slope: */
+	double  slope[2]={0.0,0.0};
+	double  slope_magnitude;
+
+	/*input parameters for structural analysis (diagnostic): */
+	double* velocity_param=NULL;
+	double  vxvy_list[numgrids][2];
+
+	/*friction: */
+	double alpha2_list[numgrids]={0.0,0.0,0.0};
+	double alpha2;
+
+	double MAXSLOPE=.06; // 6 %
+	double MOUNTAINKEXPONENT=10;
+
+	/*recover extra inputs from users, at current convergence iteration: */
+	velocity_param=ParameterInputsRecover(inputs,"velocity"); 
+
+	/* Get node coordinates and dof list: */
+	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	/* Set Ke_gg to 0: */
+	for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
+
+	if (shelf){
+		/*no friction, do nothing*/
+		return;
+	}
+
+	if (friction_type!=2)throw ErrorException(__FUNCT__," non-viscous friction not supported yet!");
+
+	/*Initialize vxvy_list: */
+	for(i=0;i<numgrids;i++){
+		if(velocity_param){
+			vxvy_list[i][0]=velocity_param[doflist[i*numberofdofspernode+0]];
+			vxvy_list[i][1]=velocity_param[doflist[i*numberofdofspernode+1]];
+		}
+		else{
+			vxvy_list[i][0]=0;
+			vxvy_list[i][1]=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");
+	
+	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=&vxvy_list[0][0];
+	friction->p=p;
+	friction->q=q;
+
+	/*Compute alpha2_list: */
+	FrictionGetAlpha2(&alpha2_list[0],friction);
+
+	/*Erase friction object: */
+	DeleteFriction(&friction);
+
+	#ifdef _DEBUGELEMENTS_
+	if(my_rank==RANK && id==ELID){ 
+		printf("   alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]);
+	}
+	#endif
+
+	/* 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, 2);
+
+	#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);
+
+
+		// If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
+		//velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
+		GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3);
+		slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
+
+		if (slope_magnitude>MAXSLOPE){
+			alpha2_list[0]=pow(10,MOUNTAINKEXPONENT);
+			alpha2_list[1]=pow(10,MOUNTAINKEXPONENT);
+			alpha2_list[2]=pow(10,MOUNTAINKEXPONENT);
+		}
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+
+		/*Get L matrix: */
+		GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
+
+		/*Now, take care of the basal friction if there is any: */
+		GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3);
+
+		if (velocity_param){
+			DL_scalar=alpha2*gauss_weight*Jdet;
+		}
+		else{
+			DL_scalar=0;
+		}
+			
+		for (i=0;i<2;i++){
+			DL[i][i]=DL_scalar;
+		}
+		
+		/*  Do the triple producte tL*D*L: */
+		TripleMultiply( &L[0][0],2,numdofs,1,
+					&DL[0][0],2,2,0,
+					&L[0][0],2,numdofs,0,
+					&Ke_gg_gaussian[0][0],0);
+
+		for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+
+	} // for (ig=0; ig<num_gauss; ig++)
+
+	/*Add Ke_gg to global matrix Kgg: */
+	MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,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__ 
@@ -1928,2 +2013,31 @@
 	return Jelem;
 }
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::NodeConfiguration"
+void  Tria::NodeConfiguration(int* tria_node_ids,Node* tria_nodes[3],int* tria_node_offsets){
+
+	int i;
+	for(i=0;i<3;i++){
+		node_ids[i]=tria_node_ids[i];
+		nodes[i]=tria_nodes[i];
+		node_offsets[i]=tria_node_offsets[i];
+	}
+
+}
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::MaticeConfiguration"
+void  Tria::MaticeConfiguration(Matice* tria_matice,int tria_matice_offset){
+	matice=tria_matice;
+	matice_offset=tria_matice_offset;
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::MaticeConfiguration"
+void  Tria::MatparConfiguration(Matpar* tria_matpar,int tria_matpar_offset){
+
+	matpar=tria_matpar;
+	matpar_offset=tria_matpar_offset;
+
+}
+