Index: issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15511)
+++ issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15517)
@@ -6949,5 +6949,5 @@
 
 	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussPenta(3,2);
+	gauss=new GaussPenta(5,5);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
@@ -7927,5 +7927,5 @@
 
 	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussPenta(3,2);
+	gauss=new GaussPenta(5,5);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
Index: issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15511)
+++ issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15517)
@@ -692,4 +692,40 @@
 }
 /*}}}*/
+/*FUNCTION Tria::GetAreaCoordinates{{{*/
+void Tria::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints){
+	/*Computeportion of the element that is grounded*/ 
+
+	int         i,j,k;
+	IssmDouble  area_init,area_portion;
+	IssmDouble  xyz_bis[3][3];
+
+	GetJacobianDeterminant(&area_init, &xyz_list[0][0],NULL);
+
+	/*Initialize xyz_list with original xyz_list of triangle coordinates*/
+	for(j=0;j<3;j++){ 
+		for(k=0;k<3;j++){
+			xyz_bis[j][k]=xyz_list[j][k];
+		}
+	}
+	for(i=0;i<numpoints;i++){
+		for(j=0;j<3;j++){ 
+			for(k=0;k<3;j++){
+				/*Change appropriate line*/
+				xyz_bis[j][k]=xyz_zero[i][k];
+			}
+
+			/*Compute area fraction*/
+			GetJacobianDeterminant(&area_portion, &xyz_bis[0][0],NULL);
+			*(area_coordinates+3*i+j)=area_portion/area_init;
+
+			/*Reinitialize xyz_list*/
+			for(k=0;k<3;j++){
+				/*Reinitialize xyz_list with original coordinates*/
+				xyz_bis[j][k]=xyz_list[j][k];
+			}
+		}
+	}
+}
+/*}}}*/
 /*FUNCTION Tria::GetDofList {{{*/
 void  Tria::GetDofList(int** pdoflist, int approximation_enum,int setenum){
@@ -724,4 +760,58 @@
 }
 /*}}}*/
