Index: /issm/trunk-jpl/src/c/classes/objects/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/TriaRef.cpp	(revision 14766)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/TriaRef.cpp	(revision 14767)
@@ -21,5 +21,6 @@
 
 /*Element macros*/
-#define NUMNODES 3
+#define NUMNODESP1 3
+#define NUMNODESP2 6
 
 /*Object constructors and destructor*/
@@ -65,9 +66,9 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES)
+	 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
 	 */
 
 	int i;
-	IssmDouble dbasis[NDOF2][NUMNODES];
+	IssmDouble dbasis[NDOF2][NUMNODESP1];
 
 	/*Get dh1dh2dh3 in actual coordinate system: */
@@ -75,7 +76,7 @@
 
 	/*Build B: */
-	for (i=0;i<NUMNODES;i++){
-		B[NDOF1*NUMNODES*0+NDOF1*i]=dbasis[0][i]; 
-		B[NDOF1*NUMNODES*1+NDOF1*i]=dbasis[1][i]; 
+	for (i=0;i<NUMNODESP1;i++){
+		B[NDOF1*NUMNODESP1*0+NDOF1*i]=dbasis[0][i]; 
+		B[NDOF1*NUMNODESP1*1+NDOF1*i]=dbasis[1][i]; 
 	}
 }
@@ -91,9 +92,9 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES)
+	 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
 	 */
 
 	int i;
-	IssmDouble dbasis[NDOF2][NUMNODES];
+	IssmDouble dbasis[NDOF2][NUMNODESP1];
 
 	/*Get dh1dh2dh3 in actual coordinate system: */
@@ -101,11 +102,11 @@
 
 	/*Build B: */
-	for (i=0;i<NUMNODES;i++){
-		*(B+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i]; //B[0][NDOF2*i]=dbasis[0][i];
-		*(B+NDOF2*NUMNODES*0+NDOF2*i+1)=0;
-		*(B+NDOF2*NUMNODES*1+NDOF2*i)=0;
-		*(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i];
-		*(B+NDOF2*NUMNODES*2+NDOF2*i)=(float).5*dbasis[1][i]; 
-		*(B+NDOF2*NUMNODES*2+NDOF2*i+1)=(float).5*dbasis[0][i]; 
+	for (i=0;i<NUMNODESP1;i++){
+		*(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i]; //B[0][NDOF2*i]=dbasis[0][i];
+		*(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0;
+		*(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0;
+		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dbasis[1][i]; 
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][i]; 
 	}
 }
@@ -122,9 +123,9 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES)
+	 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
 	 */
 
 	/*Same thing in the actual coordinate system: */
-	IssmDouble dbasis[NDOF2][NUMNODES];
+	IssmDouble dbasis[NDOF2][NUMNODESP1];
 
 	/*Get dh1dh2dh3 in actual coordinates system : */
@@ -132,11 +133,11 @@
 
 	/*Build B': */
-	for (int i=0;i<NUMNODES;i++){
-		*(B+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i]; 
-		*(B+NDOF2*NUMNODES*0+NDOF2*i+1)=0; 
-		*(B+NDOF2*NUMNODES*1+NDOF2*i)=0; 
-		*(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i]; 
-		*(B+NDOF2*NUMNODES*2+NDOF2*i)=0.5*dbasis[1][i]; 
-		*(B+NDOF2*NUMNODES*2+NDOF2*i+1)=0.5*dbasis[0][i]; 
+	for (int i=0;i<NUMNODESP1;i++){
+		*(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i]; 
+		*(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0; 
+		*(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0; 
+		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i]; 
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i)=0.5*dbasis[1][i]; 
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=0.5*dbasis[0][i]; 
 	}
 }
@@ -151,5 +152,5 @@
 	 */
 
-	IssmDouble l1l3[NUMNODES];
+	IssmDouble l1l3[NUMNODESP1];
 
 	GetNodalFunctions(&l1l3[0],gauss);
@@ -170,5 +171,5 @@
 	 */
 
-	IssmDouble l1l3[NUMNODES];
+	IssmDouble l1l3[NUMNODESP1];
 
 	GetNodalFunctions(&l1l3[0],gauss);
