Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 8116)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 8117)
@@ -436,6 +436,9 @@
 	RheologyLawEnum,
 	PatersonEnum,
-	ArrheniusEnum
-	/*}}}*/
+	ArrheniusEnum,
+	/*}}}*/
+	VxMeshEnum,
+	VyMeshEnum,
+	VzMeshEnum
 };
 
Index: /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp	(revision 8116)
+++ /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp	(revision 8117)
@@ -384,4 +384,7 @@
 		case PatersonEnum : return "Paterson";
 		case ArrheniusEnum : return "Arrhenius";
+		case VxMeshEnum : return "VxMesh";
+		case VyMeshEnum : return "VyMesh";
+		case VzMeshEnum : return "VzMesh";
 		default : return "unknown";
 
Index: /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp	(revision 8116)
+++ /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp	(revision 8117)
@@ -382,4 +382,7 @@
 	else if (strcmp(name,"Paterson")==0) return PatersonEnum;
 	else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
+	else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
+	else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
+	else if (strcmp(name,"VzMesh")==0) return VzMeshEnum;
 	else _error_("Enum %s not found",name);
 
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 8116)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 8117)
@@ -1703,5 +1703,5 @@
 	int        artdiff;
 	int        i,j,ig,found=0;
-	double     Jdet,u,v,w,epsvel;
+	double     Jdet,u,v,w,um,vm,wm,epsvel;
 	double     gravity,rho_ice,rho_water;
 	double     heatcapacity,thermalconductivity,dt;
@@ -1740,7 +1740,10 @@
 	this->parameters->FindParam(&artdiff,ArtDiffEnum);
 	this->parameters->FindParam(&epsvel,EpsVelEnum);
-	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
-	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
-	Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
+	Input* vx_input=inputs->GetInput(VxEnum);      _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);      _assert_(vy_input);
+	Input* vz_input=inputs->GetInput(VzEnum);      _assert_(vz_input);
+	Input* vxm_input=inputs->GetInput(VxMeshEnum); _assert_(vxm_input);
+	Input* vym_input=inputs->GetInput(VyMeshEnum); _assert_(vym_input);
+	Input* vzm_input=inputs->GetInput(VzMeshEnum); _assert_(vzm_input);
 	if (artdiff==2) diameter=MinEdgeLength(xyz_list);
 
@@ -1777,11 +1780,14 @@
 		vy_input->GetParameterValue(&v, gauss);
 		vz_input->GetParameterValue(&w, gauss);
+		vxm_input->GetParameterValue(&um,gauss);
+		vym_input->GetParameterValue(&vm,gauss);
+		vzm_input->GetParameterValue(&wm,gauss);
 
 		D_scalar_advec=gauss->weight*Jdet;
 		if(dt) D_scalar_advec=D_scalar_advec*dt;
 
-		D[0][0]=D_scalar_advec*u;D[0][1]=0;         D[0][2]=0;
-		D[1][0]=0;         D[1][1]=D_scalar_advec*v;D[1][2]=0;
-		D[2][0]=0;         D[2][1]=0;         D[2][2]=D_scalar_advec*w;
+		D[0][0]=D_scalar_advec*(um-u);D[0][1]=0;                    D[0][2]=0;
+		D[1][0]=0;                    D[1][1]=D_scalar_advec*(v-vm);D[1][2]=0;
+		D[2][0]=0;                    D[2][1]=0;                    D[2][2]=D_scalar_advec*(w-wm);
 
 		TripleMultiply(&B_advec[0][0],3,numdof,1,
@@ -1810,5 +1816,5 @@
 		if(artdiff==1){
 			/*Build K: */
-			D_scalar_artdiff=gauss->weight*Jdet/(pow(u,2)+pow(v,2)+epsvel);
+			D_scalar_artdiff=gauss->weight*Jdet/(pow(u-um,2)+pow(v-vm,2)+epsvel);
 			if(dt) D_scalar_artdiff=D_scalar_artdiff*dt;
 			K[0][0]=D_scalar_artdiff*pow(u,2);       K[0][1]=D_scalar_artdiff*fabs(u)*fabs(v);
@@ -1826,9 +1832,10 @@
 			GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
 
-			tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
+			tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
 
 			for(i=0;i<numdof;i++){
 				for(j=0;j<numdof;j++){
-					Ke_gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i])*(u*dh1dh6[0][j]+v*dh1dh6[1][j]+w*dh1dh6[2][j]);
+					Ke_gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*
+					  ((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i])*((u-um)*dh1dh6[0][j]+(v-vm)*dh1dh6[1][j]+(w-wm)*dh1dh6[2][j]);
 				}
 			}
@@ -1836,5 +1843,5 @@
 				for(i=0;i<numdof;i++){
 					for(j=0;j<numdof;j++){
-						Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
+						Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i]);
 					}
 				}
@@ -6238,4 +6245,12 @@
 			break;
 
+		case ThermalAnalysisEnum:
+			/*Initialize mesh velocity*/
+			for(i=0;i<6;i++)nodeinputs[i]=0;
+			this->inputs->AddInput(new PentaVertexInput(VxMeshEnum,nodeinputs));
+			this->inputs->AddInput(new PentaVertexInput(VyMeshEnum,nodeinputs));
+			this->inputs->AddInput(new PentaVertexInput(VzMeshEnum,nodeinputs));
+			break;
+
 		default:
 			/*No update for other solution types*/