Index: /issm/trunk/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk/src/c/classes/Elements/PentaRef.cpp	(revision 17015)
+++ /issm/trunk/src/c/classes/Elements/PentaRef.cpp	(revision 17016)
@@ -22,4 +22,6 @@
 #define NUMNODESP2xP1 12
 #define NUMNODESP2    15
+#define NUMNODESP2xP4 30
+#define NUMNODESP1xP3 12
 
 /*Object constructors and destructor*/
@@ -1333,4 +1335,69 @@
 			basis[14]=2.*gauss->coord1*gauss->coord2*zeta*(zeta+1.);
 			return;
+			
+			/*
+						MODIFICATION
+			*/
+		case P2xP4Enum :
+			/*Corner nodes*/
+			basis[ 0]=(2./3.)*gauss->coord1*(2.*gauss->coord1-1.)*(zeta-1.)*(zeta-0.5 )*(zeta)*(zeta+0.5);
+			basis[ 1]=(2./3.)*gauss->coord2*(2.*gauss->coord2-1.)*(zeta-1.)*(zeta-0.5 )*(zeta)*(zeta+0.5);
+			basis[ 2]=(2./3.)*gauss->coord3*(2.*gauss->coord3-1.)*(zeta-1.)*(zeta-0.5 )*(zeta)*(zeta+0.5);
+			basis[ 3]=(2./3.)*gauss->coord1*(2.*gauss->coord1-1.)*(zeta-0.5)*(zeta)*(zeta+0.5)*(zeta +1.);
+			basis[ 4]=(2./3.)*gauss->coord2*(2.*gauss->coord2-1.)*(zeta-0.5)*(zeta)*(zeta+0.5)*(zeta +1.);
+			basis[ 5]=(2./3.)*gauss->coord3*(2.*gauss->coord3-1.)*(zeta-0.5)*(zeta)*(zeta+0.5)*(zeta +1.);
+		
+			/*mid-sides of quads*/
+			basis[ 6]=4.*gauss->coord1*(zeta-1.)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
+			basis[ 7]=4.*gauss->coord2*(zeta-1.)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
+			basis[ 8]=4.*gauss->coord3*(zeta-1.)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
+			
+			/*mid-sides of triangles*/
+			basis[ 9]=(2./3.)*gauss->coord2*gauss->coord3*(zeta-1.)*(zeta-0.5)*zeta*(zeta+0.5);
+			basis[10]=(2./3.)*gauss->coord1*gauss->coord3*(zeta-1.)*(zeta-0.5)*zeta*(zeta+0.5);
+			basis[11]=(2./3.)*gauss->coord1*gauss->coord2*(zeta-1.)*(zeta-0.5)*zeta*(zeta+0.5);
+			basis[12]=(2./3.)*gauss->coord2*gauss->coord3*(zeta-0.5)*zeta*(zeta+0.5)*(zeta+1.);
+			basis[13]=(2./3.)*gauss->coord1*gauss->coord3*(zeta-0.5)*zeta*(zeta+0.5)*(zeta+1.);
+			basis[14]=(2./3.)*gauss->coord1*gauss->coord2*(zeta-0.5)*zeta*(zeta+0.5)*(zeta+1.);
+		
+			/*quarter-sides of quads*/
+			basis[15]=-(2./3.)*gauss->coord1*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
+			basis[16]=-(2./3.)*gauss->coord2*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
+			basis[17]=-(2./3.)*gauss->coord3*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
+			basis[18]=-(2./3.)*gauss->coord1*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
+			basis[19]=-(2./3.)*gauss->coord2*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
+			basis[20]=-(2./3.)*gauss->coord3*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
+
+			/* mid-sides of interior triangles*/
+			basis[21]=-(8./3.)*gauss->coord2*gauss->coord3*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
+			basis[22]=-(8./3.)*gauss->coord1*gauss->coord3*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
+			basis[23]=-(8./3.)*gauss->coord1*gauss->coord2*(zeta-1.0)*(zeta-0.5)*zeta*(zeta+1.);
+			basis[24]=16.*gauss->coord2*gauss->coord3*(zeta-1.0)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
+			basis[25]=16.*gauss->coord1*gauss->coord3*(zeta-1.0)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
+			basis[26]=16.*gauss->coord1*gauss->coord2*(zeta-1.0)*(zeta-0.5)*(zeta+0.5)*(zeta+1.);
+			basis[27]=-(8./3.)*gauss->coord2*gauss->coord3*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
+			basis[28]=-(8./3.)*gauss->coord1*gauss->coord3*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
+			basis[29]=-(8./3.)*gauss->coord1*gauss->coord2*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
+			return;
+			
+			
+		case P1xP3Enum :
+			/*Corner nodes*/
+			basis[ 0]=-(9.)/(16.)*gauss->coord1*(zeta-1)*(zeta-1./3.)*(zeta+1./3.);	
+			basis[ 1]=-(9.)/(16.)*gauss->coord2*(zeta-1)*(zeta-1./3.)*(zeta+1./3.);	
+			basis[ 2]=-(9.)/(16.)*gauss->coord3*(zeta-1)*(zeta-1./3.)*(zeta+1./3.);	
+			basis[ 3]=(9.)/(16.)*gauss->coord1*(zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
+			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*/
+			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.);
+			basis[ 8]=(27.)/(16.)*gauss->coord3*(zeta-1)*(zeta-1./3.)*(zeta+1.);
+			basis[ 9]=-(27.)/(16.)*gauss->coord1*(zeta-1)*(zeta+1./3.)*(zeta+1.);
+			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;
+			
 		default:
 			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
@@ -1731,4 +1798,182 @@
 			dbasis[NUMNODESP2*2+14] = 2.*gauss->coord1*gauss->coord2*(2.*zeta+1.);
 			return;
+			
+			/*
+			Modification
+			*/
+		case P2xP4Enum :
+			/*Nodal function 1*/
+			dbasis[NUMNODESP2xP4*0+0 ] = (-2* gauss->coord1 + 0.5 ) *(2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ; 
+			dbasis[NUMNODESP2xP4*1+0 ] = (-((2.*SQRT3)/(3.))*gauss->coord1 + (SQRT3/6.) )*(2./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ; 
+			dbasis[NUMNODESP2xP4*2+0 ] = (2./3.)* gauss->coord1 *(2.* gauss->coord1 -1)*( (2.*zeta-1)*(zeta -0.5)*(zeta +0.5) + 2.* zeta *zeta *(zeta -1.)); 
+      /*Nodal function 2*/
+			dbasis[NUMNODESP2xP4*0+1 ] = (2.*gauss->coord2 - 0.5 ) *(2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ; 
+			dbasis[NUMNODESP2xP4*1+1 ] = (-((2.*SQRT3)/(3.))*gauss->coord2 + (SQRT3/6.) )*(2./3.)*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ; 
+      dbasis[NUMNODESP2xP4*2+1 ] = (2./3.)*gauss->coord2*(2.*gauss->coord2 -1.)*((2.*zeta-1.)*(zeta -0.5)*(zeta +0.5) + 2. * zeta *zeta*(zeta -1.)); 
+      /*Nodal function 3*/
+			dbasis[NUMNODESP2xP4*0+2 ] = 0. ; 
+			dbasis[NUMNODESP2xP4*1+2 ] = (((4.*SQRT3)/(3.))*gauss->coord3 - (SQRT3)/(3.))*(2./3.)*(zeta -1.)*(zeta-0.5)*(zeta)*(zeta+0.5); 
+      dbasis[NUMNODESP2xP4*2+2 ] = (2./3.)*gauss->coord3*(2.* gauss->coord3 -1.)*( (2.*zeta-1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta -1.)); 
+      /*Nodal function 4*/
+			dbasis[NUMNODESP2xP4*0+3 ] = (-2.* gauss->coord1 + 0.5 ) *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. ); 
+			dbasis[NUMNODESP2xP4*1+3 ] = (-((2.*SQRT3)/(3.)) *gauss->coord1 + (SQRT3)/(6.) ) *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. ); 
+      dbasis[NUMNODESP2xP4*2+3 ] = (2./3.)*gauss->coord1*(2.*gauss->coord1 -1.)*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta +1.)); 
+      /*Nodal function 5*/
+			dbasis[NUMNODESP2xP4*0+4 ] = (2*gauss->coord2 - 0.5 ) *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. ); 
+			dbasis[NUMNODESP2xP4*1+4 ] = -((2.*SQRT3)/(3.)) *gauss->coord2 + (SQRT3/(6.))*(2./3.)*(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
+      dbasis[NUMNODESP2xP4*2+4 ] = (2./3.)*gauss->coord2 *(2.*gauss->coord2 -1.)*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta +1.)); 
+      /*Nodal function 6*/
+			dbasis[NUMNODESP2xP4*0+5 ] = 0. ; 
+			dbasis[NUMNODESP2xP4*1+5 ] = (((4.*SQRT3)/(3.))*gauss->coord3 - SQRT3/3. ) *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. ); 
+      dbasis[NUMNODESP2xP4*2+5 ] = (2./3.)*gauss->coord3 *(2.*gauss->coord3 -1.)*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta +1)); 
+      /*Nodal function 7*/
+			dbasis[NUMNODESP2xP4*0+6 ] =  -2. *(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*1+6 ] =  (-4.*SQRT3)/(6.) *(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1.)  ; 
+      dbasis[NUMNODESP2xP4*2+6 ] =  4.*gauss->coord1*( 4.*zeta *zeta*zeta - (5./2.)*zeta ); 
+      /*Nodal function 8*/
+			dbasis[NUMNODESP2xP4*0+7 ] =  2.*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*1+7 ] =  ((-4.*SQRT3)/(6.)) *(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. )  ; 
+      dbasis[NUMNODESP2xP4*2+7 ] =  4.* gauss->coord2*( 4.*zeta *zeta*zeta - (5./2.)*zeta ); 
+      /*Nodal function 9*/
+			dbasis[NUMNODESP2xP4*0+8 ] = 0. ; 
+			dbasis[NUMNODESP2xP4*1+8 ] = ((4.*SQRT3)/(3.)) *(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. )  ; 
+			dbasis[NUMNODESP2xP4*2+8 ] =  4.*gauss->coord3*( 4.*zeta *zeta*zeta - (5./2.)*zeta ); 
+      /*Nodal function 10*/
+			dbasis[NUMNODESP2xP4*0+9 ] = (4./3.)*(gauss->coord3 )*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ; 
+			dbasis[NUMNODESP2xP4*1+9 ] =  4.* (((- SQRT3)/6.)* gauss->coord3+ (SQRT3/3.) *gauss->coord2) *(2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5);
+			dbasis[NUMNODESP2xP4*2+9 ] = 4.* gauss->coord2 * gauss->coord3 *(2./3.)*((2.*zeta-1.)*(zeta -0.5)*(zeta +0.5) + 2.* zeta*zeta*(zeta -1.)); 
+      /*Nodal function 11*/
+			dbasis[NUMNODESP2xP4*0+10] = -(4./3.)* gauss->coord3*(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ; 
+			dbasis[NUMNODESP2xP4*1+10] =  4.* ((SQRT3/3.)*gauss->coord1- (SQRT3/6.)*gauss->coord3) *(2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5);
+			dbasis[NUMNODESP2xP4*2+10] = 4.* gauss->coord3*gauss->coord1 *(2./3.)*( (2*zeta-1.)*(zeta -0.5)*(zeta +0.5) + 2* zeta * zeta*(zeta -1));
+      /*Nodal function 12*/
+			dbasis[NUMNODESP2xP4*0+11] = 4.* (gauss->coord1*0.5- gauss->coord2 *0.5 )*(2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ; 
+			dbasis[NUMNODESP2xP4*1+11] =  4.* (- (SQRT3/6.) *gauss->coord2 - (SQRT3/6.) *gauss->coord1) *(2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5);
+			dbasis[NUMNODESP2xP4*2+11] = 4.* gauss->coord1*gauss->coord2 *(2./3.) *( (2.*zeta-1)*(zeta -0.5)*(zeta +0.5) + 2* zeta* zeta*(zeta -1));
+      /*Nodal function 13*/
+			dbasis[NUMNODESP2xP4*0+12] = (4./3.) * gauss->coord3 *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*1+12] =  4.* ( (- SQRT3/6. ) * gauss->coord3 + (SQRT3/3.) *gauss->coord2) *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
+			dbasis[NUMNODESP2xP4*2+12] = 4.* gauss->coord2*gauss->coord3 *(2./3.)*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2*zeta*zeta*(zeta +1.)); 
+      /*Nodal function 14*/
+			dbasis[NUMNODESP2xP4*0+13] = (- 4./3.) *gauss->coord3*(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*1+13] =  4.* ((SQRT3/3.)*gauss->coord1- (SQRT3/6.)*gauss->coord3) *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
+			dbasis[NUMNODESP2xP4*2+13] = 4.* gauss->coord3*gauss->coord1 *(2./3.)*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2.*zeta *zeta*(zeta +1.)); 
+      /*Nodal function 15*/
+			dbasis[NUMNODESP2xP4*0+14] = 4.* ( gauss->coord1*0.5- gauss->coord2*0.5) *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*1+14] =  4.* (- (SQRT3/6.)*gauss->coord2 - (SQRT3/6.)*gauss->coord1) *(2./3.) *(zeta - 0.5)*(zeta)*(zeta+0.5)*(zeta +1. );
+			dbasis[NUMNODESP2xP4*2+14] = 4.* gauss->coord1*gauss->coord2 *(2./3.)*( (2.*zeta+1.)*(zeta -0.5)*(zeta +0.5) + 2* zeta *zeta*(zeta +1.)); 
+      /*Nodal function 16*/
+			dbasis[NUMNODESP2xP4*0+15] =  (8./6.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*1+15] =  (8./3.) *(SQRT3/6.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+15] =  -(8./3.) *gauss->coord1 *((2.*zeta -1. )*(zeta-0.5)*(zeta +1.) +  zeta*(zeta -1.)*( 2.*zeta + 0.5)); 
+      /*Nodal function 17*/
+			dbasis[NUMNODESP2xP4*0+16] = - (8./6.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*1+16] =  (8./3.) * (SQRT3/6.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+16] =  -(8./3.)*gauss->coord2 *((2.*zeta -1. )*(zeta-0.5)*(zeta +1.) +  zeta *(zeta -1.)*( 2.*zeta + 0.5)); 
+      /*Nodal function 18*/
+			dbasis[NUMNODESP2xP4*0+17] = 0. ; 
+			dbasis[NUMNODESP2xP4*1+17] =  (8./3.) *(SQRT3/3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+17] =  -(8./3.)*gauss->coord3 *((2.* zeta -1. ) *(zeta-0.5)* (zeta +1.) +  zeta *(zeta -1.)*( 2.*zeta + 0.5));
+      /*Nodal function 19*/
+			dbasis[NUMNODESP2xP4*0+18] =  (4./3.) *(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*1+18] =  (8./3.)* (SQRT3/6.)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+18] =  -(8./3.)*gauss->coord1*((2.*zeta -1. ) *(zeta+0.5)* (zeta +1.) +  zeta* (zeta -1.)*( 2.*zeta + 3./2.));
+      /*Nodal function 20*/
+			dbasis[NUMNODESP2xP4*0+19] =  (8./6.) *(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*1+19] =  (8./3.)*(SQRT3/6) *(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+19] =  -(8./3.)*gauss->coord2*((2.*zeta -1. )*(zeta+0.5)*(zeta +1.) +  zeta*(zeta -1.)*( 2.*zeta + 3./2.)); 
+      /*Nodal function 21*/
+			dbasis[NUMNODESP2xP4*0+20] = 0 ; 
+			dbasis[NUMNODESP2xP4*1+20] =  -(8./3.)* (SQRT3/3) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+20] =  -(8./3.)*gauss->coord3 *(2. *zeta -1. )*(zeta+0.5)*(zeta +1.) +  zeta*(zeta -1.)*( 2.*zeta + 3./2.);
+      /*Nodal function 22*/
+			dbasis[NUMNODESP2xP4*0+21] = -(32./6.) *gauss->coord3 *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ); 
+			dbasis[NUMNODESP2xP4*1+21] = -(32./3.)*(- (SQRT3/6.)*gauss->coord3+ (SQRT3/3.)*gauss->coord2 ) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+21] = -(32./3.)*gauss->coord2 *gauss->coord3*((2.*zeta -1. )*(zeta-0.5)*(zeta +1.) +  zeta*(zeta -1.)*( 2.*zeta + 0.5)); 
+      /*Nodal function 23*/
+			dbasis[NUMNODESP2xP4*0+22] = (32./6.) *gauss->coord3 *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ); 
+			dbasis[NUMNODESP2xP4*1+22] = -(32./3.)*( (SQRT3/3.)*gauss->coord1- (SQRT3/6.)*gauss->coord3 )*(zeta -1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+22] = -(32./3.) *gauss->coord1*gauss->coord3*((2.*zeta -1. )*(zeta-0.5)*(zeta +1.) +  zeta*(zeta -1.)*( 2.*zeta + 0.5)); 
+      /*Nodal function 24*/
+			dbasis[NUMNODESP2xP4*0+23] = -(32./3.)*( gauss->coord1*0.5- gauss->coord2*0.5) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ); 
+			dbasis[NUMNODESP2xP4*1+23] = -(32./3.)*(- (SQRT3/6.)*gauss->coord2- (SQRT3/6.)*gauss->coord1 ) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+23] = -(32./3.)*gauss->coord1*gauss->coord2*((2.*zeta -1. )* (zeta-0.5) *(zeta +1.) +  zeta* (zeta -1.)*( 2.*zeta + 0.5));
+      /*Nodal function 25*/
+			dbasis[NUMNODESP2xP4*0+24] = 8. *gauss->coord3 *(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ); 
+			dbasis[NUMNODESP2xP4*1+24] = 16.*(- (SQRT3/6.)*gauss->coord3 + (SQRT3/3.) *gauss->coord2) *(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+24] = 16.* gauss->coord2 * gauss->coord3* ( 4.* zeta*zeta*zeta - (5./2.) *zeta ); 
+      /*Nodal function 26*/
+			dbasis[NUMNODESP2xP4*0+25] = 8. *gauss->coord3*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ); 
+			dbasis[NUMNODESP2xP4*1+25] = -16.*((SQRT3/3.)*gauss->coord1- (SQRT3/6.)*gauss->coord3 )*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+25] = 16. * gauss->coord1 * gauss->coord3 *( 4. *zeta*zeta*zeta - (5./2.)* zeta );
+      /*Nodal function 27*/
+			dbasis[NUMNODESP2xP4*0+26] = 16.*( gauss->coord1*0.5-gauss->coord2*0.5) *(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ); 
+			dbasis[NUMNODESP2xP4*1+26] = 16.*(- (SQRT3/6.)*gauss->coord2- (SQRT3/6.)*gauss->coord1 )*(zeta - 1.)*(zeta - 0.5)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+26] = 16. *gauss->coord1 *gauss->coord2 *( 4.* zeta*zeta*zeta - (5./2.)*zeta ); 
+      /*Nodal function 28*/
+			dbasis[NUMNODESP2xP4*0+27] = -(32./6.)* gauss->coord3*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ); 
+			dbasis[NUMNODESP2xP4*1+27] = -(32./3.)*(- (SQRT3/6.)*gauss->coord3+ (SQRT3/3.)*gauss->coord2 ) *(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+27] = -(32./3.)* gauss->coord2*gauss->coord3* ((2.*zeta -1. )*(zeta+0.5)*(zeta +1.) +  zeta*(zeta -1.)*( 2.*zeta + (3./2.))); 
+      /*Nodal function 29*/
+			dbasis[NUMNODESP2xP4*0+28] = (32./6.) *gauss->coord3 *(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ); 
+			dbasis[NUMNODESP2xP4*1+28] = -(32./3.)*((SQRT3/3.)*gauss->coord1- (SQRT3/6.)*gauss->coord3 ) *(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+28] = -(32./3.)* gauss->coord1*gauss->coord3* ((2.*zeta -1. )*(zeta+0.5)*(zeta +1.) +  zeta*(zeta -1.)*( 2.*zeta + (3./2.))); 
+      /*Nodal function 30*/
+			dbasis[NUMNODESP2xP4*0+29] = -(32./3.)*( gauss->coord1*0.5- gauss->coord2*0.5)*(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ); 
+			dbasis[NUMNODESP2xP4*1+29] = -(32./3.)*(- (SQRT3/6.)*gauss->coord2- (SQRT3/6.)*gauss->coord1 ) *(zeta - 1.)*(zeta)*(zeta+0.5)*(zeta +1. ) ; 
+			dbasis[NUMNODESP2xP4*2+29] = -(32./3.)* gauss->coord1*gauss->coord2 *((2.*zeta -1. )*(zeta+0.5)*(zeta +1) +   zeta*(zeta -1.)*( 2.*zeta + (3./2.))); 
+			return;
+			
+		case P1xP3Enum :
+			/*Nodal function 1*/
+	    dbasis[NUMNODESP1xP3*0+0 ] =  (9./32.)*(zeta-1)*(zeta-1./3.)*(zeta+1./3.);
+	    dbasis[NUMNODESP1xP3*1+0 ] = ((3.*SQRT3)/32.)*(zeta-1)*(zeta-1./3.)*(zeta+1./3.);
+	    dbasis[NUMNODESP1xP3*2+0 ] =- (9./16.)* gauss->coord1 *( 2. *zeta *( zeta -1. ) + ( zeta - (1./3.) )*( zeta + (1./3.) ));
+			/*Nodal function 2*/
+	    dbasis[NUMNODESP1xP3*0+1 ] = - (9./32.)*(zeta-1)*(zeta-1./3.)*(zeta+1./3.);
+	    dbasis[NUMNODESP1xP3*1+1 ] = ((3.*SQRT3)/32.) *(zeta-1)*(zeta-1./3.)*(zeta+1./3.);
+	    dbasis[NUMNODESP1xP3*2+1 ] =- (9./16.)*gauss->coord2 *( 2.* zeta* ( zeta -1. ) + ( zeta - (1./3.) )*( zeta + (1./3.) ));
+		  /*Nodal function 3*/
+	    dbasis[NUMNODESP1xP3*0+2 ] =  0.;
+	    dbasis[NUMNODESP1xP3*1+2 ] = - ((3.*SQRT3)/32.)*(zeta-1)*(zeta-1./3.)*(zeta+1./3.);
+	    dbasis[NUMNODESP1xP3*2+2 ] = - (9./16.)* gauss->coord3* ( 2. *zeta *( zeta -1. ) + ( zeta - (1./3.) )*( zeta + (1./3.) ));
+			/*Nodal function 4*/	 
+	    dbasis[NUMNODESP1xP3*0+3 ] = -  (9./32.)*(zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
+	    dbasis[NUMNODESP1xP3*1+3 ] =  -((3.*SQRT3)/32.) *(zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
+	    dbasis[NUMNODESP1xP3*2+3 ] = (9./16.)* gauss->coord1*( 2.* zeta *( zeta +1. ) + ( zeta - (1./3.) )*( zeta + (1./3.) ));
+			/*Nodal function 5*/	
+	    dbasis[NUMNODESP1xP3*0+4 ] =   (9./32.)* (zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
+	    dbasis[NUMNODESP1xP3*1+4 ] = - ((3.*SQRT3)/32.) *(zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
+	    dbasis[NUMNODESP1xP3*2+4 ] = (9./16.)* gauss->coord2* ( 2.* zeta *( zeta +1. ) + ( zeta - (1./3.) )*( zeta + (1./3.) ));
+			/*Nodal function 6*/	
+	    dbasis[NUMNODESP1xP3*0+5 ] =  0.;
+	    dbasis[NUMNODESP1xP3*1+5 ] =  ((3.*SQRT3)/16.)  *(zeta-1./3.)*(zeta+1./3.)*(zeta+1.);
+	    dbasis[NUMNODESP1xP3*2+5 ] =  (9./16.)* gauss->coord3 *( 2.* zeta * ( zeta  + 1. ) + ( zeta - (1./3.) )*( zeta + (1./3.) ));
+			/*Nodal function 7*/	
+	    dbasis[NUMNODESP1xP3*0+6 ] = -  (27./32.) *(zeta-1)*(zeta-1./3.)*(zeta+1.);
+	    dbasis[NUMNODESP1xP3*1+6 ] = -  (27./32.) *(zeta-1)*(zeta-1./3.)*(zeta+1.);
+	    dbasis[NUMNODESP1xP3*2+6 ] =  gauss->coord1*(27./16.)*( 2.* zeta *( zeta - (1./3.)) + ( zeta - 1. )*( zeta + 1. ));
+			/*Nodal function 8*/	
+	    dbasis[NUMNODESP1xP3*0+7 ] =  (27./32.) *(zeta-1)*(zeta-1./3.)*(zeta+1.);
+	    dbasis[NUMNODESP1xP3*1+7 ] = -(27./32.) *(zeta-1)*(zeta-1./3.)*(zeta+1.);
+	    dbasis[NUMNODESP1xP3*2+7 ] =  gauss->coord2*(27./16.)*( 2.* zeta *( zeta - (1./3.)) + ( zeta - 1. )*( zeta + 1. ));
+			/*Nodal function 9*/	
+	    dbasis[NUMNODESP1xP3*0+8 ] = 0.;
+	    dbasis[NUMNODESP1xP3*1+8 ] =  ((9.*SQRT3)/16.) *(zeta-1.)*(zeta-1./3.)*(zeta+1.);
+	    dbasis[NUMNODESP1xP3*2+8 ] =  gauss->coord3*(27./16.)*( 2. *zeta *( zeta - (1./3.)) + ( zeta - 1. )*( zeta + 1. ));
+			/*Nodal function 10*/	
+	    dbasis[NUMNODESP1xP3*0+9 ] = - (27./32.) *(zeta-1.)*(zeta+1./3.)*(zeta+1.);
+	    dbasis[NUMNODESP1xP3*1+9 ] = - ((9.*SQRT3)/32.) *(zeta-1.)*(zeta+1./3.)*(zeta+1.);
+	    dbasis[NUMNODESP1xP3*2+9 ] =  gauss->coord1 *(27./16.)*( 2* zeta *( zeta + (1./3.)) + ( zeta - 1. )*( zeta + 1. ));
+			/*Nodal function 11*/	
+	    dbasis[NUMNODESP1xP3*0+10] =  (27./32.) *(zeta-1)*(zeta+1./3.)*(zeta+1);
+	    dbasis[NUMNODESP1xP3*1+10] = - ((9.*SQRT3)/32.)  *(zeta-1.)*(zeta+1./3.)*(zeta+1);
+	    dbasis[NUMNODESP1xP3*2+10] =  gauss->coord2 *(27./16.) * 2.* zeta *( zeta + (1./3.) + ( zeta - 1. )*( zeta + 1. ));
+			/*Nodal function 12*/	
+	    dbasis[NUMNODESP1xP3*0+11] = 0.;
+	    dbasis[NUMNODESP1xP3*1+11] =  ((9.*SQRT3)/16.) *(zeta-1.)*(zeta+1./3.)*(zeta+1);
+	    dbasis[NUMNODESP1xP3*2+11] =  gauss->coord3 *(27./16.)* 2.* zeta *( zeta + (1./3.) + ( zeta - 1. )*( zeta + 1. ));
+			return;
+			
 		default:
 			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
@@ -1973,4 +2218,6 @@
 		case MINIEnum:              return NUMNODESP1b+NUMNODESP1;
 		case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
+		case P2xP4Enum:				return NUMNODESP2xP4;
+		case P1xP3Enum:				return NUMNODESP1xP3;
 		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
 	}
Index: /issm/trunk/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/shared/Enum/EnumDefinitions.h	(revision 17015)
+++ /issm/trunk/src/c/shared/Enum/EnumDefinitions.h	(revision 17016)
@@ -524,4 +524,6 @@
 	MINIcondensedEnum,
 	TaylorHoodEnum,
+	P2xP4Enum,
+	P1xP3Enum,
 	/*}}}*/
 	/*Results{{{*/
Index: /issm/trunk/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk/src/c/shared/Enum/EnumToStringx.cpp	(revision 17015)
+++ /issm/trunk/src/c/shared/Enum/EnumToStringx.cpp	(revision 17016)
@@ -516,4 +516,6 @@
 		case MINIcondensedEnum : return "MINIcondensed";
 		case TaylorHoodEnum : return "TaylorHood";
+		case P2xP4Enum : return "P2xP4";
+		case P1xP3Enum : return "P1xP3";
 		case SaveResultsEnum : return "SaveResults";
 		case BoolExternalResultEnum : return "BoolExternalResult";
Index: /issm/trunk/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk/src/c/shared/Enum/StringToEnumx.cpp	(revision 17015)
+++ /issm/trunk/src/c/shared/Enum/StringToEnumx.cpp	(revision 17016)
@@ -528,4 +528,6 @@
 	      else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
 	      else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
+	      else if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
+	      else if (strcmp(name,"P1xP3")==0) return P1xP3Enum;
 	      else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
 	      else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
Index: /issm/trunk/src/m/enum/EnumDefinitions.py
===================================================================
--- /issm/trunk/src/m/enum/EnumDefinitions.py	(revision 17015)
+++ /issm/trunk/src/m/enum/EnumDefinitions.py	(revision 17016)
@@ -508,4 +508,6 @@
 def MINIcondensedEnum(): return StringToEnum("MINIcondensed")[0]
 def TaylorHoodEnum(): return StringToEnum("TaylorHood")[0]
+def P2xP4Enum(): return StringToEnum("P2xP4")[0]
+def P1xP3Enum(): return StringToEnum("P1xP3")[0]
 def SaveResultsEnum(): return StringToEnum("SaveResults")[0]
 def BoolExternalResultEnum(): return StringToEnum("BoolExternalResult")[0]