Index: issm/trunk-jpl/src/m/mesh/bamg.py
===================================================================
--- issm/trunk-jpl/src/m/mesh/bamg.py	(revision 13493)
+++ issm/trunk-jpl/src/m/mesh/bamg.py	(revision 13495)
@@ -7,4 +7,5 @@
 from expread import *
 from expwrite import *
+from SegIntersect import *
 from MatlabFuncs import *
 from BamgMesher import *
@@ -121,91 +122,91 @@
 
 				elif numpy.any(numpy.logical_not(flags)):
-					raise RuntimeError("bamg.py for rifts is not complete.")
 					#We LOTS of work to do
 					print("Rift tip outside of or on the domain has been detected and is being processed...")
-					"""
-
-					%check that only one point is outside (for now)
-					if sum(~flags)~=1,
-						error('bamg error message: only one point outside of the domain is supported yet');
-					end
-
-					%Move tip outside to the first position
-					if flags(1)==0,
-						%OK, first point is outside (do nothing),
-					elseif (flags(end)==0),
-						rift(i).x=flipud(rift(i).x);
-						rift(i).y=flipud(rift(i).y);
-					else
-						error('bamg error message: only a rift tip can be outside of the domain');
-					end
-
-					%Get cordinate of intersection point
-					x1=rift(i).x(1); y1=rift(i).y(1);
-					x2=rift(i).x(2); y2=rift(i).y(2);
-					for j=1:length(domain[0]['x'])-1;
-						if SegIntersect([x1 y1; x2 y2],[domain[0]['x'](j) domain[0]['y'](j); domain[0]['x'](j+1) domain[0]['y'](j+1)]),
-
-							%Get position of the two nodes of the edge in domain
-							i1=j;
-							i2=mod(j+1,domain[0]['nods']);
-
-							%rift is crossing edge [i1 i2] of the domain
-							%Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
-							x3=domain[0]['x'](i1); y3=domain[0]['y'](i1);
-							x4=domain[0]['x'](i2); y4=domain[0]['y'](i2);
-							x=det([det([x1 y1; x2 y2])  x1-x2;det([x3 y3; x4 y4])  x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
-							y=det([det([x1 y1; x2 y2])  y1-y2;det([x3 y3; x4 y4])  y3-y4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
-
-							segdis= sqrt((x4-x3)**2+(y4-y3)**2);
-							tipdis=[sqrt((x-x3)**2+(y-y3)**2)  sqrt((x-x4)**2+(y-y4)**2)];
-
-							if (min(tipdis)/segdis) < getfieldvalue(options,'toltip',0),
-								disp('moving tip-domain intersection point');
-
-								%Get position of the closest point
-								if tipdis(1)>tipdis(2),
-									pos=i2;
-								else
-									pos=i1;
-								end
-
-								%This point is only in Vertices (number pos).
-								%OK, now we can add our own rift
-								nods=rift(i).nods-1;
-								bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(2:end) rift(i).y(2:end) ones(nods,1)]];
-								bamg_geometry.Edges=[bamg_geometry.Edges;...
-									pos count+1  (1+i);...
-									[transpose(count+1:count+nods-1) transpose([count+2:count+nods])  (1+i)*ones(nods-1,1)]];
+
+					#check that only one point is outside (for now)
+					if numpy.sum(numpy.logical_not(flags).astype(int))!=1:
+						raise RuntimeError("bamg error message: only one point outside of the domain is supported yet")
+
+					#Move tip outside to the first position
+					if   not flags[0]:
+						#OK, first point is outside (do nothing),
+						pass
+					elif not flags[-1]:
+						rifti['x']=numpy.flipud(rifti['x'])
+						rifti['y']=numpy.flipud(rifti['y'])
+					else:
+						raise RuntimeError("bamg error message: only a rift tip can be outside of the domain")
+
+					#Get cordinate of intersection point
+					x1=rifti['x'][0]
+					y1=rifti['y'][0]
+					x2=rifti['x'][1]
+					y2=rifti['y'][1]
+					for j in xrange(0,numpy.size(domain[0]['x'])-1):
+						if SegIntersect(numpy.array([[x1,y1],[x2,y2]]),numpy.array([[domain[0]['x'][j],domain[0]['y'][j]],[domain[0]['x'][j+1],domain[0]['y'][j+1]]])):
+
+							#Get position of the two nodes of the edge in domain
+							i1=j
+							i2=j+1
+
+							#rift is crossing edge [i1 i2] of the domain
+							#Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
+							x3=domain[0]['x'][i1]
+							y3=domain[0]['y'][i1]
+							x4=domain[0]['x'][i2]
+							y4=domain[0]['y'][i2]
+#							x=det([det([x1 y1; x2 y2])  x1-x2;det([x3 y3; x4 y4])  x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
+#							y=det([det([x1 y1; x2 y2])  y1-y2;det([x3 y3; x4 y4])  y3-y4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
+							x=numpy.linalg.det(numpy.array([[numpy.linalg.det(numpy.array([[x1,y1],[x2,y2]])),x1-x2],[numpy.linalg.det(numpy.array([[x3,y3],[x4,y4]])),x3-x4]]))/numpy.linalg.det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]]));
+							y=numpy.linalg.det(numpy.array([[numpy.linalg.det(numpy.array([[x1,y1],[x2,y2]])),y1-y2],[numpy.linalg.det(numpy.array([[x3,y3],[x4,y4]])),y3-y4]]))/numpy.linalg.det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]]));
+
+							segdis= sqrt((x4-x3)**2+(y4-y3)**2)
+							tipdis=numpy.array([sqrt((x-x3)**2+(y-y3)**2),sqrt((x-x4)**2+(y-y4)**2)])
+
+							if numpy.min(tipdis)/segdis < options.getfieldvalue('toltip',0):
+								print("moving tip-domain intersection point")
+
+								#Get position of the closer point
+								if tipdis[0]>tipdis[1]:
+									pos=i2
+								else:
+									pos=i1
+
+								#This point is only in Vertices (number pos).
+								#OK, now we can add our own rift
+								nods=rifti['nods']-1
+								bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.hstack((rifti['x'][1:].reshape(-1,1),rifti['y'][1:].reshape(-1,1),numpy.ones((nods,1))))))
+								bamg_geometry.Edges=numpy.vstack((bamg_geometry.Edges,\
+									numpy.array([[pos,count+1,(1+i)]]),\
+									numpy.hstack((numpy.arange(count+1,count+nods).reshape(-1,1),numpy.arange(count+2,count+nods+1).reshape(-1,1),(1+i)*numpy.ones((nods-1,1))))))
+								count+=nods
+
+								break
+
+							else:
+								#Add intersection point to Vertices
+								bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.array([[x,y,1]])))
+								count+=1
+
+								#Decompose the crossing edge into 2 subedges
+								pos=numpy.nonzero(numpy.logical_and(bamg_geometry.Edges[:,0]==i1,bamg_geometry.Edges[:,1]==i2))[0]
+								if not pos:
+									raise RuntimeError("bamg error message: a problem occurred...")
+								bamg_geometry.Edges=numpy.vstack((bamg_geometry.Edges[0:pos-1,:],\
+									numpy.array([[bamg_geometry.Edges[pos,0],count                     ,bamg_geometry.Edges[pos,2]]]),\
+									numpy.array([[count                     ,bamg_geometry.Edges[pos,1],bamg_geometry.Edges[pos,2]]]),\
+									bamg_geometry.Edges[pos+1:,:]))
+
+								#OK, now we can add our own rift
+								nods=rifti['nods']-1
+								bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.hstack((rifti['x'][1:].reshape(-1,1),rifti['y'][1:].reshape(-1,1),numpy.ones((nods,1))))))
+								bamg_geometry.Edges=numpy.vstack((bamg_geometry.Edges,\
+									numpy.array([[count,count+1,2]]),\
+									numpy.hstack((numpy.arange(count+1,count+nods).reshape(-1,1),numpy.arange(count+2,count+nods+1).reshape(-1,1),(1+i)*numpy.ones((nods-1,1))))))
 								count=count+nods;
 
-								break;
-
-							else
-								%Add intersection point to Vertices
-								bamg_geometry.Vertices=[bamg_geometry.Vertices; x y 1];
-								count=count+1;
-
-								%Decompose the crossing edge in 2 subedges
-								pos=find(bamg_geometry.Edges(:,1)==i1 & bamg_geometry.Edges(:,2)==i2);
-								if isempty(pos) error('bamg error message: a problem occured...'); end
-								bamg_geometry.Edges=[bamg_geometry.Edges(1:pos-1,:);...
-									bamg_geometry.Edges(pos,1) count                           bamg_geometry.Edges(pos,3);...
-									count                      bamg_geometry.Edges(pos,2)   bamg_geometry.Edges(pos,3);...
-									bamg_geometry.Edges(pos+1:end,:)];
-
-								%OK, now we can add our own rift
-								nods=rift(i).nods-1;
-								bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(2:end) rift(i).y(2:end) ones(nods,1)]];
-								bamg_geometry.Edges=[bamg_geometry.Edges;...
-									count  count+1  2 ;...
-									[transpose(count+1:count+nods-1) transpose([count+2:count+nods])  (1+i)*ones(nods-1,1)]];
-								count=count+nods;
-
-								break;
-							end
-						end
-					end
-					"""
+								break
+
 				else:
 					nods=rifti['nods']-1
@@ -384,6 +385,6 @@
 
 				#Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
-				x=det(numpy.array([det(numpy.array([[x1,y1],[x2,y2]])),x1-x2],[det(numpy.array([[x3,y3],[x4,y4]])),x3-x4])/det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])))
-				y=det(numpy.array([det(numpy.array([[x1,y1],[x2,y2]])),y1-y2],[det(numpy.array([[x3,y3],[x4,y4]])),y3-y4])/det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])))
+				x=numpy.linalg.det(numpy.array([numpy.linalg.det(numpy.array([[x1,y1],[x2,y2]])),x1-x2],[numpy.linalg.det(numpy.array([[x3,y3],[x4,y4]])),x3-x4])/numpy.linalg.det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])))
+				y=numpy.linalg.det(numpy.array([numpy.linalg.det(numpy.array([[x1,y1],[x2,y2]])),y1-y2],[numpy.linalg.det(numpy.array([[x3,y3],[x4,y4]])),y3-y4])/numpy.linalg.det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])))
 
 				#Add vertex to the list of vertices
