/* pde3b.c  3D non equal h boundary value PDE, iterative solution  */
/*  u(x,y,z) = x^3  + y^3  +z^3  + 2x^2y + 3xy^2 + 4z^2x + 5y + 6  */
/*     du/dx = 3x^2 + 0    + 0   + 4xy   + 3y^2  + 4z^2  + 0  + 0  */
/* d^2u/dx^2 = 6x   + 0    + 0   + 4y    + 0     + 0     + 0  + 0  */
/*     du/dy = 0    + 3y^2 + 0   + 2x^2  + 6xy   + 0     + 5  + 0  */
/* d^2u/dy^2 = 0    + 6y   + 0   + 0     + 6x    + 0     + 0  + 0  */
/*     du/dz = 0    + 0    +3z^2 + 0     + 0     + 8zx   + 0  + 0  */
/* d^2u/dz^2 = 0    + 0    +6z   + 0     + 0     + 8x    + 0  + 0  */
/*                                                        */
/* solve d^2u/dx^2 + d^2u/dy^2 + d^2u/dz^2 = c(x,y,z)     */
/* by second order method on cube -1,-1,-1 to 1.2,1,1     */

/* second order numerical difference                      */
/* d^2u/dx^2 + d^2u/dy^2 + d^2u/dz^2 =                    */
/* 1/hx^2(u(x-hx,y,z)+u(x+hx,y,z)-2u(x,y,z)) +            */
/* 1/hy^2(u(x,y-hy,z)+u(x,y+hy,z)-2u(x,y,z)) +            */
/* 1/hz^2(u(x,y,z-hz)+u(x,y,z+hz)-2u(x,y,z)) + O(h^2)     */
/*  = c(x,y,z) = 20x + 10y + 6z                           */
/* solve for new u(x,y,z) =                               */
/* ((u(x-hx,y,z)+u(x+hx,y,z))/hx^2 +                      */
/*  (u(x,y-hy,z)+u(x,y+hy,z))/hy^2 +                      */
/*  (u(x,y,z-hz)+u(x,y,z+hz))/hz^2 - c(x,y,z)) /          */
/*  (2/hx^2 + 2/hy^2 + 2/hz^2)                            */

#include <stdio.h>
#undef  abs
#define abs(a) ((a)<0.0?(-(a)):(a))
#undef  min
#define min(a,b) ((a)<(b)?(a):(b))

static int nx = 11;
static int ny = 11;
static int nz = 9;
static double us[40][40][40]; /* solution being iterated        */
static double ut[40][40][40]; /* temporary for partial results  */
           /* nx-2 by ny-2 by nz-2 internal, non boundary cells */
static const double xmin = -1.0;
static const double ymin = -1.0;
static const double zmin = -1.0;
static const double xmax =  1.2;
static const double ymax =  1.0;
static const double zmax =  1.0;
static double hx; /* x step = (xmax-xmin)/(nx-1) */
static double hy; /* y step = (ymax-ymin)/(ny-1) */
static double hz; /* z step = (zmax-zmin)/(nz-1) */
static double diff=1.0e100;  /* sum abs difference form last iteration */
static double error;         /* sum of absolute exact-computed         */
static double avgerror;      /* average error                          */
static double maxerror;      /* maximum cell error                     */

static double u(double x, double y, double z) /* for boundary and check */
{
  return x*x*x + y*y*y + z*z*z + 2.0*x*x*y + 3.0*x*y*y +
         4.0*z*z*x + 5.0*y +6.0;
}

static double c(double x, double y, double z) /* from solve */
{
  return 20.0*x + 10.0*y + 6.0*z;
}

static double u000(int i, int j, int k) /* new iterated value */
{
  return ( (us[i-1][j][k] + us[i+1][j][k])/(hx*hx) + 
           (us[i][j-1][k] + us[i][j+1][k])/(hy*hy) +
           (us[i][j][k-1] + us[i][j][k+1])/(hz*hz) -
           c(xmin+i*hx,ymin+j*hy,zmin+k*hz))/
           (2.0/(hx*hx) + 2.0/(hy*hy) + 2.0/(hz*hz));
}

static void init_boundary()
{
  int i, j, k;
  
  /* set boundary on six faces */
  for(i=0; i<nx; i++)
  {
    for(j=0; j<ny; j++)
    {
      us[i][j][0]    = u(xmin+i*hx, ymin+j*hy, zmin);
      us[i][j][nz-1] = u(xmin+i*hx, ymin+j*hy, zmax);
    }
  }
  for(j=0; j<ny; j++)
  {
    for(k=0; k<nz; k++)
    {
      us[0][j][k]    = u(xmin, ymin+j*hy, zmin+k*hz);
      us[nx-1][j][k] = u(xmax, ymin+j*hy, zmin+k*hz);
    }
  }
  for(i=0; i<nx; i++)
  {
    for(k=0; k<nz; k++)
    {
      us[i][0][k]    = u(xmin+i*hx, ymin, zmin+k*hz);
      us[i][ny-1][k] = u(xmin+i*hx, ymax, zmin+k*hz);
    }
  }
  /* zero internal cells */
  for(i=1; i<nx-1; i++)
    for(j=1; j<ny-1; j++)
      for(k=1; k<nz-1; k++)
        us[i][j][k] = 0.0;
} /* end init_boundary */

static void iterate()
{
  int i, j, k;
  
  diff = 0.0;
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<ny-1; j++)
    {
      for(k=1; k<nz-1; k++)
      {
        ut[i][j][k] = u000(i,j,k); 
        diff = diff + abs(ut[i][j][k]-us[i][j][k]);
      }
    }
  }
  /* copy temp to solution */
  for(i=1; i<nx-1; i++)
    for(j=1; j<ny-1; j++)
      for(k=1; k<nz-1; k++)
        us[i][j][k] = ut[i][j][k];
} /* end iterate */

static void check() /* typically not known, instrumentation */
{
  int i, j, k;
  double uexact, err;
  
  error = 0.0;
  maxerror = 0.0;
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<ny-1; j++)
    {
      for(k=1; k<nz-1; k++)
      {
        uexact = u(xmin+i*hx,ymin+j*hy,zmin+k*hz);
        err = abs(us[i][j][k]-uexact);
        error = error + err;
        if(err>maxerror) maxerror=err;
      }
    }
  }
} /* end check */

static void printus()
{
  int i, j, k;
  
  printf("\n computed solution\n");
  for(k=0; k<nz; k++)
  {
    printf("u[x,y,%d] plane \n", k);
    for(i=0; i<nx; i++)
    {
      for(j=0; j<min(ny,11); j++)
      {
        printf("%6.3f ", us[i][j][k]);
      }
      printf("\n");
    }
  }
  printf("\n");
} /* end printus */

static void printexact() /* typically not known */
{
  int i, j, k;
  
  printf("\n exact solution\n");
  for(k=0; k<nz; k++)
  {
    printf("u[x,y,%d] plane \n", k);
    for(i=0; i<nx; i++)
    {
      for(j=0; j<min(ny,11); j++)
      {
        printf("%6.3f ", u(xmin+i*hx,ymin+j*hy,zmin+k*hz));
      }
      printf("\n");
    }
  }
  if(min(ny,11)<ny) printf("Not all columns printed \n");
  printf("\n");
} /* end printexact */

static void printdiff()
{
  int i, j, k;
  
  printf("\n difference exact-computed solution\n");
  printf("\n");
  for(k=1; k<nz-1; k++)
  {
    printf("u[x,y,%d] plane \n", k);
    for(i=1; i<nx-1; i++)
    {
      printf("       ");
      for(j=1; j<min(ny-1,10); j++)
      {
        printf("%7.4f", u(xmin+i*hx,ymin+j*hy,zmin+k*hz)-us[i][j][k]);
      }
      printf("\n");
    }
    printf("\n");
  }
} /* end printdiff */

int main(int argc, char * argv[])
{
  int iter=0;
  
  printf("pde3b.c running\n");
  hx = (xmax-xmin)/(nx-1); /* step size in x */
  hy = (ymax-ymin)/(ny-1); /* step size in y */
  hz = (zmax-zmin)/(nz-1); /* step size in z */
  printf("xmin=%g, xmax=%g, nx=%d, hx=%g \n", xmin, xmax, nx, hx);
  printf("ymin=%g, ymax=%g, ny=%d, hy=%g \n", ymin, ymax, ny, hy);
  printf("zmin=%g, zmax=%g, nz=%d, hz=%g \n", zmin, zmax, nz, hz);
  
  init_boundary();
  printf(" initial conditions set\n");
  
  while(diff>0.00001 && iter<149)
  {
    iterate();
    check();
    avgerror = error/((nx-2)*(ny-2)*(nz-2));
    printf("iter=%d, diff=%g, error=%g, avg=%g, max=%g\n",
           iter, diff, error, avgerror, maxerror);
    iter++;
    if(iter%50==0) { printus(); printdiff(); }
  }
  printexact();
  printus();
  printdiff();
  return 0;
} /* end main of pde3b.c */