Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 22074)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 22075)
@@ -22,4 +22,5 @@
 #define NUMNODESP1xP2 9
 #define NUMNODESP1xP3 12
+#define NUMNODESP1xP4 15
 #define NUMNODESP2xP1 12
 #define NUMNODESP2    18
@@ -76,4 +77,11 @@
 			break;
 		case P1xP3Enum:
+			numindices = 3;
+			indices    = xNew<int>(numindices);
+			indices[0] = 0;
+			indices[1] = 1;
+			indices[2] = 2;
+			break;
+		case P1xP4Enum:
 			numindices = 3;
 			indices    = xNew<int>(numindices);
@@ -423,5 +431,5 @@
 			basis[ 4]=(9.)/(16.)*gauss->coord2*(zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
 			basis[ 5]=(9.)/(16.)*gauss->coord3*(zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
-			/*quarter-sides of quads*/
+			/*third-sides of quads*/
 			basis[ 6]=(27.)/(16.)*gauss->coord1*(zeta-1)*(zeta-1./3.)*(zeta+1.);
 			basis[ 7]=(27.)/(16.)*gauss->coord2*(zeta-1)*(zeta-1./3.)*(zeta+1.);
@@ -430,4 +438,24 @@
 			basis[10]=-(27.)/(16.)*gauss->coord2*(zeta-1)*(zeta+1./3.)*(zeta+1.);
 			basis[11]=-(27.)/(16.)*gauss->coord3*(zeta-1)*(zeta+1./3.)*(zeta+1.);
+			return;
+		case P1xP4Enum :
+			/*Corner nodes*/
+			basis[ 0]=gauss->coord1*(2./3.)*(zeta-1.)*(zeta-0.5 )*(zeta)*(zeta+0.5);
+			basis[ 1]=gauss->coord2*(2./3.)*(zeta-1.)*(zeta-0.5 )*(zeta)*(zeta+0.5);
+			basis[ 2]=gauss->coord3*(2./3.)*(zeta-1.)*(zeta-0.5 )*(zeta)*(zeta+0.5);
+			basis[ 3]=gauss->coord1*(2./3.)*(zeta-0.5)*(zeta)*(zeta+0.5)*(zeta +1.);
+			basis[ 4]=gauss->coord2*(2./3.)*(zeta-0.5)*(zeta)*(zeta+0.5)*(zeta +1.);
+			basis[ 5]=gauss->coord3*(2./3.)*(zeta-0.5)*(zeta)*(zeta+0.5)*(zeta +1.);
+			/*mid-sides of quads (center of vertical edges)*/
+			basis[ 6]=gauss->coord1*4.*(zeta-1.)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
+			basis[ 7]=gauss->coord2*4.*(zeta-1.)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
+			basis[ 8]=gauss->coord3*4.*(zeta-1.)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
+			/*quarter-sides of quads (-0.5 and +0.5 of vertical edges)*/
+			basis[ 9]=gauss->coord1*(-8./3.)*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
+			basis[10]=gauss->coord2*(-8./3.)*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
+			basis[11]=gauss->coord3*(-8./3.)*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
+			basis[12]=gauss->coord1*(-8./3.)*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
+			basis[13]=gauss->coord2*(-8./3.)*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
+			basis[14]=gauss->coord3*(-8./3.)*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
 			return;
 		default:
@@ -964,4 +992,69 @@
 			dbasis[NUMNODESP1xP3*2+11] = -gauss->coord3 *(27./16.)*( 2.* zeta *( zeta + (1./3.)) + ( zeta - 1. )*( zeta + 1. ));
 			return;
+		case P1xP4Enum :
+			/*Nodal function 1*/
+			dbasis[NUMNODESP1xP4*0+0 ] = -0.5*(2./3.)     *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5);
+			dbasis[NUMNODESP1xP4*1+0 ] = -SQRT3/6.*(2./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5); 
+			dbasis[NUMNODESP1xP4*2+0 ] =  gauss->coord1 * 2./3.*( (2.*zeta-1)*(zeta -0.5)*(zeta +0.5) + 2.* zeta *zeta *(zeta -1.)); 
+			/*Nodal function 2*/
+			dbasis[NUMNODESP1xP4*0+1 ] = +0.5*(2./3.)     *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5);
+			dbasis[NUMNODESP1xP4*1+1 ] = -SQRT3/6.*(2./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5);
+			dbasis[NUMNODESP1xP4*2+1 ] = gauss->coord2* 2./3.* ((2.*zeta-1.)*(zeta -0.5)*(zeta +0.5) + 2. * zeta *zeta*(zeta -1.)); 
+			/*Nodal function 3*/
+			dbasis[NUMNODESP1xP4*0+2 ] = 0. ; 
+			dbasis[NUMNODESP1xP4*1+2 ] = SQRT3/3.*(2./3.)*(zeta -1.)*(zeta-0.5)*(zeta)*(zeta+0.5); 
+			dbasis[NUMNODESP1xP4*2+2 ] = gauss->coord3* 2./3.*( (2.*zeta-1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta -1.)); 
+			/*Nodal function 4*/
+			dbasis[NUMNODESP1xP4*0+3 ] = -0.5 *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. ); 
+			dbasis[NUMNODESP1xP4*1+3 ] = -SQRT3/6.*(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. ); 
+			dbasis[NUMNODESP1xP4*2+3 ] = gauss->coord1* 2./3.*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta +1.)); 
+			/*Nodal function 5*/
+			dbasis[NUMNODESP1xP4*0+4 ] = +0.5 *    (2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. ); 
+			dbasis[NUMNODESP1xP4*1+4 ] = -SQRT3/6.*(2./3.)*(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
+			dbasis[NUMNODESP1xP4*2+4 ] = gauss->coord2 * 2./3.*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta +1.)); 
+			/*Nodal function 6*/
+			dbasis[NUMNODESP1xP4*0+5 ] = 0. ; 
+			dbasis[NUMNODESP1xP4*1+5 ] = SQRT3/3.*(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. ); 
+			dbasis[NUMNODESP1xP4*2+5 ] = gauss->coord3 * 2./3.*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta +1)); 
+			/*Nodal function 7*/
+			dbasis[NUMNODESP1xP4*0+6 ] = -0.5 * 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP1xP4*1+6 ] = -SQRT3/6.* 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1.)  ; 
+			dbasis[NUMNODESP1xP4*2+6 ] = gauss->coord1* 4.*( 4.*zeta *zeta*zeta - (5./2.)*zeta ); 
+			/*Nodal function 8*/
+			dbasis[NUMNODESP1xP4*0+7 ] = +0.5* 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP1xP4*1+7 ] = -SQRT3/6.* 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. )  ; 
+			dbasis[NUMNODESP1xP4*2+7 ] = gauss->coord2* 4.*( 4.*zeta *zeta*zeta - (5./2.)*zeta ); 
+
+			/*Nodal function 9*/
+			dbasis[NUMNODESP1xP4*0+8 ] = 0. ; 
+			dbasis[NUMNODESP1xP4*1+8 ] = SQRT3/3. * 4.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. )  ; 
+			dbasis[NUMNODESP1xP4*2+8 ] = gauss->coord3* 4.*( 4.*zeta *zeta*zeta - (5./2.)*zeta ); 
+
+			/*Nodal function 10*/
+			dbasis[NUMNODESP1xP4*0+9 ] = -0.5* (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP1xP4*1+9 ] = -SQRT3/3.* (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP1xP4*2+9 ] = gauss->coord1* (-8./3.)*((2.*zeta -1.)*(zeta-0.5)*(zeta +1.) +zeta*(zeta -1.)*( 2.*zeta + 0.5)); 
+			/*Nodal function 11*/
+			dbasis[NUMNODESP1xP4*0+10] = +0.5* (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP1xP4*1+10] = -SQRT3/6.* (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP1xP4*2+10] = gauss->coord2* (-8./3.)*((2.*zeta -1.)*(zeta-0.5)*(zeta +1.) +zeta *(zeta -1.)*( 2.*zeta + 0.5)); 
+			/*Nodal function 12*/
+			dbasis[NUMNODESP1xP4*0+11] = 0. ; 
+			dbasis[NUMNODESP1xP4*1+11] = SQRT3/6.* (-8./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP1xP4*2+11] = gauss->coord3* (-8./3.)*((2.*zeta-1.)*(zeta-0.5)*(zeta +1.) +zeta *(zeta -1.)*( 2.*zeta + 0.5));
+			/*Nodal function 13*/
+			dbasis[NUMNODESP1xP4*0+12] = -0.5* (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP1xP4*1+12] = -SQRT3/6.* (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP1xP4*2+12] = gauss->coord1* (-8./3.)*((2.*zeta -1. ) *(zeta+0.5)* (zeta +1.) +  zeta* (zeta -1.)*( 2.*zeta + 3./2.));
+			/*Nodal function 14*/
+			dbasis[NUMNODESP1xP4*0+13] = +0.5* (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP1xP4*1+13] = -SQRT3/6.* (-8./3.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP1xP4*2+13] = gauss->coord2* (-8./3.)*((2.*zeta -1. )*(zeta+0.5)*(zeta +1.) +  zeta*(zeta -1.)*( 2.*zeta + 3./2.)); 
+			/*Nodal function 15*/
+			dbasis[NUMNODESP1xP4*0+14] = 0 ; 
+			dbasis[NUMNODESP1xP4*1+14] = SQRT3/6.* (-8./3.)*(zeta - 1.)*(zeta + 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP1xP4*2+14] = gauss->coord3* (-8./3.)*((2. *zeta -1. )*(zeta+0.5)*(zeta +1.) +  zeta*(zeta -1.)*( 2.*zeta + 3./2.));
+
+			return;
 		default:
 			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
