from SimNDT.engine.engineBase import EngineBase

#
# Setup and control of numerical simulation
# Input: Simulation.Im derived from Scenario.Iabs 
#
try:
    import SimNDT.engine.efit2dcython as EFIT2DCython
except:
    print("error in import for Serial EFIT2D (efit2dcython)")

from SimNDT.core.material import Material
from SimNDT.core.constants import BC

import copy
import numpy as np

global ErrorImportCL

try:
    import pyopencl as cl

    ErrorImportCL = False
except ImportError:
    ErrorImportCL = True

c11 = 4350 * 4350 * 7750
c12 = 4350 * 4350 * 7750 - 2260 * 2260 * 7750
c22 = c11
c44 = 2260 * 2260 * 7750
pzt = Material("pzt", 7750.0, c11, c12, c22, c44)


class EFIT2D(EngineBase):
    def __init__(self, simPack, Platform):
        EngineBase.__init__(self, simPack, Platform)
        self.max_value = 0.0
        self._Receiver = False

    def setup_CL(self):
        if self.Platform == "OpenCL":
            self.initFieldsCL()

    # Main setup of numerical model
    def materialSetup(self):
        print("EFIT2D: Creating material matrix from material library.")
        Materials = copy.deepcopy(self.simPack.Materials)
        MRI, NRI = np.shape(self.simPack.Simulation.Im)

        # Condicion de vacio si uno de los materiales es aire
        for n in range(len(Materials)):
            if Materials[n].Rho < 2.0:
                Materials[n].Rho = 10e23
                Materials[n].C11 = 1e-20
                Materials[n].C12 = 1e-20
                Materials[n].C22 = 1e-20
                Materials[n].C44 = 1e-20
                Materials[n].Eta_v = 1e-20
                Materials[n].Eta_s = 1e-20
            # verificar que no sea cero exactamente
            if Materials[n].C44 == 0:
                Materials[n].C44 = 1e-30
            if Materials[n].Eta_s == 0:
                Materials[n].Eta_s = 1e-30

        self.Rho = np.zeros((MRI, NRI), dtype=np.float32)
        self.C11 = np.zeros((MRI, NRI), dtype=np.float32)
        self.C12 = np.zeros((MRI, NRI), dtype=np.float32)
        self.C22 = np.zeros((MRI, NRI), dtype=np.float32)
        self.C44 = np.ones((MRI, NRI), dtype=np.float32) * 1e-30
        self.Eta_v = np.zeros((MRI, NRI), dtype=np.float32)
        self.Eta_s = np.zeros((MRI, NRI), dtype=np.float32)

        for n in range(len(Materials)):
            ind = np.nonzero(self.simPack.Simulation.Im == Materials[n].Label)
            self.Rho[ind] = Materials[n].Rho
            self.C11[ind] = Materials[n].C11
            self.C12[ind] = Materials[n].C12
            self.C22[ind] = Materials[n].C22
            self.C44[ind] = Materials[n].C44
            self.Eta_v[ind] = Materials[n].Eta_v
            self.Eta_s[ind] = Materials[n].Eta_s

        # Find the base materials (always label 0)
        for i, item in enumerate(Materials):
            if item.Label == 0:
                base_material = i

        # Set material in boundaries
        ind = np.nonzero(self.simPack.Simulation.Im == 255.0)
        self.Rho[ind] = Materials[base_material].Rho
        self.C11[ind] = Materials[base_material].C11
        self.C12[ind] = Materials[base_material].C12
        self.C22[ind] = Materials[base_material].C22
        self.C44[ind] = Materials[base_material].C44
        self.Eta_v[ind] = Materials[base_material].Eta_v
        self.Eta_s[ind] = Materials[base_material].Eta_s

    # Direct import of Rho, VL and VT material matrix
    def materialImport(self,Rho,VL,VT,Eta_v,Eta_s):
        print("EFIT2D: Importing material matrix.")
        eta_v=1e-30 
        eta_s = 1e-30
        MRI, NRI = np.shape(Rho)
        self.Rho = Rho.astype('float32')
        ind = np.nonzero(Rho < 2)
        self.Rho[ind] = 10e23
        self.Eta_v = np.ones((MRI, NRI), dtype=np.float32)*eta_v
        self.Eta_s = np.ones((MRI, NRI), dtype=np.float32)*eta_s
        Lam = Rho*(VL*VL-2*VT*VT)
        Mu  = Rho*(VT*VT)
        C11 = Lam+2*Mu
        C12 = Lam
        C22 = C11
        C44 = Mu
        self.C11 = C11.astype('float32')
        self.C12 = C12.astype('float32')
        self.C22 = C22.astype('float32')
        self.C44 = C44.astype('float32')
        self.C11[ind] = 1e-20
        self.C12[ind] = 1e-20
        self.C22[ind] = 1e-20
        self.C44[ind] = 1e-20
        print("Rho=["+str(self.Rho.min())+","+str(self.Rho.max())+"]")
        print("C11=["+str(self.C11.min())+","+str(self.C11.max())+"]")
        print("C12=["+str(self.C12.min())+","+str(self.C12.max())+"]")
        print("C22=["+str(self.C22.min())+","+str(self.C22.max())+"]")
        print("C44=["+str(self.C44.min())+","+str(self.C44.max())+"]")
        # Warning: If boundary is air the boundary must be set to base material, but applyBoundaries sets boundary again to 0; don't care about it here?
        
        

    def initFields(self):

        MRI, NRI = np.shape(self.simPack.Simulation.Im)
        Signal = self.simPack.Signal
        t = self.simPack.Simulation.t
        self.source = Signal.generate(t)

        self.Vx = np.zeros((MRI, NRI), dtype=np.float32)
        self.Vy = np.zeros((MRI, NRI), dtype=np.float32)
        self.Txx = np.zeros((MRI, NRI), dtype=np.float32)
        self.Txy = np.zeros((MRI, NRI), dtype=np.float32)
        self.Tyy = np.zeros((MRI, NRI), dtype=np.float32)
        self.DVx = np.zeros((MRI, NRI), dtype=np.float32)
        self.DVy = np.zeros((MRI, NRI), dtype=np.float32)
        self.ABS = np.ones((MRI, NRI), dtype=np.float32)
        self.SV = np.zeros((MRI, NRI), dtype=np.float32)

    def staggeredProp(self):

        MRI, NRI = np.shape(self.simPack.Simulation.Im)

        BXtemp = np.zeros((MRI, NRI), dtype=np.float32)
        BYtemp = np.zeros((MRI, NRI), dtype=np.float32)
        self.BX = np.zeros((MRI, NRI), dtype=np.float32)
        self.BY = np.zeros((MRI, NRI), dtype=np.float32)
        self.C44_Eff = np.zeros((MRI, NRI), dtype=np.float32)

        BXtemp[:, :] = 1.0 / self.Rho[:, :]
        BYtemp[:, :] = 1.0 / self.Rho[:, :]

        self.BX[:-2, :] = 0.5 * (BXtemp[1:-1, :] + BXtemp[:-2, :])
        self.BX[-2, :] = np.copy(BXtemp[-2, :])

        self.BY[:, :-2] = 0.5 * (BYtemp[:, 1:-1] + BYtemp[:, :-2])
        self.BY[:, -2] = np.copy(BYtemp[:, -2])

        self.C44_Eff[:-2, :-2] = 4. / (
            (1. / self.C44[:-2, :-2]) + (1. / self.C44[1:-1, :-2]) + (1. / self.C44[:-2, 1:-1]) + (
                1. / self.C44[1:-1, 1:-1]))

        self.Eta_vs = np.zeros((MRI, NRI), dtype=np.float32)
        self.Eta_ss = np.zeros((MRI, NRI), dtype=np.float32)
        self.Eta_vs = self.Eta_v + 2 * self.Eta_s
        self.Eta_ss[:-2, :-2] = 4. / (
            (1. / self.Eta_s[:-2, :-2]) + (1. / self.Eta_s[1:-1, :-2]) + (1. / self.Eta_s[:-2, 1:-1]) + (
                1. / self.Eta_s[1:-1, 1:-1]))

    def applyBoundaries(self):
        Materials = copy.deepcopy(self.simPack.Materials)
        MRI, NRI = np.shape(self.simPack.Simulation.Im)

        APARA = 0.015

        Top = int(self.simPack.Simulation.TapGrid[0])
        Bottom = int(self.simPack.Simulation.TapGrid[1])
        Left = int(self.simPack.Simulation.TapGrid[2])
        Right = int(self.simPack.Simulation.TapGrid[3])

        _Top = Top - np.arange(0, Top)
        _Bottom = np.arange(MRI - Bottom + 1, MRI) - MRI + Bottom - 1
        _Left = Left - np.arange(0, Left)
        _Right = np.arange(NRI - Right + 1, NRI) - NRI + Right - 1

        Boundaries = self.simPack.Boundary
        for boundary in Boundaries:
            if boundary.Name == "Top" and boundary.BC == BC.AirLayer:
                self.BX[0:Top, :] = 0.0
                self.BY[0:Top, :] = 0.0
                self.C11[0:Top, :] = 0.0
                self.C12[0:Top, :] = 0.0
                self.C22[0:Top, :] = 0.0
                self.C44_Eff[0:Top, :] = 0.0
                self.Eta_vs[0:Top, :] = 0.0
                self.Eta_ss[0:Top, :] = 0.0

            if boundary.Name == "Bottom" and boundary.BC == BC.AirLayer:
                self.BX[MRI - Bottom:, :] = 0.0
                self.BY[MRI - Bottom:, :] = 0.0
                self.C11[MRI - Bottom:, :] = 0.0
                self.C12[MRI - Bottom:, :] = 0.0
                self.C22[MRI - Bottom:, :] = 0.0
                self.C44_Eff[MRI - Bottom:, :] = 0.0
                self.Eta_vs[MRI - Bottom:, :] = 0.0
                self.Eta_ss[MRI - Bottom:, :] = 0.0

            if boundary.Name == "Left" and boundary.BC == BC.AirLayer:
                self.BX[:, 0:Left] = 0.0
                self.BY[:, 0:Left] = 0.0
                self.C11[:, 0:Left] = 0.0
                self.C12[:, 0:Left] = 0.0
                self.C22[:, 0:Left] = 0.0
                self.C44_Eff[:, 0:Left] = 0.0
                self.Eta_vs[:, 0:Left] = 0.0
                self.Eta_ss[:, 0:Left] = 0.0
            if boundary.Name == "Right" and boundary.BC == BC.AirLayer:
                self.BX[:, NRI - Right:] = 0.0
                self.BY[:, NRI - Right:] = 0.0
                self.C11[:, NRI - Right:] = 0.0
                self.C12[:, NRI - Right:] = 0.0
                self.C22[:, NRI - Right:] = 0.0
                self.C44_Eff[:, NRI - Right:] = 0.0
                self.Eta_vs[:, NRI - Right:] = 0.0
                self.Eta_ss[:, NRI - Right:] = 0.0

        for j in range(0, NRI):
            self.ABS[0:Top, j] = np.exp(- ((APARA * _Top) ** 2))
            self.ABS[MRI - Bottom + 1:, j] = np.exp(- ((APARA * _Bottom) ** 2))
        for i in range(0, MRI):
            self.ABS[i, 0:Left] = np.exp(- ((APARA * _Left) ** 2))
            self.ABS[i, NRI - Right + 1:] = np.exp(- ((APARA * _Right) ** 2))

    def sourceSetup(self):

        longitudinal = self.simPack.Source.Longitudinal
        shear = self.simPack.Source.Shear

        Pressure = self.simPack.Source.Pressure
        Displacement = self.simPack.Source.Displacement

        if Pressure and not Displacement:
            self.typeSource = 0
        elif not Pressure and Displacement:
            self.typeSource = 1
        elif Pressure and Displacement:
            self.typeSource = 2
        else:
            self.typeSource = 0

        if longitudinal and not shear:
            self.typeWave = np.int32(0)
        elif not longitudinal and shear:
            self.typeWave = np.int32(1)
        elif longitudinal and shear:
            self.typeWave = np.int32(2)
        else:
            self.typeWave = np.int32(0)

        self.Amplitude = self.simPack.Signal.Amplitude
        XL = self.simPack.Inspection.XL
        YL = self.simPack.Inspection.YL
        NX = np.size(XL, 0)

        IsPZT = self.simPack.Transducers[0].PZT

        size1 = int((self.simPack.Simulation.TapGrid[0]))
        size2 = int((self.simPack.Simulation.TapGrid[1]))

        for IT in range(-size1, 0):
            for m in range(0, NX):
                xl = int(XL[m, 0])
                yl = int(YL[m, 0])
                self.BX[xl + IT, yl] = 1.0 / pzt.Rho if IsPZT else 0
                self.BY[xl + IT, yl] = 1.0 / pzt.Rho if IsPZT else 0
                self.C11[xl + IT, yl] = pzt.C11 if IsPZT else 0
                self.C12[xl + IT, yl] = pzt.C12 if IsPZT else 0
                self.C22[xl + IT, yl] = pzt.C22 if IsPZT else 0
                self.C44[xl + IT, yl] = pzt.C44 if IsPZT else 0

        if self.simPack.Inspection.Name == "Transmission":
            for IT in range(1, size2):
                for m in range(0, NX):
                    xl = int(XL[m, 1])
                    yl = int(YL[m, 1])
                    self.BX[xl + IT, yl] = 1.0 / pzt.Rho if IsPZT else 0
                    self.BY[xl + IT, yl] = 1.0 / pzt.Rho if IsPZT else 0
                    self.C11[xl + IT, yl] = pzt.C11 if IsPZT else 0
                    self.C12[xl + IT, yl] = pzt.C12 if IsPZT else 0
                    self.C22[xl + IT, yl] = pzt.C22 if IsPZT else 0
                    self.C44[xl + IT, yl] = pzt.C44 if IsPZT else 0

    def simSetup(self):
        self.MRI, self.NRI = np.shape(self.simPack.Simulation.Im)
        XL = self.simPack.Inspection.XL
        YL = self.simPack.Inspection.YL

        self.NX = np.size(XL, 0)
        self.XL = np.copy(np.int32(XL[:, 0]))
        self.YL = np.copy(np.int32(YL[:, 0]))

        if self.simPack.Transducers[0].Window and not self.simPack.Transducers[0].PointSource:
            self.win = np.float32(1.0 - np.cos(2 * np.pi * np.arange(0, self.NX) / (self.NX + 1)))
        else:
            self.win = np.ones((self.NX,), dtype=np.float32)

    def initFieldsCL(self):

        self.Txx_buf = cl.Buffer(self.ctx, self.mf.READ_WRITE | self.mf.COPY_HOST_PTR, hostbuf=self.Txx)
        self.Tyy_buf = cl.Buffer(self.ctx, self.mf.READ_WRITE | self.mf.COPY_HOST_PTR, hostbuf=self.Tyy)
        self.Txy_buf = cl.Buffer(self.ctx, self.mf.READ_WRITE | self.mf.COPY_HOST_PTR, hostbuf=self.Txy)
        self.Vx_buf = cl.Buffer(self.ctx, self.mf.READ_WRITE | self.mf.COPY_HOST_PTR, hostbuf=self.Vx)
        self.Vy_buf = cl.Buffer(self.ctx, self.mf.READ_WRITE | self.mf.COPY_HOST_PTR, hostbuf=self.Vy)

        self.DVx_buf = cl.Buffer(self.ctx, self.mf.READ_WRITE | self.mf.COPY_HOST_PTR, hostbuf=self.DVx)
        self.DVy_buf = cl.Buffer(self.ctx, self.mf.READ_WRITE | self.mf.COPY_HOST_PTR, hostbuf=self.DVy)

        self.ABS_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.ABS)
        self.BX_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.BX)
        self.BY_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.BY)
        self.C11_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.C11)
        self.C12_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.C12)
        self.C44_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.C44_Eff)

        self.ETA_vs_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.Eta_vs)
        self.ETA_s_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.Eta_s)
        self.ETA_ss_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.Eta_ss)

        self.XL_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.XL)
        self.YL_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.YL)

        self.WIN_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.win)

        if not self.simPack.Inspection.Name == "Tomography":
            self._Receiver = True
            XL = self.simPack.Inspection.XL
            YL = self.simPack.Inspection.YL
            self.receiver_buf = cl.Buffer(self.ctx, self.mf.WRITE_ONLY | self.mf.COPY_HOST_PTR,
                                          hostbuf=self.receiver_signals)
            self.XXL = np.copy(np.int32(XL[:, 1]))
            self.YYL = np.copy(np.int32(YL[:, 1]))
            self.XXL_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.XXL)
            self.YYL_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.YYL)

        self.program = cl.Program(self.ctx, self.kernel()).build()

    def receiverSetup(self):
        TimeSteps = int(self.simPack.Simulation.TimeSteps)
        N_IR = np.size(self.simPack.Inspection.IR, 1)
        self.receiver_signals = np.zeros((TimeSteps, N_IR - 1), dtype=np.float32)

    # @blab+
    def clEnqueue(self, Var, Var_buf):
        # external access to Txx .. for snapshoting needs enqueueing of CL buffers
        cl.enqueue_copy(self.queue, Var, Var_buf)

    def receivers(self, Var, Var_buf):
        if self._Receiver:
            self.program.Receiver_EFIT2D(self.queue, (self.NX,), None, self.Txx_buf, self.receiver_buf,
                                         np.int32(self.n),
                                         self.XXL_buf, self.YYL_buf).wait()
        else:
            cl.enqueue_copy(self.queue, Var, Var_buf)
            self.receiver_signals[self.n, :] = self.simPack.Inspection.getReceivers(Var)

    def saveOutput(self):
        if self._Receiver:
            cl.enqueue_copy(self.queue, self.receiver_signals, self.receiver_buf).wait()

    def receiversSerial(self, Var):
        self.receiver_signals[self.n, :] = self.simPack.Inspection.getReceivers(Var)

    def run(self):

        y = np.float32(self.source[self.n])
        self.program.Velocity_EFIT2D_Voigt(self.queue, (self.NRI, self.MRI,), None,
                                           self.Txx_buf, self.Txy_buf, self.Tyy_buf,
                                           self.Vx_buf, self.Vy_buf, self.DVx_buf, self.DVy_buf,
                                           self.BX_buf, self.BY_buf, self.ABS_buf, self.ddx).wait()

        if self.typeSource == 1 or self.typeSource == 2:
            self.program.Source_Vel_EFIT2D(self.queue, (self.NX,), None,
                                           self.Vx_buf,
                                           self.Vy_buf,
                                           self.XL_buf,
                                           self.YL_buf,
                                           y, self.typeWave, self.WIN_buf).wait()

        self.program.Stress_EFIT2D_Voigt(self.queue, (self.NRI, self.MRI,), None,
                                         self.Txx_buf, self.Txy_buf, self.Tyy_buf,
                                         self.Vx_buf, self.Vy_buf, self.DVx_buf, self.DVy_buf,
                                         self.C11_buf, self.C12_buf, self.C44_buf, self.ETA_vs_buf,
                                         self.ETA_s_buf, self.ETA_ss_buf, self.ABS_buf).wait()

        if self.typeSource == 0 or self.typeSource == 2:
            self.program.Source_EFIT2D(self.queue, (self.NX,), None,
                                       self.Txx_buf,
                                       self.Tyy_buf,
                                       self.Txy_buf,
                                       self.XL_buf,
                                       self.YL_buf,
                                       y, self.typeWave, self.WIN_buf).wait()

        self.receivers(self.Txx, self.Txx_buf)

    def runGL(self):
        cl.enqueue_copy(self.queue, self.Vx, self.Vx_buf)
        cl.enqueue_copy(self.queue, self.Vy, self.Vy_buf)
        self.SV = np.sqrt(self.Vx ** 2 + self.Vy ** 2)
        value = np.max(np.abs(self.SV))
        self.max_value = self.max_value if self.max_value >= value else value
        self.SV = 20. * np.log10(np.abs(self.SV) / self.max_value + 1e-40)

    def runGLSerial(self, step=50):

        vx = np.copy(self.Vx)
        vy = np.copy(self.Vy)

        self.SV = np.sqrt(vx ** 2 + vy ** 2)
        value = np.max(np.abs(self.SV))
        self.max_value = self.max_value if self.max_value >= value else value
        self.SV = 20. * np.log10(np.abs(self.SV) / self.max_value + 1e-40)

    def runSerial(self):

        y = np.float32(self.source[self.n])

        self.Vx, self.Vy, self.DVx, self.DVy = EFIT2DCython.velocityVoigt(self.Txx, self.Txy, self.Tyy,
                                                                          self.Vx, self.Vy, self.DVx, self.DVy,
                                                                          self.BX, self.BY, self.ABS, self.ddx,
                                                                          self.dt)

        if self.typeSource == 1 or self.typeSource == 2:
            self.Vx, self.Vy = EFIT2DCython.sourceVel(self.Vx, self.Vy, self.XL, self.YL, y, self.typeWave, self.win,
                                                      self.dtdxx)

        self.Txx, self.Tyy, self.Txy = EFIT2DCython.stressVoigt(self.Txx, self.Txy, self.Tyy,
                                                                self.Vx, self.Vy, self.DVx, self.DVy,
                                                                self.C11, self.C12, self.C44_Eff,
                                                                self.Eta_vs, self.Eta_s, self.Eta_ss, self.ABS,
                                                                self.dtx)

        if self.typeSource == 0 or self.typeSource == 2:
            self.Txx, self.Tyy = EFIT2DCython.sourceStress(self.Txx, self.Tyy, self.XL, self.YL, y, self.typeWave,
                                                           self.win, self.dtdxx)

        self.receiversSerial(self.Txx)

    def kernel(self):

        macro = """
					 #define MRI		%s
					 #define NRI		%s
					 #define ind(i, j)	( ( (i)*NRI) + (j) )
					 #define dtx		%gf
					 #define dtdxx		%gf
					 #define Stencil	2
					 #define NX			%s
					 #define dt			%gf



		""" % (str(self.MRI), str(self.NRI), self.dtx, self.dtdxx,
               str(self.NX), self.dt)

        return macro + """


			__kernel void Receiver_EFIT2D(__global float *Buffer, __global float *receiver, const int t,
							              __global const int *XXL, __global const int *YYL){


			  __private float _tmp = 0.0f;


			   for (int i=0; i<get_global_size(0); ++i)
		       {
				 _tmp +=  Buffer[ind(XXL[i],YYL[i])];
			   }
			   receiver[t] = _tmp/(float)get_global_size(0);
            }



			__kernel void Source_Vel_EFIT2D(__global float *Vx, __global float *Vy, 
										__global const int	 *XL, __global const int *YL,  
										  const float source, const int Stype, __global const float *win){

					 uint m =  get_global_id(0);

					 int ix = XL[m];
					 int iy = YL[m];
					 if (Stype ==0){
					  Vx[ind(ix,iy)]   -= (source*dtdxx)*win[m];
					 }
					else if (Stype ==1){
					  Vy[ind(ix,iy)]   -= (source*dtdxx)*win[m];
					 }
					else if (Stype ==2){
					  Vx[ind(ix,iy)]   -= (source*dtdxx)*win[m];
					  Vy[ind(ix,iy)]   -= (source*dtdxx)*win[m];
					 }

			   }



			__kernel void Source_EFIT2D(__global float *Txx, __global float *Tyy, __global float *Txy,
										__global const int	 *XL, __global const int *YL,  
										  const float source, const int Stype,__global const float *win){

					 uint m =  get_global_id(0);

					 int ix = XL[m];
					 int iy = YL[m];
					 if (Stype ==0){
					  Txx[ind(ix,iy)]	-= (source*dtdxx)*win[m];
					  //Txy[ind(ix,iy)]	   = 0;
					 }
					else if (Stype ==1){
					  Tyy[ind(ix,iy)]	-= (source*dtdxx)*win[m];
					  //Txy[ind(ix,iy)]	   = 0;
					 }
					else if (Stype ==2){
					  Txx[ind(ix,iy)]	-= (source*dtdxx)*win[m];
					  Tyy[ind(ix,iy)]	-= (source*dtdxx)*win[m];
					  //Txy[ind(ix,iy)]	   = 0;
					 }

			   }


		  __kernel void Velocity_EFIT2D_Voigt( __global float *Txx, __global float *Txy, __global float *Tyy,
											   __global float *vx,	__global float *vy,
											   __global float *dvx, __global float *dvy,
											   __global const float *BX, __global const float *BY,
											   __global const float *ABS, const float ddx){

					int j = get_global_id(0);
					int i = get_global_id(1);


						  if (i <  MRI-1 && j>0){
								dvx[ind(i,j)]	 = (BX[ind(i,j)]*ddx)*( Txx[ind(i+1,j)] - Txx[ind(i,j)]	  + Txy[ind(i,j)] - Txy[ind(i,j-1)] );
								vx[ind(i,j)]	+= dt*dvx[ind(i,j)];
						  }

						  if (i > 0 &&	j< NRI-1){
								dvy[ind(i,j)]	 = (BY[ind(i,j)]*ddx)*( Txy[ind(i,j)]  - Txy[ind(i-1,j)] + Tyy[ind(i,j+1)] - Tyy[ind(i,j)]		 );
								vy[ind(i,j)]	+= dt*dvy[ind(i,j)];

						  }

						   barrier(CLK_GLOBAL_MEM_FENCE);
							// Apply absorbing boundary conditions
						   vx[ind(i,j)]	  *= ABS[ind(i,j)];
						   vy[ind(i,j)]	  *= ABS[ind(i,j)];


						 }


				__kernel void Stress_EFIT2D_Voigt(__global float *Txx, __global float *Txy, __global float *Tyy,
												  __global float *vx,  __global float *vy,
												  __global float *dvx, __global float *dvy,
												  __global const float *C11,   __global const float *C12, __global const float *C44,
												  __global const float *ETA_VS, __global const float *ETA_S, __global const float *ETA_SS,
												  __global const float *ABS) {


						int j = get_global_id(0);
						int i = get_global_id(1);


						if (i>0 &&	j>0 ){
								Txx[ind(i,j)] += (	( C11[ind(i,j)]* dtx )*(vx[ind(i,j)] - vx[ind(i-1,j)]) +
													( C12[ind(i,j)]* dtx )*(vy[ind(i,j)] - vy[ind(i,j-1)]) +
													( (ETA_VS[ind(i,j)])*dtx) * (dvx[ind(i,j)] - dvx[ind(i-1,j)]) +
													( (ETA_S [ind(i,j)])*dtx) * (dvy[ind(i,j)] - dvy[ind(i,j-1)]) );

								Tyy[ind(i,j)] += (	( C12[ind(i,j)]* dtx )*(vx[ind(i,j)] - vx[ind(i-1,j)]) +
													( C11[ind(i,j)]* dtx )*(vy[ind(i,j)] - vy[ind(i,j-1)]) +
													( (ETA_VS[ind(i,j)])*dtx) *(dvy[ind(i,j)] - dvy[ind(i,j-1)]) +
													( (ETA_S [ind(i,j)])*dtx) * (dvx[ind(i,j)] - dvx[ind(i-1,j)]) );
						 }

						 if (i<MRI-1  && j<NRI-1){
								Txy[ind(i,j)] +=	(  ( C44[ind(i,j)]	   * dtx )*(vx[ind(i,j+1)] - vx[ind(i,j)] + vy[ind(i+1,j)] - vy[ind(i,j)])		+
													   ( ETA_SS[ind(i,j)]  * dtx )*(dvx[ind(i,j+1)] - dvx[ind(i,j)] + dvy[ind(i+1,j)] - dvy[ind(i,j)] ) );
								 }

						  barrier(CLK_GLOBAL_MEM_FENCE);

						  Txx[ind(i,j)] *= ABS[ind(i,j)];
						  Tyy[ind(i,j)] *= ABS[ind(i,j)];
						  Txy[ind(i,j)] *= ABS[ind(i,j)];


				}



		"""
