Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18957)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18958)
@@ -1275,5 +1275,5 @@
 			}
 			break;
-		case 2: // two vertices have ice: get area of quadrangle
+		case 2: // two vertices have ice: fraction is computed from first ice vertex to last in CCW fashion
 			for(i=0;i<2;i++){
 				fraction[i]=(level-lsf[indices[i]])/(lsf[indices[numiceverts]]-lsf[indices[i]]);
@@ -3197,95 +3197,51 @@
 /*}}}*/
 void       Tria::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
-
-	int         normal_orientation=0;
-	IssmDouble  s1,s2;
-	IssmDouble  levelset[NUMVERTICES];
-
-	/*Recover parameters and values*/
-	IssmDouble* xyz_zero = xNew<IssmDouble>(2*3);
-	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=1; //orientation of quadrangle depending on distribution of levelsetfunction
-		/*New point 1*/
-		xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]);
-		xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]);
-		xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]);
-
-		/*New point 0*/
-		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]);
-		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]);
-		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+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=1;
-		/*New point 1*/
-		xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]);
-		xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]);
-		xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]);
-
-		/*New point 2*/
-		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]);
-		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]);
-		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+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=1;
-		/*New point 0*/
-		xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]);
-		xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]);
-		xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]);
-
-		/*New point 2*/
-		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]);
-		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]);
-		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]);
-	}
-	else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
-		xyz_zero[3*0+0]=xyz_list[0*3+0];
-		xyz_zero[3*0+1]=xyz_list[0*3+1];
-		xyz_zero[3*0+2]=xyz_list[0*3+2];
-
-		/*New point 2*/
-		xyz_zero[3*1+0]=xyz_list[1*3+0];
-		xyz_zero[3*1+1]=xyz_list[1*3+1];
-		xyz_zero[3*1+2]=xyz_list[1*3+2];
-	}
-	else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
-		xyz_zero[3*0+0]=xyz_list[2*3+0];
-		xyz_zero[3*0+1]=xyz_list[2*3+1];
-		xyz_zero[3*0+2]=xyz_list[2*3+2];
-
-		/*New point 2*/
-		xyz_zero[3*1+0]=xyz_list[0*3+0];
-		xyz_zero[3*1+1]=xyz_list[0*3+1];
-		xyz_zero[3*1+2]=xyz_list[0*3+2];
-	}
-	else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
-		xyz_zero[3*0+0]=xyz_list[1*3+0];
-		xyz_zero[3*0+1]=xyz_list[1*3+1];
-		xyz_zero[3*0+2]=xyz_list[1*3+2];
-
-		/*New point 2*/
-		xyz_zero[3*1+0]=xyz_list[2*3+0];
-		xyz_zero[3*1+1]=xyz_list[2*3+1];
-		xyz_zero[3*1+2]=xyz_list[2*3+2];
-	}
-	else _error_("Case not covered");
-
-	/*Assign output pointer*/
-	*pxyz_zero= xyz_zero;
+	/* Return coordinates where levelset intersects element edges.
+	 * Attention: In case that no intersection exists, NULL pointer is returned.*/
+
+	/*Intermediaries*/
+	const int dim=3;
+	int numiceverts;
+	int i, n, e, counter;
+	IssmDouble s[2];
+	int* indices=NULL;
+	IssmDouble* xyz_zero=NULL;
+
+	this->GetLevelsetIntersection(&indices, &numiceverts, s, MaskIceLevelsetEnum, 0.);
+	
+	//TODO: check if for 2 iceverts front segment is oriented in CCW way
+	
+	if(numiceverts>0) xyz_zero=xNew<IssmDouble>(2*dim);
+	if((numiceverts>0)&&(numiceverts<NUMVERTICES)){
+		counter=0;
+		for(i=0;i<numiceverts;i++){	// iterate over ice vertices
+			for(n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
+				for(e=0;e<dim;e++){ // spatial direction
+					int ind_ice		=dim*indices[i]+e;
+					int ind_noice	=dim*indices[n]+e;
+					int ind			=dim*counter+e;
+					xyz_zero[ind]=xyz_list[ind_ice]+s[counter]*(xyz_list[ind_noice]-xyz_list[ind_ice]);
+				}
+				counter++;
+			}
+		}
+	}
+	else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge
+		IssmDouble lsf[NUMVERTICES];
+		this->GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+		counter=0;
+		for(i=0;i<NUMVERTICES;i++){
+			if(lsf[indices[i]]==0.){
+				for(e=0;e<dim;e++)	xyz_zero[dim*counter+e]=xyz_list[dim*indices[i]+e];
+				counter++;
+			}
+			if(counter==2) break;
+		}
+	}
+	_assert_(counter==2);
+
+	/*Cleanup & return*/
+	xDelete<int>(indices);
+	*pxyz_zero=xyz_zero;
 }
 /*}}}*/