def __eval_fexpl(self, u, t):
        """
        Helper routine to evaluate the explicit part of the RHS

        Args:
            u: current values (not used here)
            t: current time

        Returns:
            explicit part of RHS
        """

        fexpl = mesh(self.nvars)

        # Copy values of u into pyClaw state object
        self.state.q[0, :] = u.values[0, :]

        # Evaluate right hand side
        tmp = self.solver.dqdt(self.state)
        fexpl.values[0, :] = tmp.reshape(self.nvars[1:])

        # Copy values of u into pyClaw state object
        self.state.q[0, :] = u.values[1, :]

        # Evaluate right hand side
        tmp = self.solver.dqdt(self.state)
        fexpl.values[1, :] = tmp.reshape(self.nvars[1:])

        # DEBUGGING
        # fexpl.values[0,:] = 0.0*self.mesh
        # fexpl.values[1,:] = 0.0*self.mesh
        return fexpl

    def solve_system(self,rhs,factor,u0,t):
        """
        Simple linear solver for (I-dtA)u = rhs

        Args:
            rhs: right-hand side for the nonlinear system
            factor: abbrev. for the node-to-node stepsize (or any other factor required)
            u0: initial guess for the iterative solver (not used here so far)
            t: current time (e.g. for time-dependent BCs)

        Returns:
            solution as mesh
        """

        b         = rhs.values.flatten()
        cb        = Callback()

        sol, info = LA.gmres( self.Id - factor*self.M, b, x0=u0.values.flatten(), tol=self.gmres_tol, restart=self.gmres_restart, maxiter=self.gmres_maxiter, callback=cb)
        # If this is a dummy call with factor==0.0, do not log because it should not be counted as a solver call
        if factor!=0.0:
          #print "SDC: Number of GMRES iterations: %3i --- Final residual: %6.3e" % ( cb.getcounter(), cb.getresidual() )
          self.logger.add(cb.getcounter())
        me        = mesh(self.nvars)
        me.values = unflatten(sol, 4, self.N[0], self.N[1])

        return me

    def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-dtA)u = rhs

        Args:
            rhs: right-hand side for the nonlinear system
            factor: abbrev. for the node-to-node stepsize (or any other factor required)
            u0: initial guess for the iterative solver (not used here so far)
            t: current time (e.g. for time-dependent BCs)

        Returns:
            solution as mesh
        """

        M1 = sp.hstack((sp.eye(self.nvars[1]), -factor * self.A))
        M2 = sp.hstack((-factor * self.A, sp.eye(self.nvars[1])))
        M = sp.vstack((M1, M2))

        b = np.concatenate((rhs.values[0, :], rhs.values[1, :]))

        sol = LA.spsolve(M, b)

        me = mesh(self.nvars)
        me.values[0, :], me.values[1, :] = np.split(sol, 2)

        return me

    def __eval_fexpl(self,u,t):
        """
        Helper routine to evaluate the explicit part of the RHS

        Args:
            u: current values (not used here)
            t: current time

        Returns:
            explicit part of RHS
        """


        fexpl        = mesh(self.nvars)

        # Copy values of u into pyClaw state object
        self.state.q[0,:,:] = u.values[0,:,:]

        # Evaluate right hand side
        self.solver.before_step(self.solver, self.state)
        tmp = self.solver.dqdt(self.state)
        
        fexpl.values[0,:,:] = unflatten(tmp, 1, self.nvars[1], self.nvars[2])


        # Copy values of u into pyClaw state object
        #self.state.q[0,:,:] = u.values[1,:,:]

        # Evaluate right hand side
        #tmp = self.solver.dqdt(self.state)
        #fexpl.values[1,:,:] = tmp.reshape(self.nvars[1:])

        return fexpl

    def u_exact(self,t):
        """
        Routine to compute the exact solution at time t

        Args:
            t: current time

        Returns:
            exact solution
        """
        
        dtheta = 0.01
        H      = 10.0
        a      = 5.0
        x_c    = -50.0
        
        me        = mesh(self.nvars)
        me.values[0,:,:] = 0.0*self.xx
        me.values[1,:,:] = 0.0*self.xx
        #me.values[2,:,:] = 0.0*self.xx
        #me.values[3,:,:] = np.exp(-0.5*(self.xx-0.0)**2.0/0.15**2.0)*np.exp(-0.5*(self.zz-0.5)**2/0.15**2)
        #me.values[2,:,:] = np.exp(-0.5*(self.xx-0.0)**2.0/0.05**2.0)*np.exp(-0.5*(self.zz-0.5)**2/0.2**2)
        me.values[2,:,:] = dtheta*np.sin( np.pi*self.zz/H )/( 1.0 + np.square(self.xx - x_c)/(a*a))
        me.values[3,:,:] = 0.0*self.xx
        return me

def check_datatypes_mesh(init):
    import pySDC.datatype_classes.mesh as m


    m1 = m.mesh(init)
    m2 = m.mesh(m1)


    m1.values[:] = 1.0
    m2.values[:] = 2.0

    m3 = m1 + m2
    m4 = m1 - m2
    m5 = 0.1*m1
    m6 = m1

    m7 = abs(m1)

    m8 = m.mesh(m1)

    assert isinstance(m3,type(m1))
    assert isinstance(m4,type(m1))
    assert isinstance(m5,type(m1))
    assert isinstance(m6,type(m1))
    assert isinstance(m7,float)

    assert m2 is not m1
    assert m3 is not m1
    assert m4 is not m1
    assert m5 is not m1
    assert m6 is m1

    assert np.shape(m3.values) == np.shape(m1.values)
    assert np.shape(m4.values) == np.shape(m1.values)
    assert np.shape(m5.values) == np.shape(m1.values)

    assert np.all(m1.values==1.0)
    assert np.all(m2.values==2.0)
    assert np.all(m3.values==3.0)
    assert np.all(m4.values==-1.0)
    assert np.all(m5.values==0.1)
    assert np.all(m8.values==1.0)
    assert m7 >= 0

    def u_exact(self,t):
        """
        Routine to compute the exact solution at time t

        Args:
            t: current time

        Returns:
            exact solution
        """

        me = mesh(self.nvars)
        me.values = np.cos(2.0*np.pi*(self.mesh-self.c*t))
        return me

    def u_exact(self,t):
        """
        Routine to compute the exact solution at time t

        Args:
            t: current time

        Returns:
            exact solution
        """
        
        me             = mesh(self.nvars)
        me.values[0,:] = 0.5*u_initial( self.mesh - (self.cadv + self.cs)*t , self.waveno) - 0.5*u_initial( self.mesh - (self.cadv - self.cs)*t , self.waveno)
        me.values[1,:] = 0.5*u_initial( self.mesh - (self.cadv + self.cs)*t , self.waveno) + 0.5*u_initial( self.mesh - (self.cadv - self.cs)*t , self.waveno)
        return me

    def __eval_fimpl(self,u,t):
        """
        Helper routine to evaluate the implicit part of the RHS

        Args:
            u: current values
            t: current time (not used here)

        Returns:
            implicit part of RHS
        """

        fimpl = mesh(self.nvars)
        fimpl.values = self.A.dot(u.values)
        return fimpl

    def u_exact(self, t):
        """
        Dummy routine for the exact solution, currently only passes the initial values

        Args:
            t: current time
        Returns:
            mesh type containing the initial values
        """

        # thou shall not call this at time > 0
        assert t is 0
        me = mesh(2)
        me.values = self.u0
        return me

    def prolong_space(self, G):
        """
        Dummy prolongation routine

        Args:
            G: the coarse level data (easier to access than via the coarse attribute)
        """

        if isinstance(G, mesh):
            F = mesh(G)
        elif isinstance(G, rhs_imex_mesh):
            F = rhs_imex_mesh(G)
        else:
            print("Transfer error")
            exit()
        return F

    def solve_system(self,rhs,factor,u0,t):
        """
        Simple linear solver for (I-dtA)u = rhs

        Args:
            rhs: right-hand side for the nonlinear system
            factor: abbrev. for the node-to-node stepsize (or any other factor required)
            u0: initial guess for the iterative solver (not used here so far)
            t: current time (e.g. for time-dependent BCs)

        Returns:
            solution as mesh
        """

        me = mesh(self.nvars)
        me.values = LA.spsolve(sp.eye(self.nvars)-factor*self.A,rhs.values)
        return me

    def u_exact(self,t):
        """
        Routine to compute the exact solution at time t

        Args:
            t: current time

        Returns:
            exact solution
        """
        
        xc,yc = self.state.grid.p_centers
        me        = mesh(self.nvars)
        me.values[0,:,:] = np.sin(2*np.pi*xc)*np.sin(2*np.pi*yc)
        me.values[1,:,:] = np.sin(2*np.pi*xc)*np.sin(2*np.pi*yc)

        return me

    def eval_f(self,u,t):
        """
        Routine to compute the RHS for both components simultaneously

        Args:
            t: current time (not used here)
            u: the current values
        Returns:
            RHS, 2 components
        """

        x1 = u.values[0]
        x2 = u.values[1]
        f = mesh(2)
        f.values[0] = -x2 + x1*(1 - x1**2 - x2**2)
        f.values[1] = x1 + 3*x2*(1 - x1**2 - x2**2)
        return f

    def prolong_space(self,G):
        """
        Prolongation implementation

        Args:
            G: the coarse level data (easier to access than via the coarse attribute)
        """

        if isinstance(G,mesh):
            u_fine = mesh(self.init_c,val=0)
            u_fine.values = self.Pspace.dot(G.values)
        elif isinstance(G,rhs_imex_mesh):
            u_fine = rhs_imex_mesh(self.init_c)
            u_fine.impl.values = self.Pspace.dot(G.impl.values)
            u_fine.expl.values = self.Pspace.dot(G.expl.values)

        return u_fine

    def restrict_space(self,F):
        """
        Restriction implementation

        Args:
            F: the fine level data (easier to access than via the fine attribute)
        """

        if isinstance(F,mesh):
            u_coarse = mesh(self.init_c,val=0)
            u_coarse.values = self.Rspace.dot(F.values)
        elif isinstance(F,rhs_imex_mesh):
            u_coarse = rhs_imex_mesh(self.init_c)
            u_coarse.impl.values = self.Rspace.dot(F.impl.values)
            u_coarse.expl.values = self.Rspace.dot(F.expl.values)

        return u_coarse

    def restrict_space(self, F):
        """
        Dummy restriction routine

        Args:
            F: the fine level data (easier to access than via the fine attribute)

        """

        if isinstance(F, mesh):
            G = mesh(F)
        elif isinstance(F, rhs_imex_mesh):
            G = rhs_imex_mesh(F)
        else:
            print("Transfer error")
            exit()
        return G

    def __eval_fimpl(self,u,t):
        """
        Helper routine to evaluate the implicit part of the RHS

        Args:
            u: current values
            t: current time (not used here)

        Returns:
            implicit part of RHS
        """

        temp         = u.values.flatten()
        temp         = self.M.dot(temp)
        fimpl        = mesh(self.nvars,val=0.0)
        fimpl.values = unflatten(temp, 4, self.N[0], self.N[1])
        
        return fimpl

    def __eval_fexpl(self,u,t):
        """
        Helper routine to evaluate the explicit part of the RHS

        Args:
            u: current values (not used here)
            t: current time

        Returns:
            explicit part of RHS
        """
        
        # Evaluate right hand side
        fexpl        = mesh(self.nvars,val=0.0)
        temp         = u.values.flatten()
        temp         = self.D_upwind.dot(temp)
        fexpl.values = unflatten( temp, 4, self.N[0], self.N[1])
              
        return fexpl