Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 4784)
+++ /issm/trunk/src/c/Makefile.am	(revision 4785)
@@ -506,4 +506,6 @@
 					./modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp\
 					./modules/MeshProfileIntersectionx/ElementSegment.cpp\
+					./modules/MeshProfileIntersectionx/SegmentIntersect.cpp\
+					./modules/MeshProfileIntersectionx/NodeInElement.cpp\
 					./modules/ContourToMeshx/ContourToMeshx.cpp\
 					./modules/ContourToMeshx/ContourToMeshxt.cpp\
@@ -1065,4 +1067,6 @@
 					./modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp\
 					./modules/MeshProfileIntersectionx/ElementSegment.cpp\
+					./modules/MeshProfileIntersectionx/SegmentIntersect.cpp\
+					./modules/MeshProfileIntersectionx/NodeInElement.cpp\
 					./modules/ContourToMeshx/ContourToMeshx.cpp\
 					./modules/ContourToMeshx/ContourToMeshxt.cpp\
Index: /issm/trunk/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp
===================================================================
--- /issm/trunk/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp	(revision 4784)
+++ /issm/trunk/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp	(revision 4785)
@@ -4,11 +4,79 @@
 #include "./MeshProfileIntersectionx.h"
 		
-void ElementSegment(DataSet* segments_dataset,double* xgrids,double* ygrids,double* xsegment,double* ysegment){
-
+void ElementSegment(DataSet* segments_dataset,int el,double* xgrids,double* ygrids,double* xsegment,double* ysegment){
 
 	/*We have a tria element (xgrids,ygrids) 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: */
 
-	ISSMERROR("not supported yet!");
+	int i;
+	double alpha;
+	double alpha1,alpha2;
+	double beta1,beta2;
+	double gamma1,gamma2;
+	
+	int    edge1,edge2,edge3;
 
+	double xel[2],yel[2];
+	double coord1,coord2;
+	double xfinal[2],yfinal[2];
+
+	
+	/*edge 1: */
+	xel[0]=xgrids[0];  yel[0]=ygrids[0]; xel[1]=xgrids[1];  yel[1]=ygrids[1];
+	edge1=SegmentIntersect(&alpha1,&alpha2, xel,yel,xsegment,ysegment);
+
+	/*edge 2: */
+	xel[0]=xgrids[1];  yel[0]=ygrids[1]; xel[1]=xgrids[2];  yel[1]=ygrids[2];
+	edge2=SegmentIntersect(&beta1,&beta2, xel,yel,xsegment,ysegment);
+
+	/*edge 3: */
+	xel[0]=xgrids[2];  yel[0]=ygrids[2]; xel[1]=xgrids[0];  yel[1]=ygrids[0];
+	edge3=SegmentIntersect(&gamma1,&gamma2, xel,yel,xsegment,ysegment);
+
+	if(    (edge1==IntersectEnum) && (edge2==IntersectEnum) && (edge3==IntersectEnum)   ){
+		/*This case is impossible: */
+		ISSMERROR(" error: a line cannot go through 3 different vertices!");
+	}
+	else if(    ((edge1==IntersectEnum) && (edge2==IntersectEnum)) || ((edge2==IntersectEnum) && (edge3==IntersectEnum)) || ((edge3==IntersectEnum) && (edge1==IntersectEnum))   ){
+	
+		/*segment interscts 2 opposite edges of our triangle, at 2 segment coordinates, pick up the lowest (coord1) and highest (coord2): */
+		if((edge1==IntersectEnum) && (edge2==IntersectEnum)) {coord1=min(alpha1,beta1); coord2=max(alpha1,beta1);}
+		if((edge2==IntersectEnum) && (edge3==IntersectEnum)) {coord1=min(beta1,gamma1); coord2=max(beta1,gamma1);}
+		if((edge3==IntersectEnum) && (edge1==IntersectEnum)) {coord1=min(gamma1,alpha1); coord2=max(gamma1,alpha1);}
+
+		/*check this segment did not intersect at a vertex of the tria: */
+		if(coord1!=coord2){
+
+			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]);
+			yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
+
+			segments_dataset->AddObject(new  Segment(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
+		}
+	}
+	else if(  (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum)   ){
+
+		/*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: */
+		if (NodeInElement(xgrids,ygrids,xsegment[0],ysegment[0])){
+			coord1=0;
+			if(edge1==IntersectEnum){coord2=alpha1;}
+			if(edge2==IntersectEnum){coord2=beta1;}
+			if(edge3==IntersectEnum){coord2=gamma1;}
+		}
+		else{
+			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]);
+		yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
+
+		segments_dataset->AddObject(new  Segment(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
+	}
 }
Index: /issm/trunk/src/c/modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp
===================================================================
--- /issm/trunk/src/c/modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp	(revision 4784)
+++ /issm/trunk/src/c/modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp	(revision 4785)
@@ -4,5 +4,5 @@
 #include "./MeshProfileIntersectionx.h"
 		
-void ElementSegmentsIntersection(DataSet* segments_dataset,double* xgrids,double* ygrids,double* xc,double* yc,int numgrids){
+void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xgrids,double* ygrids,double* xc,double* yc,int numgrids){
 
 	int i;
@@ -18,5 +18,5 @@
 		ysegment[1]=yc[i+1];
 
-		ElementSegment(segments_dataset,xgrids,ygrids,xsegment,ysegment);
+		ElementSegment(segments_dataset,el, xgrids,ygrids,xsegment,ysegment);
 
 	}
Index: /issm/trunk/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h
===================================================================
--- /issm/trunk/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h	(revision 4784)
+++ /issm/trunk/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h	(revision 4785)
@@ -13,6 +13,8 @@
 void MeshProfileIntersectionx( double** psegments, int* pnumseg, int* index, double* x, double* y, int nel, int nods,  Contour** contours,int numcontours);
 void MeshSegmentsIntersection(double** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numgrids);
-void ElementSegmentsIntersection(DataSet* segments_dataset,double* xgrids,double* ygrids,double* xc,double* yc,int numgrids);
-void ElementSegment(DataSet* segments_dataset,double* xgrids,double* ygrids,double* xsegment,double* ysegment);
+void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xgrids,double* ygrids,double* xc,double* yc,int numgrids);
+void ElementSegment(DataSet* segments_dataset,int el,double* xgrids,double* ygrids,double* xsegment,double* ysegment);
+int  SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2);
+bool NodeInElement(double* xgrids, double* ygrids, double x, double y);
 
 #endif /* _MESHPROFILEINTERSECTIONX_H */
