import os.path
import numpy
from collections import OrderedDict
from pairoptions import *
from bamggeom import *
from bamgmesh import *
from expread import *
from expwrite import *
from SegIntersect import *
from MatlabFuncs import *
from BamgMesher import *
from ContourToNodes import *

def bamg(md,*args):
	"""
	BAMG - mesh generation

	   Available options (for more details see ISSM website http://issm.jpl.nasa.gov/):

	   - domain :            followed by an ARGUS file that prescribes the domain outline
	   - hmin :              minimum edge length (default is 10^-100)
	   - hmax :              maximum edge length (default is 10^100)
	   - hVertices :         imposed edge length for each vertex (geometry or mesh)
	   - hminVertices :      minimum edge length for each vertex (mesh)
	   - hmaxVertices :      maximum edge length for each vertex (mesh)

	   - anisomax :          maximum ratio between the smallest and largest edges (default is 10^30)
	   - coeff :             coefficient applied to the metric (2-> twice as many elements, default is 1)
	   - cutoff :            scalar used to compute the metric when metric type 2 or 3 are applied
	   - err :               error used to generate the metric from a field
	   - errg :              geometric error (default is 0.1)
	   - field :             field of the model that will be used to compute the metric
	                         to apply several fields, use one column per field
	   - gradation :         maximum ratio between two adjacent edges
	   - Hessiantype :       0 -> use double P2 projection (default)
	                         1 -> use Green formula
	   - KeepVertices :      try to keep initial vertices when adaptation is done on an existing mesh (default 1)
	   - MaxCornerAngle :    maximum angle of corners in degree (default is 10)
	   - maxnbv :            maximum number of vertices used to allocate memory (default is 10^6)
	   - maxsubdiv :         maximum subdivision of exisiting elements (default is 10)
	   - metric :            matrix (numberofnodes x 3) used as a metric
	   - Metrictype :        1 -> absolute error          c/(err coeff^2) * Abs(H)        (default)
	                         2 -> relative error          c/(err coeff^2) * Abs(H)/max(s,cutoff*max(s))
	                         3 -> rescaled absolute error c/(err coeff^2) * Abs(H)/(smax-smin)
	   - nbjacoby :          correction used by Hessiantype=1 (default is 1)
	   - nbsmooth :          number of metric smoothing procedure (default is 3)
	   - omega :             relaxation parameter of the smoothing procedure (default is 1.8)
	   - power :             power applied to the metric (default is 1)
	   - splitcorners :      split triangles whuch have 3 vertices on the outline (default is 1)
	   - geometricalmetric : take the geometry into account to generate the metric (default is 0)
	   - verbose :           level of verbosity (default is 1)

	   - rifts :             followed by an ARGUS file that prescribes the rifts
	   - toltip :            tolerance to move tip on an existing point of the domain outline
	   - tracks :            followed by an ARGUS file that prescribes the tracks that the mesh will stick to
	   - RequiredVertices :  mesh vertices that are required. [x,y,ref]; ref is optional
	   - tol :               if the distance between 2 points of the domain outline is less than tol, they
	                         will be merged

	   Examples:
	      md=bamg(md,'domain','DomainOutline.exp','hmax',3000);
	      md=bamg(md,'field',[md.inversion.vel_obs md.geometry.thickness],'hmax',20000,'hmin',1000);
	      md=bamg(md,'metric',A,'hmin',1000,'hmax',20000,'gradation',3,'anisomax',1);
	"""

	#process options
	options=pairoptions(*args)