@@ -1069,4 +1162,11 @@
 			break;
 		case P1xP3Enum:
+			numseg = 3;
+			indices = xNew<int>(numseg*2);
+			indices[0*2 + 0] = 0; indices[0*2 + 1] = 6;
+			indices[1*2 + 0] = 1; indices[1*2 + 1] = 7;
+			indices[2*2 + 0] = 2; indices[2*2 + 1] = 8;
+			break;
+		case P1xP4Enum:
 			numseg = 3;
 			indices = xNew<int>(numseg*2);
@@ -1130,4 +1230,5 @@
 		case P2xP4Enum:             return NUMNODESP2xP4;
 		case P1xP3Enum:             return NUMNODESP1xP3;
+		case P1xP4Enum:             return NUMNODESP1xP4;
 		case P1P1Enum:              return NUMNODESP1*2;
 		case P1P1GLSEnum:           return NUMNODESP1*2;
@@ -1196,4 +1297,5 @@
 		case P1xP2Enum:
 		case P1xP3Enum:
+		case P1xP4Enum:
 			numindices = 3;
 			indices    = xNew<int>(numindices);
Index: /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp	(revision 22074)
+++ /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp	(revision 22075)
@@ -254,4 +254,16 @@
 			}
 			break;
+		case P1xP4Enum:
+			if(start==-1){
+				for(int i=0+3;i<15;i+=3) this->values[i]=this->values[0];
+				for(int i=1+3;i<15;i+=3) this->values[i]=this->values[1];
+				for(int i=2+3;i<15;i+=3) this->values[i]=this->values[2];
+			}
+			else{
+				for(int i=0;i<12;i+=3) this->values[i]=this->values[12];
+				for(int i=1;i<12;i+=3) this->values[i]=this->values[13];
+				for(int i=2;i<12;i+=3) this->values[i]=this->values[14];
+			}
+			break;
 		default:
 			_error_("not supported yet for type "<<EnumToStringx(this->interpolation_type));
Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp	(revision 22074)
+++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp	(revision 22075)
@@ -601,4 +601,27 @@
 			}
 			break;
