Index: ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h
===================================================================
--- ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h	(revision 18787)
+++ ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h	(revision 18788)
@@ -12,8 +12,9 @@
 void MeshProfileIntersectionx(double** psegments, int* pnumseg, int* index, double* x, double* y, int nel, int nods,  Contour<IssmPDouble>** contours,int numcontours);
 void MeshSegmentsIntersection(double** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numnodes);
 void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xnodes,double* ynodes,double* xc,double* yc,int numnodes);
-void ElementSegment(DataSet* segments_dataset,int el,double* xnodes,double* ynodes,double* xsegment,double* ysegment);
+void ElementSegment(DataSet* segments_dataset,int el,int contouri, double* xnodes,double* ynodes,double* xsegment,double* ysegment);
 int  SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2);
 bool NodeInElement(double* xnodes, double* ynodes, double x, double y);
+bool IsIdenticalNode(double x1, double y1, double x2, double y2, double tolerance);
 
 #endif /* _MESHPROFILEINTERSECTIONX_H */
Index: ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp	(revision 18787)
+++ ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp	(revision 18788)
@@ -128,10 +128,15 @@
 		xsegment[1]=xc[i+1];
 		ysegment[0]=yc[i];
 		ysegment[1]=yc[i+1];
-		ElementSegment(segments_dataset,el, xnodes,ynodes,xsegment,ysegment);
+		/*if (el==318 && i==9){
+			_printf_("contour: " << i << " " << xsegment[0] << " " << ysegment[0] << " " << xsegment[1] << " " << ysegment[1] 
+				<< " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " << 
+				ynodes[2] << "\n");
+		}*/
+		ElementSegment(segments_dataset,el, i, xnodes,ynodes,xsegment,ysegment);
 	}
 }/*}}}*/
-void ElementSegment(DataSet* segments_dataset,int el,double* xnodes,double* ynodes,double* xsegment,double* ysegment){/*{{{*/
+void ElementSegment(DataSet* segments_dataset,int el, int contouri, double* xnodes,double* ynodes,double* xsegment,double* ysegment){/*{{{*/
 
 	/*We have a tria element (xnodes,ynodes) and a segment (xsegment,ysegment). Find whether they intersect. 
 	 * If they do, create a Segment object with the intersection, and add to segments_dataset dataset: */
@@ -158,11 +163,44 @@
 	xel[0]=xnodes[2];  yel[0]=ynodes[2]; xel[1]=xnodes[0];  yel[1]=ynodes[0];
 	edge3=SegmentIntersect(&gamma1,&gamma2, xel,yel,xsegment,ysegment);
 
-	/*edge can be either IntersectEnum (one iand only one intersection between the edge and the segment), ColinearEnum (edge and segment are collinear) and SeparateEnum (no intersection): */
+	/*edge can be either IntersectEnum (one and only one intersection between the edge and the segment), ColinearEnum (edge and segment are collinear) and SeparateEnum (no intersection): */
+		
+	/*if (el==318 && contouri==9){
+		_printf_(edge1 << " " << edge2 << " " << edge3 << " "  << alpha1 << " " << alpha2 << " " << beta1 << " " << beta2 << " " << gamma1 << " " << gamma2 << " " << xsegment[0] << " "  << xsegment[1] << " " << ysegment[0] << " " << ysegment[1] << " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " << ynodes[2]);
+	
+	_printf_("Bool" << (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum));
+	}*/
 
 	if(    (edge1==IntersectEnum) && (edge2==IntersectEnum) && (edge3==IntersectEnum)   ){
-		/*This case is impossible: */
-		_error_("error: a line cannot go through 3 different vertices!");
+
+		/*This can only be the case if the segment intersected through one vertex, meaning a pair from alpha1, beta1 or gamma1  is 0:*/
+		if (alpha1!=0 && alpha1!=1){
+			/*The vertex opposite edge 1 was intersected:*/
+			xfinal[0]=xsegment[0]+alpha1*(xsegment[1]-xsegment[0]);
+			yfinal[0]=ysegment[0]+alpha1*(ysegment[1]-ysegment[0]);
+			xfinal[1]=xnodes[2];
+			yfinal[1]=ynodes[2];
+		}
+		else if (beta1!=0 && beta1!=1){
+			/*The vertex opposite edge 2 was intersected:*/
+			xfinal[0]=xsegment[0]+beta1*(xsegment[1]-xsegment[0]);
+			yfinal[0]=ysegment[0]+beta1*(ysegment[1]-ysegment[0]);
+			xfinal[1]=xnodes[0];
+			yfinal[1]=ynodes[0];
+		}
+		else if (gamma1!=0 && gamma1!=1){
+			/*The vertex opposite edge 3 was intersected:*/
+			xfinal[0]=xsegment[0]+gamma1*(xsegment[1]-xsegment[0]);
+			yfinal[0]=ysegment[0]+gamma1*(ysegment[1]-ysegment[0]);
+			xfinal[1]=xnodes[1];
+			yfinal[1]=ynodes[1];
+		}
+		segments_dataset->AddObject(new  Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
+
+		
+		/*This case is impossible: not quite! */
+		//_printf_(alpha1 << " " << alpha2 << " " << beta1 << " " << beta2 << " " << gamma1 << " " << gamma2 << " " << xsegment[0] << " "  << xsegment[1] << " " << ysegment[0] << " " << ysegment[1] << " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " << ynodes[2]);
+		/* _error_("error: a line cannot go through 3 different vertices!");*/
 	}
 	else if(    ((edge1==IntersectEnum) && (edge2==IntersectEnum)) || ((edge2==IntersectEnum) && (edge3==IntersectEnum)) || ((edge3==IntersectEnum) && (edge1==IntersectEnum))   ){
 
@@ -186,6 +224,10 @@
 		}
 	}
 	else if(  (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum)   ){
+	
+		/*if (el==318 && contouri==9){
+			_printf_("hello" <<  " NodeInElement 0 " << (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) <<  " NodeInElement 1 " << (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])));
+		}*/
 
 		/*segment intersect only 1 edge. Figure out where the first point in the segment is, inside or outside the element, 
 		 * this will decide the coordinate: */
@@ -195,13 +237,39 @@
 			if(edge2==IntersectEnum){coord2=beta1;}
 			if(edge3==IntersectEnum){coord2=gamma1;}
 		}
-		else{
+		else if (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])){
 			if(edge1==IntersectEnum){coord1=alpha1;}
 			if(edge2==IntersectEnum){coord1=beta1;}
 			if(edge3==IntersectEnum){coord1=gamma1;}
 			coord2=1.0;
 		}
+		else{
+			double tolerance=1e-10;
+			/*Ok, we have an issue here. Probably one of the segments' end is on a vertex, within a certain tolerance!*/
+			if (IsIdenticalNode(xnodes[0],ynodes[0],xsegment[0],ysegment[0],tolerance) ||
+				IsIdenticalNode(xnodes[1],ynodes[1],xsegment[0],ysegment[0],tolerance) ||
+				IsIdenticalNode(xnodes[2],ynodes[2],xsegment[0],ysegment[0],tolerance)){
+				
+				/*ok, segments[0] is common to one of our vertices: */
+				//if (el==318 && contouri==9){ _printf_("ok1" << "\n"); }
+				coord1=0;
+				if(edge1==IntersectEnum){coord2=alpha1;}
+				if(edge2==IntersectEnum){coord2=beta1;}
+				if(edge3==IntersectEnum){coord2=gamma1;}
+			}
+			else if (IsIdenticalNode(xnodes[0],ynodes[0],xsegment[1],ysegment[1],tolerance) ||
+				     IsIdenticalNode(xnodes[1],ynodes[1],xsegment[1],ysegment[1],tolerance) ||
+				     IsIdenticalNode(xnodes[2],ynodes[2],xsegment[1],ysegment[1],tolerance)){
 
+				/*ok, segments[1] is common to one of our vertices: */
+				//if (el==318 && contouri==9){ _printf_("ok2" << "\n"); }
+				if(edge1==IntersectEnum){coord1=alpha1;}
+				if(edge2==IntersectEnum){coord1=beta1;}
+				if(edge3==IntersectEnum){coord1=gamma1;}
+				coord2=1.0;
+			}
+		}
+
 		xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]);
 		xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]);
 		yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]);
@@ -343,3 +411,9 @@
 	*palpha=alpha;*pbeta=beta;
 	return IntersectEnum;
 }/*}}}*/
+bool IsIdenticalNode(double x1, double y1, double x2, double y2, double tolerance){ /*{{{*/
+
+	if (sqrt(pow(x1-x2,2.0) + pow(y1-y2,2))<tolerance)return true;
+	else return false;
+
+}/*}}}*/