@@ -189,8 +190,8 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*NUMNODES)
-	 */
-
-	IssmDouble basis[NUMNODES];
+	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*NUMNODESP1)
+	 */
+
+	IssmDouble basis[NUMNODESP1];
 
 	/*Get dh1dh2dh3 in actual coordinate system: */
@@ -198,7 +199,7 @@
 
 	/*Build B_prog: */
-	for (int i=0;i<NUMNODES;i++){
-		*(B_prog+NDOF1*NUMNODES*0+NDOF1*i)=basis[i];
-		*(B_prog+NDOF1*NUMNODES*1+NDOF1*i)=basis[i];
+	for (int i=0;i<NUMNODESP1;i++){
+		*(B_prog+NDOF1*NUMNODESP1*0+NDOF1*i)=basis[i];
+		*(B_prog+NDOF1*NUMNODESP1*1+NDOF1*i)=basis[i];
 	}
 }
@@ -215,9 +216,9 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODES)
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
 	 */
 
 	/*Same thing in the actual coordinate system: */
-	IssmDouble dbasis[NDOF2][NUMNODES];
+	IssmDouble dbasis[NDOF2][NUMNODESP1];
 
 	/*Get dh1dh2dh3 in actual coordinates system : */
@@ -225,11 +226,11 @@
 
 	/*Build B': */
