import numpy

def GetAreas(index,x,y,z=numpy.array([])):
	"""
	GETAREAS - compute areas or volumes of elements

	   compute areas of triangular elements or volumes 
	   of pentahedrons

	   Usage:
	      areas  =GetAreas(index,x,y);
	      volumes=GetAreas(index,x,y,z);

	   Examples:
	      areas  =GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
	      volumes=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y,md.z);
	"""

	#get number of elements and number of nodes
	nels=numpy.size(index,axis=0)
	nods=numpy.size(x)

	#some checks
	if numpy.size(y)!=nods or (z and numpy.size(z)!=nods):
		raise TypeError("GetAreas error message: x,y and z do not have the same length.")
	if numpy.max(index)>nods:
		raise TypeError("GetAreas error message: index should not have values above %d." % nods)
	if (not z and numpy.size(index,axis=1)!=3):
		raise TypeError("GetAreas error message: index should have 3 columns for 2d meshes.")
	if (z and numpy.size(index,axis=1)!=6):
		raise TypeError("GetAreas error message: index should have 6 columns for 3d meshes.")

	#initialization
	areas=numpy.zeros(nels)
	x1=x[index[:,0].astype(int)-1]
	x2=x[index[:,1].astype(int)-1]
	x3=x[index[:,2].astype(int)-1]
	y1=y[index[:,0].astype(int)-1]
	y2=y[index[:,1].astype(int)-1]
	y3=y[index[:,2].astype(int)-1]

	#compute the volume of each element
	if not z:
		#compute the surface of the triangle
		areas=(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)))
	else:
		#V=area(triangle)*1/3(z1+z2+z3)
		thickness=numpy.mean(z[index[:,3:6].astype(int)-1])-numpy.mean(z[index[:,0:3].astype(int)-1])
		areas=(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)))*thickness

	return areas