Index: /issm/trunk/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp
===================================================================
--- /issm/trunk/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp	(revision 4784)
+++ /issm/trunk/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp	(revision 4785)
@@ -27,5 +27,5 @@
 			ygrids[j]=y[*(index+3*i+j)];
 		}
-		ElementSegmentsIntersection(segments_dataset,xgrids,ygrids,xc,yc,numgrids);
+		ElementSegmentsIntersection(segments_dataset,i,xgrids,ygrids,xc,yc,numgrids);
 	}
 
Index: /issm/trunk/src/c/modules/MeshProfileIntersectionx/NodeInElement.cpp
===================================================================
--- /issm/trunk/src/c/modules/MeshProfileIntersectionx/NodeInElement.cpp	(revision 4785)
+++ /issm/trunk/src/c/modules/MeshProfileIntersectionx/NodeInElement.cpp	(revision 4785)
@@ -0,0 +1,33 @@
+/*! \file  NodeInElement.cpp
+*/
+
+#include "./MeshProfileIntersectionx.h"
+
+bool NodeInElement(double* xgrids, double* ygrids, double x, double y){
+
+	double x1,y1;
+	double x2,y2;
+	double x3,y3;
+	double lambda1,lambda2,lambda3;
+	double det;
+
+	x1=xgrids[0];
+	x2=xgrids[1];
+	x3=xgrids[2];
+	y1=ygrids[0];
+	y2=ygrids[1];
+	y3=ygrids[2];
+
+
+	/*compute determinant: */
+	det=x1*y2-x1*y3-x3*y2-x2*y1+x2*y3+x3*y1;
+	
+	/*area coordinates: */
+	lambda1=((y2-y3)*(x-x3)+(x3-x2)*(y-y3))/det;
+	lambda2=((y3-y1)*(x-x3)+(x1-x3)*(y-y3))/det;
+	lambda3=1-lambda1-lambda2;
+
+	if( ((lambda1<=1) && (lambda1>=0)) && ((lambda2<=1) && (lambda2>=0)) && ((lambda3<=1) && (lambda3>=0))  )return true;
+	else return false;
+
+}
Index: /issm/trunk/src/c/modules/MeshProfileIntersectionx/SegmentIntersect.cpp
===================================================================
--- /issm/trunk/src/c/modules/MeshProfileIntersectionx/SegmentIntersect.cpp	(revision 4785)
+++ /issm/trunk/src/c/modules/MeshProfileIntersectionx/SegmentIntersect.cpp	(revision 4785)
@@ -0,0 +1,106 @@
+/*! \file  SegmentIntersect.cpp
+*/
+
+#include "./MeshProfileIntersectionx.h"
+
+int SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2){
+
+	/*See ISSM_DIR/src/m/utils/Geometry/SegIntersect.m for matlab routine from which we take this routine: */
+
+	/*output: */
+	double alpha=-1;
+	double beta=-1;
+	
+	double xA,xB,xC,xD,yA,yB,yC,yD;
+	double O2A[2],O2B[2],O1C[2],O1D[2];
+	double n1[2],n2[2];
+	double test1, test2, test3, test4;
+	double det;
+	double O2O1[2];
+	double pO1A,pO1B,pO1C,pO1D;
+
+	xA=x1[0]; yA=y1[0];
+	xB=x1[1]; yB=y1[1];
+	xC=x2[0]; yC=y2[0];
+	xD=x2[1]; yD=y2[1];
+
+	O2A[0]=xA -(xD/2+xC/2); O2A[1]=yA -(yD/2+yC/2);
+	O2B[0]=xB -(xD/2+xC/2); O2B[1]=yB -(yD/2+yC/2);
+	O1C[0]=xC -(xA/2+xB/2); O1C[1]=yC -(yA/2+yB/2);
+	O1D[0]=xD -(xA/2+xB/2); O1D[1]=yD -(yA/2+yB/2);
+
+
+	n1[0]=yA-yB; n1[1]=xB-xA;  //normal vector to segA
+	n2[0]=yC-yD; n2[1]=xD-xC;  //normal vector to segB
+
+	test1=n2[0]*O2A[0]+n2[1]*O2A[1];
+	test2=n2[0]*O2B[0]+n2[1]*O2B[1];
+
+	if (test1*test2>0){
+		return SeparateEnum;
+	}
+
+	test3=n1[0]*O1C[0]+n1[1]*O1C[1];
+	test4=n1[0]*O1D[0]+n1[1]*O1D[1];
+
+	if (test3*test4>0){
+		return SeparateEnum;
+	}
+
+	/*If colinear: */
+	det=n1[0]*n2[1]-n2[0]*y1[1];
+
+	if(test1*test2==0 & test3*test4==0 & det==0){
+
+		//projection on the axis O1O2
+		O2O1[0]=(xA/2+xB/2)-(xD/2+xC/2);
+		O2O1[1]=(yA/2+yB/2)-(yD/2+yC/2);
+
+		pO1A=O2O1[0]*(O2A[0]-O2O1[0])+O2O1[1]*(O2A[1]-O2O1[1]);
+		pO1B=O2O1[0]*(O2B[0]-O2O1[0])+O2O1[1]*(O2B[1]-O2O1[1]);
+		pO1C=O2O1[0]*O1C[0]+O2O1[1]*O1C[1];
+		pO1D=O2O1[0]*O1D[0]+O2O1[1]*O1D[1];
+
+		//test if one point is included in the other segment (->intersects=true)
+		if ((pO1C-pO1A)*(pO1D-pO1A)<0){
+			alpha=0; beta=0;
+			*palpha=alpha;*pbeta=beta;
+			return ColinearEnum;
+		}
+		if ((pO1C-pO1B)*(pO1D-pO1B)<0){
+			alpha=0; beta=0;
+			*palpha=alpha;*pbeta=beta;
+			return ColinearEnum;
+		}
+		if ((pO1A-pO1C)*(pO1B-pO1C)<0){
+			alpha=0; beta=0;
+			*palpha=alpha;*pbeta=beta;
+			return ColinearEnum;
+		}
+		if ((pO1A-pO1D)*(pO1B-pO1D)<0){
+			alpha=0; beta=0;
+			*palpha=alpha;*pbeta=beta;
+			return ColinearEnum;
+		}
+
+		//test if the 2 segments have the same middle (->intersects=true)
+		if (O2O1==0){
+			alpha=0; beta=0;
+			*palpha=alpha;*pbeta=beta;
+			return ColinearEnum;
+		}
+
+		//if we are here, both segments are colinear, but do not interset:
+		alpha=-1; beta=-1;
+		*palpha=alpha;*pbeta=beta;
+		return SeparateEnum;
+	}
+
+	/*if we are here, both segments intersect. Determine where in the segment coordinate 
+	 * system: */
+	beta=-1;
+	alpha=-(xA*yB-xC*yB+yC*xB-yC*xA+xC*yA-yA*xB)/(-xD*yB+xD*yA+xC*yB-xC*yA-yD*xA+yD*xB+yC*xA-yC*xB); //from intersect.m in formal calculus
+
+	*palpha=alpha;*pbeta=beta;
+	return IntersectEnum;
+}
Index: /issm/trunk/src/c/modules/MeshProfileIntersectionx/intersect.m
===================================================================
--- /issm/trunk/src/c/modules/MeshProfileIntersectionx/intersect.m	(revision 4785)
+++ /issm/trunk/src/c/modules/MeshProfileIntersectionx/intersect.m	(revision 4785)
@@ -0,0 +1,12 @@
+syms xA yA xB yB xC yC xD yD alpha beta x y
+
+A=[xA;yA];
+B=[xB;yB];
+C=[xC;yC];
+D=[xD;yD];
+
+
+Eq=C+alpha*(D-C)-A+beta*(B-A); 
+
+%from Eq, we specify the system to solve: 
+S=solve( xC+alpha*(xD-xC)-xA+beta*(xB-xA), yC+alpha*(yD-yC)-yA+beta*(yB-yA),alpha,beta);
Index: /issm/trunk/src/m/enum/ColinearEnum.m
===================================================================
--- /issm/trunk/src/m/enum/ColinearEnum.m	(revision 4785)
+++ /issm/trunk/src/m/enum/ColinearEnum.m	(revision 4785)
@@ -0,0 +1,11 @@
+function macro=ColinearEnum()
+%COLINEARENUM - Enum of Colinear
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/SynchronizeMatlabEnum
+%            Please read src/c/README for more information
+%
+%   Usage:
+%      macro=ColinearEnum()
+
+macro=234;
Index: /issm/trunk/src/m/enum/IntersectEnum.m
===================================================================
--- /issm/trunk/src/m/enum/IntersectEnum.m	(revision 4785)
+++ /issm/trunk/src/m/enum/IntersectEnum.m	(revision 4785)
@@ -0,0 +1,11 @@
+function macro=IntersectEnum()
+%INTERSECTENUM - Enum of Intersect
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/SynchronizeMatlabEnum
+%            Please read src/c/README for more information
+%
+%   Usage:
+%      macro=IntersectEnum()
+
+macro=233;
Index: /issm/trunk/src/m/enum/SeparateEnum.m
===================================================================
--- /issm/trunk/src/m/enum/SeparateEnum.m	(revision 4785)
+++ /issm/trunk/src/m/enum/SeparateEnum.m	(revision 4785)
@@ -0,0 +1,11 @@
+function macro=SeparateEnum()
+%SEPARATEENUM - Enum of Separate
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/SynchronizeMatlabEnum
+%            Please read src/c/README for more information
+%
+%   Usage:
+%      macro=SeparateEnum()
+
+macro=235;
Index: /issm/trunk/src/m/utils/Geometry/SegIntersect.m
===================================================================
--- /issm/trunk/src/m/utils/Geometry/SegIntersect.m	(revision 4784)
+++ /issm/trunk/src/m/utils/Geometry/SegIntersect.m	(revision 4785)
@@ -4,4 +4,5 @@
 %   return 1 if the two segments intersect
 %   seg1=[x1 y1; x2 y2]
+%   seg2=[x1 y1; x2 y2]
 %
 %   Usage: