import numpy
from ElementsFromEdge import *
from MatlabFuncs import *

def meshprocessoutsiderifts(md,domainoutline):
	"""
	MESHPROCESSOUTSIDERIFTS - process rifts when they touch the domain outline

	   Usage:
	      md=meshprocessoutsiderifts(md,domain)

	"""

	#go through rifts, and figure out which ones touch the domain outline
	for rift in md.rifts.riftstruct:
	
		#first, flag nodes that belong to the domain outline
		flags=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),domainoutline,'node',0)

		tips=rift.tips
		outsidetips=tips[numpy.nonzero(flags[rift.tips-1])[0]]

		#we have found outsidetips, tips that touch the domain outline. go through them
		for tip in outsidetips:
		
			#find tip in the segments, take first segment (there should be 2) that holds tip, 
			#and node_connected_to_tip is the other node on this segment:
			tipindex=numpy.nonzero(rift.segments[:,0]==tip)[0]
			if tipindex:
				tipindex=tipindex[0]
				node_connected_to_tip=rift.segments[tipindex,1]
			else:
				tipindex=numpy.nonzero(rift.segments[:,1]==tip)[0]
				tipindex=tipindex[0]
				node_connected_to_tip=rift.segments[tipindex,1]

			#ok, we have the tip node, and the first node connected to it, on the rift. Now, 
			#identify all the elements that are connected to the tip, and that are on the same 
			#side of the rift.
			A=tip
			B=node_connected_to_tip

			elements=numpy.empty(0)

			while flags(B):    #as long as B does not belong to the domain outline, keep looking.
				#detect elements on edge A,B:
				edgeelements=ElementsFromEdge(md.mesh.elements,A,B)
				#rule out those we already detected
				already_detected=ismember(edgeelements,elements)
				nextelement=edgeelements(numpy.nonzero(numpy.logical_not(already_detected))[0])
				#add new detected element to the list of elements we are looking for.
				elements=numpy.concatenate((elements,nextelement))
				#new B:
				B=md.mesh.elements[nextelement-1,numpy.nonzero(numpy.logical_not(ismember(md.mesh.elements[nextelement-1,:],numpy.array([A,B]))))]
		
			#take the list of elements on one side of the rift that connect to the tip, 
			#and duplicate the tip on them, so as to open the rift to the outside.
			num=numpy.size(md.mesh.x)+1
			md.mesh.x=numpy.concatenate((md.mesh.x,md.mesh.x[tip]))
			md.mesh.y=numpy.concatenate((md.mesh.y,md.mesh.y[tip]))
			md.mesh.numberofvertices=num
		
			#replace tip in elements
			newelements=md.mesh.elements[elements.astype(int)-1,:]
			pos=numpy.nonzero(newelements==tip)
			newelements[pos]=num
			md.mesh.elements[elements.astype(int)-1,:]=newelements
			rift.tips=numpy.concatenate((rift.tips,num))

			#deal with segments
			tipsegments=numpy.nonzero(numpy.logical_or(md.mesh.segments[:,0]==tip,md.mesh.segments[:,1]==tip))[0]
			for segment_index in tipsegments:
				pos=numpy.nonzero(md.mesh.segments[segment_index,0:2]!=tip)[0]
				other_node=md.mesh.segments[segment_index,pos]
				if not isconnected(md.mesh.elements,other_node,tip):
					pos=numpy.nonzero(md.mesh.segments[segment_index,0:2]==tip)[0]
					md.mesh.segments[segment_index,pos]=num

	#Fill in rest of fields:
	md.mesh.numberofelements=numpy.size(md.mesh.elements,axis=0)
	md.mesh.numberofvertices=numpy.size(md.mesh.x)
	md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
	md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x))
	md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1
	md.rifts.numrifts=length(md.rifts.riftstruct)
	md.flowequation.element_equation=3.*numpy.ones(md.mesh.numberofelements)
	md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices)
	md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices)
	md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements)
	md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements)

	return md

def isconnected(elements,A,B):    # {{{
	"""
	ISCONNECTED: are two nodes connected by a triangulation?

	   Usage: flag=isconnected(elements,A,B)

	"""

	elements=ElementsFromEdge(elements,A,B)
	if not elements:
		flag=0
	else:
		flag=1

	return flag
	# }}}

