Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 6211)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 6212)
@@ -2858,9 +2858,11 @@
 
 	/*Intermediaries */
-	bool       artdiff;
+	int        artdiff;
 	int        i,j,ig,found=0;
 	double     Jdet,u,v,w,epsvel;
 	double     gravity,rho_ice,rho_water;
 	double     heatcapacity,thermalconductivity,dt;
+	double     scalar_artdiff;
+	double     tau_parameter,diameter,normu;
 	double     xyz_list[NUMVERTICES][3];
 	double     B[3][numdof];
@@ -2871,4 +2873,5 @@
 	double     Bprime_advec[3][numdof];
 	double     L[numdof];
+	double     dh1dh6[3][6];
 	double     D_scalar;
 	double     D[3][3];
@@ -2903,4 +2906,5 @@
 	Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
 	Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input);
+	if (artdiff==2) diameter=MinEdgeLength(xyz_list);
 
 	/* Start  looping on the number of gaussian points: */
@@ -2961,5 +2965,5 @@
 		/*Artifficial diffusivity*/
 
-		if(artdiff){
+		if(artdiff==1){
 			/*Build K: */
 			D_scalar=gauss->weight*Jdet/(pow(u,2)+pow(v,2)+epsvel);
@@ -2972,4 +2976,21 @@
 			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 if(artdiff==2){
+
+			GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
+
+			normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
+			if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){
+				tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity));
+			}
+			else tau_parameter=diameter/(2*normu);
+			scalar_artdiff=tau_parameter*gauss->weight*Jdet;
+
+			for(i=0;i<numdof;i++){
+				for(j=0;j<numdof;j++){
+					Ke_gaussian_artdiff[i][j]=scalar_artdiff*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i])*(u*dh1dh6[0][j]+v*dh1dh6[1][j]+w*dh1dh6[2][j]);
+				}
+			}
 		}
 		else{
@@ -3746,12 +3767,16 @@
 	/*Intermediaries*/
 	int    i,j,ig,found=0;
-	int    friction_type;
+	int    friction_type,artdiff;
 	double Jdet,phi,dt;
 	double rho_ice,heatcapacity;
+	double thermalconductivity;
 	double viscosity,temperature;
+	double tau_parameter,diameter,normu;
+	double u,v,w,scalar_artdiff;
 	double scalar_def,scalar_transient;
 	double temperature_list[NUMVERTICES];
 	double xyz_list[NUMVERTICES][3];
 	double L[numdof];
+	double dh1dh6[3][6];
 	double epsilon[6];
 	GaussPenta *gauss=NULL;
@@ -3763,7 +3788,9 @@
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-	this->parameters->FindParam(&dt,DtEnum);
 	rho_ice=matpar->GetRhoIce();
 	heatcapacity=matpar->GetHeatCapacity();
+	thermalconductivity=matpar->GetThermalConductivity();
+	this->parameters->FindParam(&dt,DtEnum);
+	this->parameters->FindParam(&artdiff,ArtDiffEnum);
 	Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
@@ -3771,4 +3798,5 @@
 	Input* temperature_input=NULL;
 	if (dt) temperature_input=inputs->GetInput(TemperatureEnum); ISSMASSERT(inputs);
+	if (artdiff==2) diameter=MinEdgeLength(xyz_list);
 
 	/* Start  looping on the number of gaussian points: */
@@ -3795,4 +3823,21 @@
 			scalar_transient=temperature*Jdet*gauss->weight;
 			for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_transient*L[i];
+		}
+
+		if(artdiff==2){
+			GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
+
+			vx_input->GetParameterValue(&u, gauss);
+			vy_input->GetParameterValue(&v, gauss);
+			vz_input->GetParameterValue(&w, gauss);
+
+			normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
+			if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){
+				tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity));
+			}
+			else tau_parameter=diameter/(2*normu);
+
+			scalar_artdiff=tau_parameter*gauss->weight*Jdet*phi/(rho_ice*heatcapacity);
+			for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_artdiff*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
 		}
 	}
@@ -5441,4 +5486,26 @@
 }
 /*}}}*/
+/*FUNCTION Penta::MinEdgeLength{{{1*/
+double Penta::MinEdgeLength(double xyz_list[6][3]){
+	/*Return the minimum lenght of the nine egdes of the penta*/
+
+	int    i,node0,node1;
+	int    edges[9][2]={{0,1},{0,2},{1,2},{3,4},{3,5},{4,5},{0,3},{1,4},{2,5}}; //list of the nine edges
+	double minlength=-1;
+	double length;
+	
+	for(i=0;i<9;i++){
+		/*Find the two nodes for this edge*/
+		node0=edges[i][0];
+		node1=edges[i][1];
+
+		/*Compute the length of this edge and compare it to the minimal length*/
+		length=pow(pow(xyz_list[node0][0]-xyz_list[node1][0],2.0)+pow(xyz_list[node0][1]-xyz_list[node1][1],2.0)+pow(xyz_list[node0][2]-xyz_list[node1][2],2.0),0.5);
+		if(length<minlength || minlength<0) minlength=length;
+	}
+
+	return minlength;
+}
+/*}}}*/
 /*FUNCTION Penta::ReduceMatrixStokes {{{1*/
 void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){
@@ -5604,5 +5671,4 @@
 	*(surface_normal+1)=normal[1]/normal_norm;
 	*(surface_normal+2)=normal[2]/normal_norm;
-
-}
-/*}}}*/
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 6211)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 6212)
@@ -218,5 +218,6 @@
 		bool    IsOnShelf(void); 
 		bool    IsOnWater(void); 
-		void	  ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp);
+		double  MinEdgeLength(double xyz_list[6][3]);
+	   void	  ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp);
 		void	  ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp);
 		void	  SetClone(int* minranks);
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 6211)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 6212)
@@ -2415,5 +2415,5 @@
 
 	/*Intermediaries */
-	bool       artdiff;
+	int        artdiff;
 	int        i,j,ig,dim;
 	double     Jdettria ,vx,vy,dvxdx,dvydy;
@@ -2586,5 +2586,5 @@
 
 	/*Intermediaries */
-	bool    artdiff;
+	int     artdiff;
 	int     i,j,ig,dim;
 	double  nx,ny,norm,Jdettria;
@@ -3157,5 +3157,5 @@
 
 	/*Intermediaries */
-	bool       artdiff;
+	int        artdiff;
 	int        i,j,ig,dim;
 	double     Jdettria,DL_scalar,dt;