from scipy.io.netcdf import netcdf_file
import os
import numpy as np

class ncdfReader:
    def __init__(self,file):
        self.file = file
        self.data = netcdf_file(file,'r')
        vars = self.data.variables.keys()
        for key in vars:
            if self.data.variables[key].shape != ():
                setattr(self,key,1*self.data.variables[key][:])
            else:
                setattr(self,key,1*self.data.variables[key].getValue())

    def reset(self):
        '''Reset attributes to original values'''
        self.__init__(self.file)

    def write(self,file):
        '''Write to new file'''
        f = netcdf_file(file,'w',version=2)
        dim = ['tetinteriorsize', 'tetinterior', 'tetexteriorsize','tetexterior','ncoords', 'coord_size', 'nsurface_midpoint', 'surface_midpoint_size']
        v = ['tetrahedron_interior', 'tetrahedron_exterior', 'surface_midpoint','coords']
        #for i in self.data.dimensions.keys():
        for i in dim:
            f.createDimension(i, self.data.dimensions[i])
        #for i in self.data.variables.keys(): 
        for i in v:
            vars()[i] = f.createVariable(i, self.data.variables[i][:].dtype.str, self.data.variables[i].dimensions)
            vars()[i][:] = self.__dict__[i][:]
        f.close()

def convertParallel(mesh):
    os.system('ncdump '+str(mesh)+ ' > ' + mesh.replace('ncdf','cdl'))
    os.system('/Users/cho/work/Software/parallel-netcdf/bin/ncmpigen -o '+mesh + ' ' + mesh.replace('ncdf','cdl') )
    os.system('rm '+mesh.replace('ncdf','cdl'))
    os.system('/Users/cho/work/ace3p/bin/acdtool mesh fix '+mesh+' '+mesh)

def sidesetSurfaceCoords(meshName,tempMod,sideset):
    mesh = ncdfReader(meshName)
    temp = ncdfReader(tempMod)
    m1=mesh.coords[mesh.tetrahedron_exterior[mesh.tetrahedron_exterior[:,5]==sideset,1:6]][:,np.array([0,2,1]),:]
    m2=mesh.coords[mesh.tetrahedron_exterior[mesh.tetrahedron_exterior[:,6]==sideset,1:6]][:,np.array([0,3,2]),:]
    m3=mesh.coords[mesh.tetrahedron_exterior[mesh.tetrahedron_exterior[:,7]==sideset,1:6]][:,np.array([0,1,3]),:]
    m4=mesh.coords[mesh.tetrahedron_exterior[mesh.tetrahedron_exterior[:,8]==sideset,1:6]][:,np.array([1,2,3]),:]
    t1=temp.Temperature[mesh.tetrahedron_exterior[mesh.tetrahedron_exterior[:,5]==sideset,1:5]][:,np.array([0,2,1]),:][:,:,0]
    t2=temp.Temperature[mesh.tetrahedron_exterior[mesh.tetrahedron_exterior[:,6]==sideset,1:5]][:,np.array([0,3,2]),:][:,:,0]
    t3=temp.Temperature[mesh.tetrahedron_exterior[mesh.tetrahedron_exterior[:,7]==sideset,1:5]][:,np.array([0,1,3]),:][:,:,0]
    t4=temp.Temperature[mesh.tetrahedron_exterior[mesh.tetrahedron_exterior[:,8]==sideset,1:5]][:,np.array([1,2,3]),:][:,:,0]
    coords=np.concatenate([m1,m2,m3,m4])
    t=np.concatenate([t1,t2,t3,t4])
    return coords,t

def volume(mesh):
    vi=[]
    ve=[]
    for i in xrange(mesh.tetrahedron_interior.shape[0]):
        vi.append(1/6.0*np.abs(np.dot((mesh.coords[mesh.tetrahedron_interior[i][2]] - mesh.coords[mesh.tetrahedron_interior[i][1]]),np.cross((mesh.coords[mesh.tetrahedron_interior[i][3]] - mesh.coords[mesh.tetrahedron_interior[i][1]]),(mesh.coords[mesh.tetrahedron_interior[i][4]] - mesh.coords[mesh.tetrahedron_interior[i][1]])))))
    for i in xrange(mesh.tetrahedron_exterior.shape[0]):
        ve.append(1/6.0*np.abs(np.dot((mesh.coords[mesh.tetrahedron_exterior[i][2]] - mesh.coords[mesh.tetrahedron_exterior[i][1]]),np.cross((mesh.coords[mesh.tetrahedron_exterior[i][3]] - mesh.coords[mesh.tetrahedron_exterior[i][1]]),(mesh.coords[mesh.tetrahedron_exterior[i][4]] - mesh.coords[mesh.tetrahedron_exterior[i][1]])))))
    return vi,ve