#include "bt.h"

void exact_solution(double xi, double eta, double zeta, double dtemp[5],
                    __global double *g_ce);

//---------------------------------------------------------------------
// first store the "interpolated" values everywhere on the grid    
//---------------------------------------------------------------------
__kernel void initialize2(__global double *g_u,
                          __global double *g_ce,
                          int gp0,
                          int gp1,
                          int gp2)
{
  int i, j, k, m, ix, iy, iz;
  double xi, eta, zeta, Pface[2][3][5], Pxi, Peta, Pzeta;

  k = get_global_id(2);
  j = get_global_id(1);
  i = get_global_id(0);
  if (k >= gp2 || j >= gp1 || i >= gp0) return;

  __global double (*u)[JMAXP+1][IMAXP+1][5] =
    (__global double (*)[JMAXP+1][IMAXP+1][5])g_u;

  zeta = (double)(k) * dnzm1;
  eta = (double)(j) * dnym1;
  xi = (double)(i) * dnxm1;

  for (ix = 0; ix < 2; ix++) {
    exact_solution((double)ix, eta, zeta, &Pface[ix][0][0], g_ce);
  }

  for (iy = 0; iy < 2; iy++) {
    exact_solution(xi, (double)iy , zeta, &Pface[iy][1][0], g_ce);
  }

  for (iz = 0; iz < 2; iz++) {
    exact_solution(xi, eta, (double)iz, &Pface[iz][2][0], g_ce);
  }

  for (m = 0; m < 5; m++) {
    Pxi   = xi   * Pface[1][0][m] + (1.0-xi)   * Pface[0][0][m];
    Peta  = eta  * Pface[1][1][m] + (1.0-eta)  * Pface[0][1][m];
    Pzeta = zeta * Pface[1][2][m] + (1.0-zeta) * Pface[0][2][m];

    u[k][j][i][m] = Pxi + Peta + Pzeta - 
                    Pxi*Peta - Pxi*Pzeta - Peta*Pzeta + 
                    Pxi*Peta*Pzeta;
  }
}