+		case P1xP4Enum: 
+			switch(iv){
+				case 0 : coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
+				case 1 : coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
+				case 2 : coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
+				case 3 : coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
+				case 4 : coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
+				case 5 : coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
+
+				case 6 : coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
+				case 7 : coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
+				case 8 : coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
+
+				case 9 : coord1=1.; coord2=0.; coord3=0.; coord4=-0.5; break;
+				case 10: coord1=0.; coord2=1.; coord3=0.; coord4=-0.5; break;
+				case 11: coord1=0.; coord2=0.; coord3=1.; coord4=-0.5; break;
+
+				case 12: coord1=1.; coord2=0.; coord3=0.; coord4=+0.5; break;
+				case 13: coord1=0.; coord2=1.; coord3=0.; coord4=+0.5; break;
+				case 14: coord1=0.; coord2=0.; coord3=1.; coord4=+0.5; break;
+				default: _error_("node index should be in [0 14]");
+			}
+			break;
 		case P2xP1Enum: 
 			switch(iv){
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 22074)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 22075)
@@ -1067,4 +1067,5 @@
 	P1xP2Enum,
 	P1xP3Enum,
+	P1xP4Enum,
 	P2xP4Enum,
 	P1P1Enum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 22074)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 22075)
@@ -1022,4 +1022,5 @@
 		case P1xP2Enum : return "P1xP2";
 		case P1xP3Enum : return "P1xP3";
+		case P1xP4Enum : return "P1xP4";
 		case P2xP4Enum : return "P2xP4";
 		case P1P1Enum : return "P1P1";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 22074)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 22075)
@@ -1046,4 +1046,5 @@
 	      else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
 	      else if (strcmp(name,"P1xP3")==0) return P1xP3Enum;
+	      else if (strcmp(name,"P1xP4")==0) return P1xP4Enum;
 	      else if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
 	      else if (strcmp(name,"P1P1")==0) return P1P1Enum;