from numpy import *
import TriMesh as tm
import NodeConnectivity as nc
import ElementConnectivity as ec

def triangle(md, domainname, resolution,riftname=''):
	#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). 

	#Check that mesh was not already run, and warn user: 
	if md.mesh.numberofelements != 0.:
		choice = input('This model already has a mesh. Are you sure you want to go ahead? (y/n)')
		if choice != 'y':
			print 'no meshing done ... exiting'
			return []
		
	area = resolution**2.

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

	#Fill in rest of fields:
	md.mesh.numberofelements = size(md.mesh.elements)
	md.mesh.numberofvertices = size(md.mesh.x)
	md.mesh.z = zeros(md.mesh.numberofvertices)
	md.mesh.vertexonboundary = zeros(md.mesh.numberofvertices)
	md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)] = 1.
	md.mesh.vertexonbed = ones(md.mesh.numberofvertices)
	md.mesh.vertexonsurface = ones(md.mesh.numberofvertices)
	md.mesh.elementonbed = ones(md.mesh.numberofelements)
	md.mesh.elementonsurface = ones(md.mesh.numberofelements)

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

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