Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15521)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15522)
@@ -700,9 +700,9 @@
 	IssmDouble  xyz_bis[3][3];
 
-	GetJacobianDeterminant(&area_init, &xyz_list[0][0],NULL);
+	area_init=GetArea();
 
 	/*Initialize xyz_list with original xyz_list of triangle coordinates*/
 	for(j=0;j<3;j++){ 
-		for(k=0;k<3;j++){
+		for(k=0;k<3;k++){
 			xyz_bis[j][k]=xyz_list[j][k];
 		}
@@ -710,5 +710,5 @@
 	for(i=0;i<numpoints;i++){
 		for(j=0;j<3;j++){ 
-			for(k=0;k<3;j++){
+			for(k=0;k<3;k++){
 				/*Change appropriate line*/
 				xyz_bis[j][k]=xyz_zero[i][k];
@@ -716,9 +716,9 @@
 
 			/*Compute area fraction*/
-			GetJacobianDeterminant(&area_portion, &xyz_bis[0][0],NULL);
+			area_portion=fabs(xyz_bis[1][0]*xyz_bis[2][1] - xyz_bis[1][1]*xyz_bis[2][0] + xyz_bis[0][0]*xyz_bis[1][1] - xyz_bis[0][1]*xyz_bis[1][0] + xyz_bis[2][0]*xyz_bis[0][1] - xyz_bis[2][1]*xyz_bis[0][0])/2.;
 			*(area_coordinates+3*i+j)=area_portion/area_init;
 
 			/*Reinitialize xyz_list*/
-			for(k=0;k<3;j++){
+			for(k=0;k<3;k++){
 				/*Reinitialize xyz_list with original coordinates*/
 				xyz_bis[j][k]=xyz_list[j][k];
@@ -948,5 +948,5 @@
 		s2=levelset[2]/(levelset[2]-levelset[0]);
 
-		if(levelset[2]>0) normal_orientation=0;
+		if(levelset[2]>0) normal_orientation=1;
 		/*New point 1*/
 		xyz_zero[3*normal_orientation+0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
@@ -964,5 +964,5 @@
 		s2=levelset[0]/(levelset[0]-levelset[1]);
 
-		if(levelset[0]>0) normal_orientation=0;
+		if(levelset[0]>0) normal_orientation=1;
 		/*New point 1*/
 		xyz_zero[3*normal_orientation+0]=xyz_list[0][0]+s1*(xyz_list[2][0]-xyz_list[0][0]);
@@ -980,5 +980,5 @@
 		s2=levelset[1]/(levelset[1]-levelset[2]);
 
-		if(levelset[1]>0) normal_orientation=0;
+		if(levelset[1]>0) normal_orientation=1;
 		/*New point 0*/
 		xyz_zero[3*normal_orientation+0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
@@ -990,6 +990,36 @@
 		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]);
-		}
-
+	}
+	else if(levelset[0]==0 && levelset[1]==0){ //front is on point 0 and 1
+		xyz_zero[3*0+0]=xyz_list[0][0];
+		xyz_zero[3*0+1]=xyz_list[0][1];
+		xyz_zero[3*0+2]=xyz_list[0][2];
+
+		/*New point 2*/
+		xyz_zero[3*1+0]=xyz_list[1][0];
+		xyz_zero[3*1+1]=xyz_list[1][1];
+		xyz_zero[3*1+2]=xyz_list[1][2];
+	}
+	else if(levelset[0]==0 && levelset[2]==0){ //front is on point 0 and 1
+		xyz_zero[3*0+0]=xyz_list[0][0];
+		xyz_zero[3*0+1]=xyz_list[0][1];
+		xyz_zero[3*0+2]=xyz_list[0][2];
+
+		/*New point 2*/
+		xyz_zero[3*1+0]=xyz_list[2][0];
+		xyz_zero[3*1+1]=xyz_list[2][1];
+		xyz_zero[3*1+2]=xyz_list[2][2];
+	}
+	else if(levelset[1]==0 && levelset[2]==0){ //front is on point 0 and 1
+		xyz_zero[3*0+0]=xyz_list[1][0];
+		xyz_zero[3*0+1]=xyz_list[1][1];
+		xyz_zero[3*0+2]=xyz_list[1][2];
+
+		/*New point 2*/
+		xyz_zero[3*1+0]=xyz_list[2][0];
+		xyz_zero[3*1+1]=xyz_list[2][1];
+		xyz_zero[3*1+2]=xyz_list[2][2];
+	}
+	else _error_("Case not covered");
 }
 /*}}}*/
@@ -3124,6 +3154,4 @@
 	bool        isfront;
 
-	return NULL;
-
 	/*Retrieve all inputs and parameters*/
 	GetInputListOnVertices(&ls[0],IcelevelsetEnum);
@@ -3137,4 +3165,5 @@
 	}
 
+	return NULL;
 	/*If no front, return NULL*/
 	if(!isfront) return NULL;
@@ -3196,5 +3225,4 @@
 	delete gauss;
 	return pe;
-
 }
 /*}}}*/