import numpy
import sys
from model import *
from triangle import *
from setmask import *
from parameterize import *
from setflowequation import *
from EnumDefinitions import *
from solve import *
from MatlabFuncs import *

md=triangle(model(),'../Exp/Square.exp',150000.)
md=setmask(md,'../Exp/SquareShelf.exp','')
md=parameterize(md,'../Par/SquareSheetShelf.py')
md.extrude(5,1.)
md=setflowequation(md,'HO','all')
md.stressbalance.spcvx[numpy.nonzero(md.mesh.y>0.)]=float('NaN')
md.initialization.vx[:]=0.
md.initialization.vy[:]=0.
md.initialization.vel=numpy.zeros_like(md.initialization.vx)

md.cluster=generic('name',oshostname(),'np',3)
md=solve(md,StressbalanceSolutionEnum())
vel0=md.results.StressbalanceSolution.Vel

theta=30.*numpy.pi/180.
x=md.mesh.x
y=md.mesh.y
md.mesh.x=numpy.cos(theta)*x-numpy.sin(theta)*y
md.mesh.y=numpy.sin(theta)*x+numpy.cos(theta)*y

md.stressbalance.referential[:,0:3]=numpy.tile([numpy.cos(theta),numpy.sin(theta),0],(md.mesh.numberofvertices,1))
md.stressbalance.referential[:,3: ]=numpy.tile([0,0,1],(md.mesh.numberofvertices,1))
md=solve(md,StressbalanceSolutionEnum())
vel1=md.results.StressbalanceSolution.Vel

#plotmodel(md,'data',vel0,'data',vel1,'data',vel1-vel0,'title','Cartesian CS','title','Rotated CS','title','difference','view#all',2)
print "Error between Cartesian and rotated CS: %g" % (numpy.max(numpy.abs(vel0-vel1))/(numpy.max(numpy.abs(vel0))+sys.float_info.epsilon))

#Fields and tolerances to track changes
field_names     =['vel1']
field_tolerances=[1e-9]
field_values=[\
	vel1, \
	];