// Gmsh project created on Mon Mar 23 07:38:59 2009

//the 4 points should exist in the .geo file
Point(1) = {0, 0,0};
Point(2) = {4, 0,0};
Point(3) = {1, 1, 1};
Point(4) = {4, 1, 1};
//the two lines should exist in the .geo file
Line(1) = {1, 2};
Line(2) = {4, 3};
// there seems to be no way to retrieve directly yhe points of a given straight line

Function intersectTwoLines
p11[]=Point{1};
p12[]=Point{2};
p21[]=Point{3};
p22[]=Point{4};
//in plane xy y=ax+b
l1lx=p12[0]-p11[0];
If (l1lx==0)
	l1lx=0.00001;
EndIf
l1ly=p12[1]-p11[1];
a1=l1ly/l1lx;
b1=p11[1]-a1*p11[0];
l2lx=p22[0]-p21[0];
If (l2lx==0)
	l2lx=0.00001;
EndIf
l2ly=p22[1]-p21[1];
a2=l2ly/l2lx;
b2=p21[1]-a2*p21[0];
If (a1==a2)
	Printf ("the two lines are parallel, no intesection") /*info in console*/;
EndIf
If (a1!=a2)
	itx=(b2-b1)/(a1-a2);
	ity=a1*itx+b1;
//equations in plane xz z=ax+b
	l1lz=p12[2]-p11[2];
	a1=l1lz/l1lx;
	b1=p11[2]-a1*p11[0];
	l2lz=p22[2]-p21[2];
	a2=l2lz/l2lx;
	b2=p21[2]-a2*p21[0];
//unless the two line lies in the same plane
//there two "intersections"
//one on each line!!
//intestion maye be beyond the limits of the segment
//on the extended infinite lineTwoLines
	itz1=a1*itx+b1;
	itz2=a2*itx+b2;
	p1=newp ; Point(p1)={itx,ity,itz1};
	If (itz2!=itz1)  
		p2=newp ; 
		Point(p2)={itx,ity,itz2};
	EndIf	
EndIf
Return

Call intersectTwoLines;