+/*FUNCTION Tria::GetGroundedPart{{{*/
+void Tria::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating){
+	/*Computeportion of the element that is grounded*/ 
+
+	bool               floating=true;
+	int                point;
+	const IssmPDouble  epsilon= 1.e-15;
+	IssmDouble         gl[3];
+	IssmDouble         f1,f2;
+
+	/*Recover parameters and values*/
+	GetInputListOnVertices(&gl[0],GLlevelsetEnum);
+
+	/*Be sure that values are not zero*/
+	if(gl[0]==0) gl[0]=gl[0]+epsilon;
+	if(gl[1]==0) gl[1]=gl[1]+epsilon;
+	if(gl[2]==0) gl[2]=gl[2]+epsilon;
+
+	/*Check that not all nodes are grounded or floating*/
+	if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
+		point=0;
+		f1=1.;
+		f2=1.;
+	}
+	else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
+		point=0;
+		f1=0.;
+		f2=0.;
+	}
+	else{
+		if(gl[0]*gl[1]*gl[2]<0) floating=false;
+
+		if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+			point=2;
+			f1=gl[2]/(gl[2]-gl[0]);
+			f2=gl[2]/(gl[2]-gl[1]);
+		}
+		else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+			point=0;
+			f1=gl[0]/(gl[0]-gl[1]);
+			f2=gl[0]/(gl[0]-gl[2]);
+		}
+		else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+			point=1;
+			f1=gl[1]/(gl[1]-gl[2]);
+			f2=gl[1]/(gl[1]-gl[0]);
+		}
+	}
+	*point1=point;
+	*fraction1=f1;
+	*fraction2=f2;
+	*mainlyfloating=floating;
+}
+/*}}}*/
 /*FUNCTION Tria::GetGroundedPortion{{{*/
 IssmDouble Tria::GetGroundedPortion(IssmDouble* xyz_list){
@@ -826,56 +916,80 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetGroundedPart{{{*/
-void Tria::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating){
+/*FUNCTION Tria::GetSegmentNormal {{{*/
+void Tria:: GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[2][3]){
+
+	/*Build unit outward pointing vector*/
+	IssmDouble vector[2];
+	IssmDouble norm;
+
+	vector[0]=xyz_list[1][0] - xyz_list[0][0];
+	vector[1]=xyz_list[1][1] - xyz_list[0][1];
+
+	norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
+
+	normal[0]= + vector[1]/norm;
+	normal[1]= - vector[0]/norm;
+}
+/*}}}*/
+/*FUNCTION Tria::GetZeroLevelsetCoordinates{{{*/
+void Tria::GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum){
 	/*Computeportion of the element that is grounded*/ 
 
-	bool               floating=true;
-	int                point;
-	const IssmPDouble  epsilon= 1.e-15;
-	IssmDouble         gl[3];
-	IssmDouble         f1,f2;
+	int         normal_orientation;
+	IssmDouble  s1,s2;
+	IssmDouble  levelset[3];
 
 	/*Recover parameters and values*/
-	GetInputListOnVertices(&gl[0],GLlevelsetEnum);
-
-	/*Be sure that values are not zero*/
-	if(gl[0]==0) gl[0]=gl[0]+epsilon;
-	if(gl[1]==0) gl[1]=gl[1]+epsilon;
-	if(gl[2]==0) gl[2]=gl[2]+epsilon;
-
-	/*Check that not all nodes are grounded or floating*/
-	if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
-		point=0;
-		f1=1.;
-		f2=1.;
-	}
-	else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
-		point=0;
-		f1=0.;
-		f2=0.;
-	}
-	else{
-		if(gl[0]*gl[1]*gl[2]<0) floating=false;
-
-		if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
-			point=2;
-			f1=gl[2]/(gl[2]-gl[0]);
-			f2=gl[2]/(gl[2]-gl[1]);
-		}
-		else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
-			point=0;
-			f1=gl[0]/(gl[0]-gl[1]);
-			f2=gl[0]/(gl[0]-gl[2]);
-		}
-		else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
-			point=1;
-			f1=gl[1]/(gl[1]-gl[2]);
-			f2=gl[1]/(gl[1]-gl[0]);
-		}
-	}
-	*point1=point;
-	*fraction1=f1;
-	*fraction2=f2;
-	*mainlyfloating=floating;
+	GetInputListOnVertices(&levelset[0],levelsetenum);
+
+	if(levelset[0]*levelset[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+		/*Portion of the segments*/
+		s1=levelset[2]/(levelset[2]-levelset[1]);
+		s2=levelset[2]/(levelset[2]-levelset[0]);
+
+		if(levelset[2]>0) normal_orientation=0;
+		/*New point 1*/
+		*(xyz_zero+3*normal_orientation+0)=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
+		*(xyz_zero+3*normal_orientation+1)=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]);
+		*(xyz_zero+3*normal_orientation+2)=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]);
+
+		/*New point 0*/
+		*(xyz_zero+3*(1-normal_orientation)+0)=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]);
+		*(xyz_zero+3*(1-normal_orientation)+1)=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]);
+		*(xyz_zero+3*(1-normal_orientation)+2)=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]);
+	}
+	else if(levelset[1]*levelset[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+		/*Portion of the segments*/
+		s1=levelset[0]/(levelset[0]-levelset[2]);
+		s2=levelset[0]/(levelset[0]-levelset[1]);
+
+		if(levelset[0]>0) normal_orientation=0;
+		/*New point 1*/
+		*(xyz_zero+3*normal_orientation+0)=xyz_list[0][0]+s1*(xyz_list[2][0]-xyz_list[0][0]);
+		*(xyz_zero+3*normal_orientation+1)=xyz_list[0][1]+s1*(xyz_list[2][1]-xyz_list[0][1]);
+		*(xyz_zero+3*normal_orientation+2)=xyz_list[0][2]+s1*(xyz_list[2][2]-xyz_list[0][2]);
+
+		/*New point 2*/
+		*(xyz_zero+3*(1-normal_orientation)+0)=xyz_list[0][0]+s2*(xyz_list[1][0]-xyz_list[0][0]);
+		*(xyz_zero+3*(1-normal_orientation)+1)=xyz_list[0][1]+s2*(xyz_list[1][1]-xyz_list[0][1]);
+		*(xyz_zero+3*(1-normal_orientation)+2)=xyz_list[0][2]+s2*(xyz_list[1][2]-xyz_list[0][2]);
+	}
+	else if(levelset[0]*levelset[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+		/*Portion of the segments*/
+		s1=levelset[1]/(levelset[1]-levelset[0]);
+		s2=levelset[1]/(levelset[1]-levelset[2]);
+
+		if(levelset[1]>0) normal_orientation=0;
+		/*New point 0*/
+		*(xyz_zero+3*normal_orientation+0)=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
+		*(xyz_zero+3*normal_orientation+1)=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]);
+		*(xyz_zero+3*normal_orientation+2)=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]);
+
+		/*New point 2*/
+		*(xyz_zero+3*(1-normal_orientation)+0)=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]);
+		*(xyz_zero+3*(1-normal_orientation)+1)=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]);
+		*(xyz_zero+3*(1-normal_orientation)+2)=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]);
+		}
+
 }
 /*}}}*/