#	options=deleteduplicates(options,1);

	#initialize the structures required as input of Bamg
	bamg_options=OrderedDict()
	bamg_geometry=bamggeom()
	bamg_mesh=bamgmesh()

	# Bamg Geometry parameters {{{
	if options.exist('domain'):

		#Check that file exists
		domainfile=options.getfieldvalue('domain')
		if not os.path.exists(domainfile):
			raise IOError("bamg error message: file '%s' not found" % domainfile)
		domain=expread(domainfile)

		#Build geometry 
		count=0
		for i,domaini in enumerate(domain):

			#Check that the domain is closed
			if (domaini['x'][0]!=domaini['x'][-1] or domaini['y'][0]!=domaini['y'][-1]):
				raise RuntimeError("bamg error message: all contours provided in ''domain'' should be closed")

			#Checks that all holes are INSIDE the principle domain outline
			if i:
				flags=ContourToNodes(domaini['x'],domaini['y'],domainfile,0)
				if numpy.any(numpy.logical_not(flags)):
					raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain")

			#Add all points to bamg_geometry
			nods=domaini['nods']-1    #the domain are closed 0=end
			bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.hstack((domaini['x'][0:nods].reshape(-1,1),domaini['y'][0:nods].reshape(-1,1),numpy.ones((nods,1))))))
			bamg_geometry.Edges   =numpy.vstack((bamg_geometry.Edges,   numpy.hstack((numpy.arange(count+1,count+nods+1).reshape(-1,1),numpy.hstack((numpy.arange(count+2,count+nods+1),count+1)).reshape(-1,1),1.*numpy.ones((nods,1))))))
			if i:
				bamg_geometry.SubDomains=numpy.vstack((bamg_geometry.SubDomains,[2,count+1,1,1]))

			#update counter
			count+=nods

		#take care of rifts
		if options.exist('rifts'):

			#Check that file exists
			riftfile=options.getfieldvalue('rifts')
			if not os.path.exists(riftfile):
				raise IOError("bamg error message: file '%s' not found" % riftfile)
			rift=expread(riftfile)

			for i,rifti in enumerate(rift):

				#detect whether all points of the rift are inside the domain
				flags=ContourToNodes(rifti['x'],rifti['y'],domain[0],0)
				if numpy.all(numpy.logical_not(flags)):
					raise RuntimeError("one rift has all its points outside of the domain outline")

				elif numpy.any(numpy.logical_not(flags)):
					#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 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+=nods

								break

				else:
					nods=rifti['nods']-1
					bamg_geometry.Vertices=numpy.vstack(bamg_geometry.Vertices, numpy.hstack(rifti['x'][:],rifti['y'][:],numpy.ones((nods+1,1))))
					bamg_geometry.Edges   =numpy.vstack(bamg_geometry.Edges, numpy.hstack(numpy.arange(count+1,count+nods).reshape(-1,1),numpy.arange(count+2,count+nods+1).reshape(-1,1),i*numpy.ones((nods,1))))
					count=+nods+1

		#Deal with tracks
		if options.exist('tracks'):

			#read tracks
			track=options.getfieldvalue('tracks')
			if all(isinstance(track,(str,unicode))):
				A=expread(track)
				track=numpy.hstack((A.x.reshape(-1,1),A.y.reshape(-1,1)))
			else:
				track=float(track)    #for some reason, it is of class "single"
			if numpy.size(track,axis=1)==2:
				track=numpy.hstack((track,3.*numpy.ones((size(track,axis=0),1))))

			#only keep those inside
			flags=ContourToNodes(track[:,0],track[:,1],domainfile,0)
			track=track[numpy.nonzero(flags),:]

			#Add all points to bamg_geometry
			nods=numpy.size(track,axis=0)
			bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,track))
			bamg_geometry.Edges   =numpy.vstack((bamg_geometry.Edges,numpy.hstack((numpy.arange(count+1,count+nods).reshape(-1,1),numpy.arange(count+2,count+nods+1).reshape(-1,1),3.*numpy.ones((nods-1,1))))))

			#update counter
			count+=nods

		#Deal with vertices that need to be kept by mesher
		if options.exist('RequiredVertices'):

			#recover RequiredVertices
			requiredvertices=float(options.getfieldvalue('RequiredVertices'))    #for some reason, it is of class "single"
			if numpy.size(requiredvertices,axis=1)==2:
				requiredvertices=numpy.hstack((requiredvertices,4.*numpy.ones((numpy.size(requiredvertices,axis=0),1))))
	
			#only keep those inside
			flags=ContourToNodes(requiredvertices[:,0],domain[0],0)
			requiredvertices=requiredvertices[numpy.nonzero(flags),:]

			#Add all points to bamg_geometry
			nods=numpy.size(requiredvertices,axis=0)
			bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,requiredvertices))

			#update counter
			count+=nods

		#process geom
		#bamg_geometry=processgeometry(bamg_geometry,options.getfieldvalue('tol',float(nan)),domain[0])

	elif isinstance(md.private.bamg,dict) and 'geometry' in md.private.bamg:
		bamg_geometry=bamggeom(md.private.bamg['geometry'].__dict__) 
	else:
		#do nothing...
		pass
	#}}}
	# Bamg Mesh parameters {{{
	if not options.exist('domain') and md.mesh.numberofvertices and md.mesh.dimension==2:

		if isinstance(md.private.bamg,dict) and 'mesh' in md.private.bamg:
			bamg_mesh=bamgmesh(md.private.bamg['mesh'].__dict__)
		else:
			bamg_mesh.Vertices=numpy.hstack((md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),numpy.ones((md.mesh.numberofvertices,1))))
			bamg_mesh.Triangles=numpy.hstack((md.mesh.elements,numpy.ones((md.mesh.numberofelements,1))))

		if isinstance(md.rifts.riftstruct,dict):
			raise TypeError("bamg error message: rifts not supported yet. Do meshprocessrift AFTER bamg")
	#}}}
	# Bamg Options {{{
	bamg_options['Crack']=options.getfieldvalue('Crack',0)
	bamg_options['anisomax']=options.getfieldvalue('anisomax',10.**30)
	bamg_options['coeff']=options.getfieldvalue('coeff',1.)
	bamg_options['cutoff']=options.getfieldvalue('cutoff',10.**-5)
	bamg_options['err']=options.getfieldvalue('err',numpy.array([[0.01]]))
	bamg_options['errg']=options.getfieldvalue('errg',0.1)
	bamg_options['field']=options.getfieldvalue('field',numpy.empty((0,1)))
	bamg_options['gradation']=options.getfieldvalue('gradation',1.5)
	bamg_options['Hessiantype']=options.getfieldvalue('Hessiantype',0)
	bamg_options['hmin']=options.getfieldvalue('hmin',10.**-100)
	bamg_options['hmax']=options.getfieldvalue('hmax',10.**100)
	bamg_options['hminVertices']=options.getfieldvalue('hminVertices',numpy.empty((0,1)))
	bamg_options['hmaxVertices']=options.getfieldvalue('hmaxVertices',numpy.empty((0,1)))
	bamg_options['hVertices']=options.getfieldvalue('hVertices',numpy.empty((0,1)))
	bamg_options['KeepVertices']=options.getfieldvalue('KeepVertices',1)
	bamg_options['MaxCornerAngle']=options.getfieldvalue('MaxCornerAngle',10.)
	bamg_options['maxnbv']=options.getfieldvalue('maxnbv',10**6)
	bamg_options['maxsubdiv']=options.getfieldvalue('maxsubdiv',10.)
	bamg_options['metric']=options.getfieldvalue('metric',numpy.empty((0,1)))
	bamg_options['Metrictype']=options.getfieldvalue('Metrictype',0)
	bamg_options['nbjacobi']=options.getfieldvalue('nbjacobi',1)
	bamg_options['nbsmooth']=options.getfieldvalue('nbsmooth',3)
	bamg_options['omega']=options.getfieldvalue('omega',1.8)
	bamg_options['power']=options.getfieldvalue('power',1.)
	bamg_options['splitcorners']=options.getfieldvalue('splitcorners',1)
	bamg_options['geometricalmetric']=options.getfieldvalue('geometricalmetric',0)
	bamg_options['verbose']=options.getfieldvalue('verbose',1)
	#}}}

	#call Bamg
	[bamgmesh_out,bamggeom_out]=BamgMesher(bamg_mesh.__dict__,bamg_geometry.__dict__,bamg_options)

	# plug results onto model
	md.private.bamg=OrderedDict()
	md.private.bamg['mesh']=bamgmesh(bamgmesh_out)
	md.private.bamg['geometry']=bamggeom(bamggeom_out)
	md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
	md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
	md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
	md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
	md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
	md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)

	#Fill in rest of fields:
	md.mesh.dimension=2
	md.mesh.numberofelements=numpy.size(md.mesh.elements,axis=0)
	md.mesh.numberofvertices=numpy.size(md.mesh.x)
	md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0)
	md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
	md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
	md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices,bool)
	md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
	md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
	md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
	md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices,bool)
	md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
	md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity
	md.mesh.elementconnectivity[numpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=0
	md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int)

	#Check for orphan
	if numpy.any(numpy.logical_not(numpy.in1d(numpy.arange(1,md.mesh.numberofvertices+1),md.mesh.elements.flat))):
		raise RuntimeError("Output mesh has orphans. Decrease MaxCornerAngle to prevent outside points (ex: 0.01)")

	return md

def processgeometry(geom,tol,outline):    # {{{

	raise RuntimeError("bamg.py/processgeometry is not complete.")
	#Deal with edges
	print "Checking Edge crossing..."
	i=0
	while (i<numpy.size(geom.Edges,axis=0)):

		#edge counter
		i+=1

		#Get coordinates
		x1=geom.Vertices[geom.Edges[i,0],0]
		y1=geom.Vertices[geom.Edges[i,0],1]
		x2=geom.Vertices[geom.Edges[i,1],0]
		y2=geom.Vertices[geom.Edges[i,1],1]
		color1=geom.Edges[i,2]

		j=i    #test edges located AFTER i only
		while (j<numpy.size(geom.Edges,axis=0)):

			#edge counter
			j+=1

			#Skip if the two edges already have a vertex in common
			if any(ismember(geom.Edges[i,0:2],geom.Edges[j,0:2])):
				continue

			#Get coordinates
			x3=geom.Vertices[geom.Edges[j,0],0]
			y3=geom.Vertices[geom.Edges[j,0],1]
			x4=geom.Vertices[geom.Edges[j,1],0]
			y4=geom.Vertices[geom.Edges[j,1],1]
			color2=geom.Edges[j,2]

			#Check if the two edges are crossing one another
			if SegIntersect(numpy.array([[x1,y1],[x2,y2]]),numpy.array([[x3,y3],[x4,y4]])):

				#Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
				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
				geom.Vertices=numpy.vstack((geom.Vertices,[x,y,min(color1,color2)]))
				id=numpy.size(geom.Vertices,axis=0)

				#Update edges i and j
				edgei=geom.Edges[i,:].copy()
				edgej=geom.Edges[j,:].copy()
				geom.Edges[i,:]    =[edgei(0),id      ,edgei(2)]
				geom.Edges=numpy.vstack((geom.Edges,[id      ,edgei(1),edgei(2)]))
				geom.Edges[j,:]    =[edgej(0),id      ,edgej(2)]
				geom.Edges=numpy.vstack((geom.Edges,[id      ,edgej(1),edgej(2)]))

				#update current edge second tip
				x2=x
				y2=y

	#Check point outside
	print "Checking for points outside the domain..."
	i=0
	num=0
	while (i<numpy.size(geom.Vertices,axis=0)):

		#vertex counter
		i+=1

		#Get coordinates
		x=geom.Vertices[i,0]
		y=geom.Vertices[i,1]
		color=geom.Vertices[i,2]

		#Check that the point is inside the domain
		if color!=1 and not ContourToNodes(x,y,outline[0],1):

			#Remove points from list of Vertices
			num+=1
			geom.Vertices[i,:]=[]

			#update edges
			posedges=numpy.nonzero(geom.Edges==i)
			geom.Edges[posedges[0],:]=[]
			posedges=numpy.nonzero(geom.Edges>i)
			geom.Edges[posedges]=geom.Edges[posedges]-1

			#update counter
			i-=1

	if num:
		print "WARNING: %d points outside the domain outline have been removed" % num

	"""
	%Check point spacing
	if ~isnan(tol),
		disp('Checking point spacing...');
		i=0;
		while (i<size(geom.Vertices,1)),

			%vertex counter
			i=i+1;

			%Get coordinates
			x1=geom.Vertices(i,1);
			y1=geom.Vertices(i,2);

			j=i; %test edges located AFTER i only
			while (j<size(geom.Vertices,1)),

				%vertex counter
				j=j+1;

				%Get coordinates
				x2=geom.Vertices(j,1);
				y2=geom.Vertices(j,2);

				%Check whether the two vertices are too close
				if ((x2-x1)**2+(y2-y1)**2<tol**2)

					%Remove points from list of Vertices
					geom.Vertices(j,:)=[];

					%update edges
					posedges=find(ismember(geom.Edges,j));
					geom.Edges(posedges)=i;
					posedges=find(geom.Edges>j);
					geom.Edges(posedges)=geom.Edges(posedges)-1;

					%update counter
					j=j-1;

				end
			end
		end
	end
	%remove empty edges
	geom.Edges(find(geom.Edges(:,1)==geom.Edges(:,2)),:)=[];
	"""
	return geom
# }}}

