import numpy
from TriMesh import *
from NodeConnectivity import *
from ElementConnectivity import *
from MatlabFuncs import *

def triangle(md,domainname,*args):
	"""
	TRIANGLE - create model mesh using the triangle package

	   This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
	   where md is a @model object, domainname is the name of an Argus domain outline file, 
	   and resolution is a characteristic length for the mesh (same unit as the domain outline
	   unit). Riftname is an optional argument (Argus domain outline) describing rifts.

	   Usage:
	      md=triangle(md,domainname,resolution)
	   or md=triangle(md,domainname, resolution, riftname)

	   Examples:
	      md=triangle(md,'DomainOutline.exp',1000);
	      md=triangle(md,'DomainOutline.exp',1000,'Rifts.exp');
	"""

	#Figure out a characteristic area. Resolution is a node oriented concept (ex a 1000m  resolution node would 
	#be made of 1000*1000 area squares). 

	if len(args)==1:
		resolution=args[0]
		riftname=''
	if len(args)==2:
		riftname=args[0]
		resolution=args[1]

	#Check that mesh was not already run, and warn user: 
	if md.mesh.numberofelements:
		choice = raw_input('This model already has a mesh. Are you sure you want to go ahead? (y/n)')
		if not strcmp(choice,'y'):
			print 'no meshing done ... exiting'
			return None

	area = resolution**2

	#Mesh using TriMesh
	[md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMesh(domainname,riftname,area)

	#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(md.mesh.numberofvertices)
	md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1] = 1.
	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)

	#Now, build the connectivity tables for this mesh.
	[md.mesh.vertexconnectivity]= NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
	[md.mesh.elementconnectivity] = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)

	#type of model
	md.mesh.dimension = 2.

	return md

