import numpy
from fielddisplay import fielddisplay
from EnumDefinitions import *
from checkfield import checkfield
from WriteData import WriteData
import MatlabFuncs as m

class initialization(object):
	"""
	INITIALIZATION class definition

	   Usage:
	      initialization=initialization();
	"""

	def __init__(self): # {{{
		
		self.vx            = float('NaN')
		self.vy            = float('NaN')
		self.vz            = float('NaN')
		self.vel           = float('NaN')
		self.pressure      = float('NaN')
		self.temperature   = float('NaN')
		self.waterfraction = float('NaN')
		self.watercolumn   = float('NaN')
		self.sediment_head  = float('NaN')

		#set defaults
		self.setdefaultparameters()

		#}}}
	def __repr__(self): # {{{
		string='   initial field values:'

		string="%s\n%s"%(string,fielddisplay(self,'vx','x component of velocity [m/yr]'))
		string="%s\n%s"%(string,fielddisplay(self,'vy','y component of velocity [m/yr]'))
		string="%s\n%s"%(string,fielddisplay(self,'vz','z component of velocity [m/yr]'))
		string="%s\n%s"%(string,fielddisplay(self,'vel','velocity norm [m/yr]'))
		string="%s\n%s"%(string,fielddisplay(self,'pressure','pressure [Pa]'))
		string="%s\n%s"%(string,fielddisplay(self,'temperature','temperature [K]'))
		string="%s\n%s"%(string,fielddisplay(self,'waterfraction','fraction of water in the ice'))
		string="%s\n%s"%(string,fielddisplay(self,'watercolumn','thickness of subglacial water [m]'))
		string="%s\n%s"%(string,fielddisplay(self,'sediment_head','sediment water head of subglacial system [m]'))

		return string
		#}}}
	def setdefaultparameters(self): # {{{
		return self
	#}}}
	def checkconsistency(self,md,solution,analyses):    # {{{
		if StressbalanceAnalysisEnum() in analyses:
			if not numpy.any(numpy.logical_or(numpy.isnan(md.initialization.vx),numpy.isnan(md.initialization.vy))):
				md = checkfield(md,'fieldname','initialization.vx','NaN',1,'size',[md.mesh.numberofvertices])
				md = checkfield(md,'fieldname','initialization.vy','NaN',1,'size',[md.mesh.numberofvertices])
		if MasstransportAnalysisEnum() in analyses:
			md = checkfield(md,'fieldname','initialization.vx','NaN',1,'size',[md.mesh.numberofvertices])
			md = checkfield(md,'fieldname','initialization.vy','NaN',1,'size',[md.mesh.numberofvertices])
		if BalancethicknessAnalysisEnum() in analyses:
			md = checkfield(md,'fieldname','initialization.vx','NaN',1,'size',[md.mesh.numberofvertices])
			md = checkfield(md,'fieldname','initialization.vy','NaN',1,'size',[md.mesh.numberofvertices])
			#Triangle with zero velocity
			if numpy.any(numpy.logical_and(numpy.sum(numpy.abs(md.initialization.vx[md.mesh.elements-1]),axis=1)==0,\
			                               numpy.sum(numpy.abs(md.initialization.vy[md.mesh.elements-1]),axis=1)==0)):
				md.checkmessage("at least one triangle has all its vertices with a zero velocity")
		if ThermalAnalysisEnum() in analyses:
			md = checkfield(md,'fieldname','initialization.vx','NaN',1,'size',[md.mesh.numberofvertices])
			md = checkfield(md,'fieldname','initialization.vy','NaN',1,'size',[md.mesh.numberofvertices])
			if md.mesh.dimension()==3:
				md = checkfield(md,'fieldname','initialization.vz','NaN',1,'size',[md.mesh.numberofvertices])
			md = checkfield(md,'fieldname','initialization.pressure','NaN',1,'size',[md.mesh.numberofvertices])
		if (EnthalpyAnalysisEnum() in analyses and md.thermal.isenthalpy):
			md = checkfield(md,'fieldname','initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices])
			md = checkfield(md,'fieldname','initialization.watercolumn'  ,'>=',0,'size',[md.mesh.numberofvertices])
		if HydrologyShreveAnalysisEnum() in analyses:
			md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'size',[md.mesh.numberofvertices])
		return md
	# }}}
	def marshall(self,md,fid):    # {{{

		yts=365.0*24.0*3600.0

		WriteData(fid,'data',self.vx,'format','DoubleMat','mattype',1,'enum',VxEnum(),'scale',1./yts)
		WriteData(fid,'data',self.vy,'format','DoubleMat','mattype',1,'enum',VyEnum(),'scale',1./yts)
		WriteData(fid,'data',self.vz,'format','DoubleMat','mattype',1,'enum',VzEnum(),'scale',1./yts)
		WriteData(fid,'data',self.pressure,'format','DoubleMat','mattype',1,'enum',PressureEnum())
		WriteData(fid,'data',self.temperature,'format','DoubleMat','mattype',1,'enum',TemperatureEnum())
		WriteData(fid,'data',self.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum())
		WriteData(fid,'data',self.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum())
		WriteData(fid,'data',self.sediment_head,'format','DoubleMat','mattype',1,'enum',SedimentHeadEnum())

		if md.thermal.isenthalpy:
			tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure;
			pos  = numpy.nonzero(md.initialization.temperature > tpmp)[0]
			enthalpy      = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature);
			enthalpy[pos] = md.materials.heatcapacity*tpmp[pos].reshape(-1,1) - md.constants.referencetemperature + md.materials.latentheat*md.initialization.waterfraction[pos].reshape(-1,1)
			WriteData(fid,'data',enthalpy,'format','DoubleMat','mattype',1,'enum',EnthalpyEnum());

	# }}}