-	for (int i=0;i<NUMNODES;i++){
-		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=2*dbasis[0][i]; 
-		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=dbasis[1][i]; 
-		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=dbasis[0][i]; 
-		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=2*dbasis[1][i]; 
-		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dbasis[1][i]; 
-		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dbasis[0][i]; 
+	for (int i=0;i<NUMNODESP1;i++){
+		*(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i)=2*dbasis[0][i]; 
+		*(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i+1)=dbasis[1][i]; 
+		*(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i)=dbasis[0][i]; 
+		*(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2*dbasis[1][i]; 
+		*(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i)=dbasis[1][i]; 
+		*(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dbasis[0][i]; 
 	}
 }
@@ -247,9 +248,9 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume Bprime has been allocated already, of size: 3x(NDOF2*NUMNODES)
+	 * We assume Bprime has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
 	 */
 
 	/*Same thing in the actual coordinate system: */
-	IssmDouble dbasis[NDOF2][NUMNODES];
+	IssmDouble dbasis[NDOF2][NUMNODESP1];
 
 	/*Get dh1dh2dh3 in actual coordinates system : */
@@ -257,13 +258,13 @@
 
 	/*Build Bprime: */
-	for (int i=0;i<NUMNODES;i++){
-		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i]; 
-		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=0; 
-		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=0; 
-		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i]; 
-		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dbasis[1][i]; 
-		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dbasis[0][i]; 
-		*(Bprime+NDOF2*NUMNODES*3+NDOF2*i)=dbasis[0][i]; 
-		*(Bprime+NDOF2*NUMNODES*3+NDOF2*i+1)=dbasis[1][i]; 
+	for (int i=0;i<NUMNODESP1;i++){
+		*(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i]; 
+		*(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0; 
+		*(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i)=0; 
+		*(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i]; 
+		*(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i)=dbasis[1][i]; 
+		*(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dbasis[0][i]; 
+		*(Bprime+NDOF2*NUMNODESP1*3+NDOF2*i)=dbasis[0][i]; 
+		*(Bprime+NDOF2*NUMNODESP1*3+NDOF2*i+1)=dbasis[1][i]; 
 	}
 }
@@ -278,9 +279,9 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODES)
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
 	 */
 
 	/*Same thing in the actual coordinate system: */
-	IssmDouble dbasis[NDOF2][NUMNODES];
+	IssmDouble dbasis[NDOF2][NUMNODESP1];
 
 	/*Get dh1dh2dh3 in actual coordinates system : */
@@ -288,7 +289,7 @@
 
 	/*Build B': */
-	for (int i=0;i<NUMNODES;i++){
-		*(Bprime_prog+NDOF1*NUMNODES*0+NDOF1*i)=dbasis[0][i]; 
-		*(Bprime_prog+NDOF1*NUMNODES*1+NDOF1*i)=dbasis[1][i]; 
+	for (int i=0;i<NUMNODESP1;i++){
+		*(Bprime_prog+NDOF1*NUMNODESP1*0+NDOF1*i)=dbasis[0][i]; 
+		*(Bprime_prog+NDOF1*NUMNODESP1*1+NDOF1*i)=dbasis[1][i]; 
 	}
 }
@@ -306,5 +307,5 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume L has been allocated already, of size: NUMNODES (numdof=1), or numdofx(numdof*NUMNODES) (numdof=2)
+	 * We assume L has been allocated already, of size: NUMNODESP1 (numdof=1), or numdofx(numdof*NUMNODESP1) (numdof=2)
 	 */
 
@@ -317,14 +318,14 @@
 	/*Build L: */
 	if(numdof==1){
-		for (i=0;i<NUMNODES;i++){
+		for (i=0;i<NUMNODESP1;i++){
 			L[i]=basis[i]; 
 		}
 	}
 	else{
-		for (i=0;i<NUMNODES;i++){
-			*(L+numdof*NUMNODES*0+numdof*i)=basis[i]; //L[0][NDOF2*i]=dbasis[0][i];
-			*(L+numdof*NUMNODES*0+numdof*i+1)=0;
-			*(L+numdof*NUMNODES*1+numdof*i)=0;
-			*(L+numdof*NUMNODES*1+numdof*i+1)=basis[i];
+		for (i=0;i<NUMNODESP1;i++){
+			*(L+numdof*NUMNODESP1*0+numdof*i)=basis[i]; //L[0][NDOF2*i]=dbasis[0][i];
+			*(L+numdof*NUMNODESP1*0+numdof*i+1)=0;
+			*(L+numdof*NUMNODESP1*1+numdof*i)=0;
+			*(L+numdof*NUMNODESP1*1+numdof*i+1)=basis[i];
 		}
 	}
@@ -337,10 +338,10 @@
 	IssmDouble x1,y1,x2,y2,x3,y3;
 
-	x1=*(xyz_list+NUMNODES*0+0);
-	y1=*(xyz_list+NUMNODES*0+1);
-	x2=*(xyz_list+NUMNODES*1+0);
-	y2=*(xyz_list+NUMNODES*1+1);
-	x3=*(xyz_list+NUMNODES*2+0);
-	y3=*(xyz_list+NUMNODES*2+1);
+	x1=*(xyz_list+NUMNODESP1*0+0);
+	y1=*(xyz_list+NUMNODESP1*0+1);
+	x2=*(xyz_list+NUMNODESP1*1+0);
+	y2=*(xyz_list+NUMNODESP1*1+1);
+	x3=*(xyz_list+NUMNODESP1*2+0);
+	y3=*(xyz_list+NUMNODESP1*2+1);
 
 	*(J+NDOF2*0+0)=0.5*(x2-x1);
@@ -405,4 +406,23 @@
 }
 /*}}}*/
+/*FUNCTION TriaRef::GetNodalFunctionsP2{{{*/
+void TriaRef::GetNodalFunctionsP2(IssmDouble** pbasis,GaussTria* gauss){
+	/*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+	/*Allocate*/
+	IssmDouble* basis =xNew<IssmDouble>(NUMNODESP2);
+
+	basis[0]=gauss->coord1*(2.*gauss->coord1-1.);
+	basis[1]=gauss->coord2*(2.*gauss->coord2-1.);
+	basis[2]=gauss->coord3*(2.*gauss->coord3-1.);
+	basis[3]=4.*gauss->coord3*gauss->coord2;
+	basis[4]=4.*gauss->coord3*gauss->coord1;
+	basis[5]=4.*gauss->coord1*gauss->coord2;
+
+	/*Assign output pointer*/
+	*pbasis = basis;
+
+}
+/*}}}*/
 /*FUNCTION TriaRef::GetSegmentNodalFunctions{{{*/
 void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss,int index1,int index2){
@@ -425,5 +445,5 @@
 	 * actual coordinate system): */
 	int       i;
-	IssmDouble    dbasis_ref[NDOF2][NUMNODES];
+	IssmDouble    dbasis_ref[NDOF2][NUMNODESP1];
 	IssmDouble    Jinv[NDOF2][NDOF2];
 
@@ -439,7 +459,7 @@
 	 * [dhi/dy]       [dhi/ds]
 	 */
-	for (i=0;i<NUMNODES;i++){
-		dbasis[NUMNODES*0+i]=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i];
-		dbasis[NUMNODES*1+i]=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i];
+	for (i=0;i<NUMNODESP1;i++){
+		dbasis[NUMNODESP1*0+i]=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i];
+		dbasis[NUMNODESP1*1+i]=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i];
 	}
 
@@ -452,15 +472,46 @@
 
 	/*First nodal function: */
-	*(dl1dl3+NUMNODES*0+0)=-0.5; 
-	*(dl1dl3+NUMNODES*1+0)=-1.0/(2.0*SQRT3);
+	*(dl1dl3+NUMNODESP1*0+0)=-0.5; 
+	*(dl1dl3+NUMNODESP1*1+0)=-1.0/(2.0*SQRT3);
 
 	/*Second nodal function: */
-	*(dl1dl3+NUMNODES*0+1)=0.5;
-	*(dl1dl3+NUMNODES*1+1)=-1.0/(2.0*SQRT3);
+	*(dl1dl3+NUMNODESP1*0+1)=0.5;
+	*(dl1dl3+NUMNODESP1*1+1)=-1.0/(2.0*SQRT3);
 
 	/*Third nodal function: */
-	*(dl1dl3+NUMNODES*0+2)=0;
-	*(dl1dl3+NUMNODES*1+2)=1.0/SQRT3;
-
+	*(dl1dl3+NUMNODESP1*0+2)=0;
+	*(dl1dl3+NUMNODESP1*1+2)=1.0/SQRT3;
+
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetNodalFunctionsDerivativesP2Reference{{{*/
+void TriaRef::GetNodalFunctionsDerivativesP2Reference(IssmDouble** pdbasis,GaussTria* gauss){
+	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
+	 * natural coordinate system) at the gaussian point. */
+
+	/*Allocate*/
+	IssmDouble* dbasis =xNew<IssmDouble>(NUMNODESP2*2);
+
+	/*Nodal function 1*/
+	dbasis[NUMNODESP2*0+0]=-2.*gauss->coord1 + 0.5;
+	dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.;
+	/*Nodal function 2*/
+	dbasis[NUMNODESP2*0+0]=+2.*gauss->coord2 + 0.5;
+	dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.;
+	/*Nodal function 3*/
+	dbasis[NUMNODESP2*0+0]=0.;
+	dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.;
+	/*Nodal function 4*/
+	dbasis[NUMNODESP2*0+0]=+2.*gauss->coord3;
+	dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3;
+	/*Nodal function 5*/
+	dbasis[NUMNODESP2*0+0]=-2.*gauss->coord3;
+	dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3;
+	/*Nodal function 6*/
+	dbasis[NUMNODESP2*0+0]=2.*(gauss->coord1-gauss->coord2);
+	dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*(gauss->coord1+gauss->coord2);
+
+	/*Assign output pointer*/
+	*pdbasis = dbasis;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/objects/Elements/TriaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/TriaRef.h	(revision 14766)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/TriaRef.h	(revision 14767)
@@ -35,10 +35,12 @@
 		void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss);
 		void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,GaussTria* gauss);
-		void GetNodalFunctions(IssmDouble* l1l2l3,GaussTria* gauss);
-		void GetSegmentNodalFunctions(IssmDouble* l1l2l3,GaussTria* gauss, int index1,int index2);
+		void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss);
+		void GetNodalFunctionsP2(IssmDouble** pbasis,GaussTria* gauss);
+		void GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss, int index1,int index2);
 		void GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2);
 		void GetSegmentBprimeFlux(IssmDouble* Bprime,GaussTria* gauss, int index1,int index2);
-		void GetNodalFunctionsDerivatives(IssmDouble* l1l2l3,IssmDouble* xyz_list, GaussTria* gauss);
+		void GetNodalFunctionsDerivatives(IssmDouble* basis,IssmDouble* xyz_list, GaussTria* gauss);
 		void GetNodalFunctionsDerivativesReference(IssmDouble* dl1dl3,GaussTria* gauss);
+		void GetNodalFunctionsDerivativesP2Reference(IssmDouble** pdbasis,GaussTria* gauss);
 		void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss);
 		void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss);