#include "lu.h"

void exact(int i, int j, int k, double u000ijk[], __global double *g_ce);

__kernel void setiv(__global double *g_u,
                    __global double *g_ce,
                    int nx,
                    int ny,
                    int nz)
{
  int i, j, k, m;
  double xi, eta, zeta;
  double pxi, peta, pzeta;
  double ue_1jk[5], ue_nx0jk[5], ue_i1k[5];
  double ue_iny0k[5], ue_ij1[5], ue_ijnz[5];

  k = get_global_id(2) + 1;
  j = get_global_id(1) + 1;
  i = get_global_id(0) + 1;
  if (k >= (nz-1) || j >= (ny-1) || i >= (nx-1)) return;

  __global double (*u)[ISIZ2/2*2+1][ISIZ1/2*2+1][5] = 
    (__global double (*)[ISIZ2/2*2+1][ISIZ1/2*2+1][5])g_u;

  zeta = ( (double)k ) / (nz-1);
  eta = ( (double)j ) / (ny0-1);
  xi = ( (double)i ) / (nx0-1);
  exact(0, j, k, ue_1jk, g_ce);
  exact(nx0-1, j, k, ue_nx0jk, g_ce);
  exact(i, 0, k, ue_i1k, g_ce);
  exact(i, ny0-1, k, ue_iny0k, g_ce);
  exact(i, j, 0, ue_ij1, g_ce);
  exact(i, j, nz-1, ue_ijnz, g_ce);

  for (m = 0; m < 5; m++) {
    pxi =   ( 1.0 - xi ) * ue_1jk[m]
                  + xi   * ue_nx0jk[m];
    peta =  ( 1.0 - eta ) * ue_i1k[m]
                  + eta   * ue_iny0k[m];
    pzeta = ( 1.0 - zeta ) * ue_ij1[m]
                  + zeta   * ue_ijnz[m];

    u[k][j][i][m] = pxi + peta + pzeta
      - pxi * peta - peta * pzeta - pzeta * pxi
      + pxi * peta * pzeta;
  }
}

