Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4732)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4733)
@@ -937,4 +937,7 @@
 	else if(analysis_type==DiagnosticStokesAnalysisEnum){
 		GetSolutionFromInputsDiagnosticStokes(solution);
+	}
+	else if(analysis_type==ThermalAnalysisEnum){
+		GetSolutionFromInputsThermal(solution);
 	}
 	else{
@@ -2874,4 +2877,36 @@
 	VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
 
+}
+/*}}}*/
+/*FUNCTION Penta::GetSolutionFromInputsThermal{{{1*/
+void  Penta::GetSolutionFromInputsThermal(Vec solution){
+
+	int i;
+
+	const int    numvertices=6;
+	const int    numdofpervertex=1;
+	const int    numdof=numdofpervertex*numvertices;
+	double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+
+	int          doflist[numdof];
+	double       values[numdof];
+	double       vz;
+
+	int          dummy;
+
+	/*Get dof list: */
+	GetDofList(&doflist[0],&dummy);
+
+	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+	/*P1 element only for now*/
+	for(i=0;i<numvertices;i++){
+
+		/*Recover vz */
+		inputs->GetParameterValue(&vz,&gauss[i][0],TemperatureEnum);
+		values[i]=vz;
+	}
+
+	/*Add value to global vector*/
+	VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4732)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4733)
@@ -167,4 +167,5 @@
 		void	  GetSolutionFromInputsDiagnosticStokes(Vec solutiong);
 		void	  GetSolutionFromInputsDiagnosticVert(Vec solutiong);
+		void	  GetSolutionFromInputsThermal(Vec solutiong);
 		void	  GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);
 		void	  GetStrainRateStokes(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4732)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4733)
@@ -4665,4 +4665,7 @@
 	double    alpha2_list[numgrids];                                       //TO BE DELETED
 	double    gauss[numgrids][numgrids] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
+	double vx_list[numgrids]; //TO BE DELETED
+	double vy_list[numgrids]; //TO BE DELETED
+	double basalfriction_list[numgrids]; //TO BE DELETED
 
 	/* gaussian points: */
@@ -4703,7 +4706,12 @@
 	friction=new Friction("3d",inputs,matpar,analysis_type);
 
-	/*COMPUT alpha2_list (TO BE DELETED)*/
+	/*COMPUT alpha2_list and basalfriction_list (TO BE DELETED)*/
 	for(i=0;i<numgrids;i++){
-		friction->GetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum);
+		friction->GetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum); //TO BE DELETED
+	}
+	inputs->GetParameterValues(&vx_list[0],&gauss[0][0],3,VxEnum); //TO BE DELETED
+	inputs->GetParameterValues(&vy_list[0],&gauss[0][0],3,VyEnum); //TO BE DELETED
+	for(i=0;i<numgrids;i++){
+		basalfriction_list[i]=alpha2_list[i]*(pow(vx_list[i],(double)2.0)+pow(vy_list[i],(double)2.0)); //TO BE DELETED
 	}
 	
@@ -4729,9 +4737,5 @@
 		/*Friction: */
 		//friction->GetAlpha2(&alpha2,&gauss_coord[0],VxEnum,VyEnum,VzEnum);
-		GetParameterValue(&alpha2,&alpha2_list[0],gauss_coord); // TO BE DELETED
-
-		inputs->GetParameterValue(&vx, &gauss_coord[0],VxEnum);
-		inputs->GetParameterValue(&vy, &gauss_coord[0],VyEnum);
-		basalfriction= alpha2*(pow(vx,(double)2.0)+pow(vy,(double)2.0));
+		GetParameterValue(&basalfriction,&basalfriction_list[0],gauss_coord); // TO BE DELETED
 		
 		/*Calculate scalar parameter*/