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,'all','')
md=parameterize(md,'../Par/SquareShelf.py')
md=setflowequation(md,'SSA','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',2)
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')
print "Error between Cartesian and rotated CS: %g" % (numpy.max(numpy.abs(vel0-vel1))/(numpy.max(numpy.abs(vel0))+sys.float_info.epsilon))

#Now, put CS back to normal except on the side where the spc are applied
pos=numpy.nonzero(numpy.logical_or(x==0.,x==1000000.))[0]
md.stressbalance.referential[:]=float('NaN')
md.stressbalance.referential[pos,0:3]=numpy.tile([numpy.cos(theta),numpy.sin(theta),0],(len(pos),1))
md.stressbalance.referential[pos,3: ]=numpy.tile([0,0,1],(len(pos),1))
md=solve(md,StressbalanceSolutionEnum())
vel2=md.results.StressbalanceSolution.Vel

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

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