@@ -3006,25 +3120,70 @@
 
 	/*Intermediaries */
-	int            i,j;
-	IssmDouble     ls[3];
-	IssmDouble     xyz_list[NUMVERTICES][3];
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF2;
-	Icefront *icefront=NULL;
+	IssmDouble  ls[3];
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	bool        isfront;
 
 	return NULL;
+
 	/*Retrieve all inputs and parameters*/
-	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 	GetInputListOnVertices(&ls[0],IcelevelsetEnum);
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
-	GaussTria*     gauss  = new GaussTria(2);
-	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
-
-	/*Create Ice Front if necessary*/
+
+	/*If the level set is awlays <0, there is no ice front here*/
+	isfront = false;
 	if(ls[0]>0. || ls[1]>0. || ls[2]>0.){
 		if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]==0.)){
-			//icefront=new Icefront("2d",inputs,matpar,MacAyealApproximationEnum,analysis_type);
+			isfront = true;
+		}
+	}
+
+	/*If no front, return NULL*/
+	if(!isfront) return NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int         numnodes = this->NumberofNodes();
+	int         numdof   = numnodes*NDOF2;
+	IssmDouble  rho_ice,rho_water,gravity;
+	IssmDouble  Jdet,thickness,bed,water_pressure,ice_pressure,air_pressure;
+	IssmDouble  surface_under_water,base_under_water,pressure;
+	GaussTria*  gauss;
+	IssmDouble* basis = xNew<IssmDouble>(numnodes);
+	IssmDouble  xyz_list_front[2][3];
+	IssmDouble  area_coordinates[2][3];
+	IssmDouble  normal[2];
+
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
+	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+	Input* bed_input      =inputs->GetInput(BedEnum);       _assert_(bed_input);
+	rho_water=matpar->GetRhoWater();
+	rho_ice  =matpar->GetRhoIce();
+	gravity  =matpar->GetG();
+	GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,IcelevelsetEnum);
+	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
+	GetSegmentNormal(&normal[0],xyz_list_front);
+
+	/*Start looping on Gaussian points*/
+	gauss=new GaussTria(area_coordinates,3);
+
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+		thickness_input->GetInputValue(&thickness,gauss);
+		bed_input->GetInputValue(&bed,gauss);
+
+		surface_under_water=min(0.,thickness+bed); // 0 if the top of the glacier is above water level
+		base_under_water=min(0.,bed);              // 0 if the bottom of the glacier is above water level
+		water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
+		ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
+		air_pressure=0;
+
+		pressure = ice_pressure + water_pressure + air_pressure;
+
+		GetSegmentJacobianDeterminant(&Jdet,&xyz_list_front[0][0],gauss);
+		GetNodalFunctions(basis,gauss);
+
+		for (int i=0;i<numnodes;i++){
+			pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
+			pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
 		}
 	}
@@ -3037,4 +3196,5 @@
 	delete gauss;
 	return pe;
+
 }
 /*}}}*/
Index: issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 15511)
+++ issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 15517)
@@ -195,4 +195,5 @@
 		ElementVector* CreatePVectorSlope(void);
 		IssmDouble     GetArea(void);
+		void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints);
 		int            GetElementType(void);
 		void	         GetDofList(int** pdoflist,int approximation_enum,int setenum);
@@ -202,4 +203,6 @@
 		IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
 		void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
+		void           GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[2][3]);
+		void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum);
 		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
 		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
