Index: /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp	(revision 1675)
+++ /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp	(revision 1676)
@@ -198,4 +198,5 @@
 		penta_epsvel=model->epsvel;
 		penta_onwater=(bool)*(model->elementonwater+i);
+		penta_artdiff=model->artificial_diffusivity;
 
 		/*We need the field collapse for transient, so that we can use compute B with the average temperature*/
Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 1675)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 1676)
@@ -3389,4 +3389,5 @@
 	double  gauss_weight,area_gauss_weight,vert_gauss_weight;
 	double  gauss_coord[4];
+	double  gauss_l1l2l3[3];
 
 	int area_order=5;
@@ -3395,4 +3396,5 @@
 	int     dofs[3]={0,1,2};
 	double  dt;
+	double  K[2][2]={0.0};
 
 	double  vxvyvz_list[numgrids][3];
@@ -3411,4 +3413,5 @@
 	double     Ke_gaussian_conduct[numdof][numdof];
 	double     Ke_gaussian_advec[numdof][numdof];
+	double     Ke_gaussian_artdiff[numdof][numdof];
 	double     Ke_gaussian_transient[numdof][numdof];
 	double     B[3][numdof];
@@ -3416,4 +3419,5 @@
 	double     B_conduct[3][numdof];
 	double     B_advec[3][numdof];
+	double     B_artdiff[2][numdof];
 	double     Bprime_advec[3][numdof];
 	double     L[numdof];
@@ -3425,4 +3429,5 @@
 	double     tBD_conduct[3][numdof];
 	double     tBD_advec[3][numdof];
+	double     tBD_artdiff[3][numdof];
 	double     tLD[numdof];
 
@@ -3466,4 +3471,5 @@
 		vz_list[i]=vxvyvz_list[i][2];
 	}
+
 
 	/* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
@@ -3556,9 +3562,34 @@
 				}
 			}
+
+			/*Artifficial diffusivity*/
+			if(artdiff){
+				/*Build K: */
+				D_scalar=gauss_weight*Jdet/(pow(u,2)+pow(v,2)+epsvel);
+				if(sub_analysis_type!=SteadyAnalysisEnum()){
+					D_scalar=D_scalar*dt;
+				}
+				K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*fabs(u)*fabs(v);
+				K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2);
+
+				/*Get B_artdiff: */
+				GetB_artdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss_coord); 
+
+				/*  Do the triple product B'*K*B: */
+				MatrixMultiply(&B_artdiff[0][0],2,numdof,1,&K[0][0],2,2,0,&tBD_artdiff[0][0],0);
+				MatrixMultiply(&tBD_artdiff[0][0],numdof,2,0,&B_artdiff[0][0],2,numdof,0,&Ke_gaussian_artdiff[0][0],0);
+			}
+			else{
+				for(i=0;i<numdof;i++){
+					for(j=0;j<numdof;j++){
+						Ke_gaussian_artdiff[i][j]=0;
+					}
+				}
+			}
 					
 			/*Add Ke_gaussian to pKe: */
 			for(i=0;i<numdof;i++){
 				for(j=0;j<numdof;j++){
-					K_terms[i][j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j];
+					K_terms[i][j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j];
 				}
 			}
@@ -3622,11 +3653,43 @@
 
 #undef __FUNCT__ 
-#define __FUNCT__ "Penta::GetB_advec"
-void Penta::GetB_advec(double* B_advec, double* xyz_list, double* gauss_coord){
+#define __FUNCT__ "Penta::GetB_artdiff"
+void Penta::GetB_artdiff(double* B_artdiff, double* xyz_list, double* gauss_coord){
 
 	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 
 	 * For grid i, Bi' can be expressed in the basic coordinate system
 	 * by: 
-	 *       Bi_conduct_basic=[ h ]
+	 *       Bi_artdiff_basic=[ dh/dx ]
+	 *                       [ dh/dy ]
+	 * where h is the interpolation function for grid i.
+	 *
+	 * We assume B has been allocated already, of size: 2x(DOFPERGRID*numgrids)
+	 */
+
+	int i;
+	const int calculationdof=3;
+	const int numgrids=6;
+	int DOFPERGRID=1;
+
+	/*Same thing in the basic coordinate system: */
+	double dh1dh6_basic[calculationdof][numgrids];
+
+	/*Get dh1dh2dh3 in basic coordinates system : */
+	GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
+
+	/*Build B': */
+	for (i=0;i<numgrids;i++){
+		*(B_artdiff+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6_basic[0][i]; 
+		*(B_artdiff+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6_basic[1][i]; 
+	}
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta::GetB_advec"
+void Penta::GetB_advec(double* B_advec, double* xyz_list, double* gauss_coord){
+
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 
+	 * For grid i, Bi' can be expressed in the basic coordinate system
+	 * by: 
+	 *       Bi_advec_basic =[ h ]
 	 *                       [ h ]
 	 *                       [ h ]
Index: /issm/trunk/src/c/objects/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Penta.h	(revision 1675)
+++ /issm/trunk/src/c/objects/Penta.h	(revision 1676)
@@ -137,4 +137,5 @@
 		void  GetB_advec(double* B_advec, double* xyz_list, double* gauss_coord);
 		void  GetBprime_advec(double* Bprime_advec, double* xyz_list, double* gauss_coord);
+		void  GetB_artdiff(double* B_artdiff, double* xyz_list, double* gauss_coord);
 		void  CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
 		void  CreatePVectorThermal( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type);
Index: /issm/trunk/src/m/classes/public/solveparallel.m
===================================================================
--- /issm/trunk/src/m/classes/public/solveparallel.m	(revision 1675)
+++ /issm/trunk/src/m/classes/public/solveparallel.m	(revision 1676)
@@ -30,4 +30,5 @@
 	waitonlock([executionpath '/' md.name '.lock']);
 	%load results
+	displaystring(md.debug,'loading results from cluster');
 	md=loadresultsfromcluster(md);
 end