/* pde_laplace.c   uniform grid boundary value PDE              */
/*              known solution for(testing method               */
/*              u(x,y) = cos(k*x)*cosh(k*y)                     */
/*                                                              */
/*              Laplace differential equation to solve          */
/*              d^2u/dx^2 + d^2u/dy^2 = c(x,y) = 0              */
/*                                                              */
/* uniform grid on rectangle 0 to 1, 0 to 1                     */
/* set up and solve system of linear equations                  */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mpf_nuderivd.h"
#define  digits 100
#include "simeq.h"

#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))

#define nx 17  /* including boundary at 0 and nx-1, 0..ny-1 */
#define ny 17  /* including boundary at 0 and ny-1, 0..ny-1 */
static double xmin = 0.0;
static double xmax = 1.0;
static double ymin = 0.0;
static double ymax = 1.0;
static double hx, hy;
static double k=6.28;

#define neqn (nx-2)*(ny-2) /* DOF number of equations */
#define cs neqn            /* RHS */

static double xg[nx];
static double yg[ny];
/* nx-2 by ny-2  internal, non boundary cells */
static double cxx[nx];
static double cyy[ny];
static double us[neqn]; /* solution vector */
static double ut[neqn*(neqn+1)]; /* includes RHS */

static double U(double x, double y)
{               /* solution for(boundary and optional check */
  return cos(k*x)*cosh(k*y);
} /* end U */

static double c(double x, double y)
{               /* RHS of differential equation to solve */
/*              d^2u/dx^2 + d^2u/dy^2 = c(x,y) = 0       */
  return 0.0;
} /*  end c */

static int S(int i, int j) /* zero based index includes boundary */
{                           /* this is row of ut[] */
  return (i-1)*(ny-2)+(j-1); /* index into us and ut, RHS, no boundary */
} /* end S */

static int SS(int i, int j, int ii, int jj)
                             /* zero based index includes boundary */
{
  return S(i,j)*(neqn+1)+S(ii,jj); /* ut[], RHS, no boundary */
} /* end SS */

static void print()
{
  double error = 0.0;
  double max_error = 0.0;
  double avg_error = 0.0;
  double val;
  int i, j;

  printf("\n");
  printf("exact solution u, computed us, error\n");
  for(i=0; i<nx; i++)
  {
    for(j=0; j<ny; j++)
    {
      if(i==0 || i==nx-1 || j==0 || j==ny-1)
      {                                    /* boundary */
        val = U(xg[i],yg[j]);
        error = 0.0;
      }
      else
      {
        val = us[S(i,j)];
        error = abs(val-U(xg[i],yg[j]));
        avg_error = avg_error + error;
        if(error > max_error)  max_error = error;
      }
      printf("xg[%2d]=%7.5f, yg[%2d]=%7.5f, u=%15.6e, us=%15.6e, err=%15.6e\n",
	     i, xg[i], j, yg[j], U(xg[i],yg[j]), val, error);
    }
  }
  printf("avg_error=%g, max_error=%g\n", avg_error/(double)neqn, max_error);
} /* end print */

static void init_matrix()
{
  int i, j, ii, jj, k, kj, ik;

  /* cs is RHS, constant column subscript */
  /* zero internal cells */
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<ny-1; j++)
    {
      for(ii=1; ii<nx-1; ii++)
      {
        for(jj=1; jj<ny-1; jj++)
	{
          ut[SS(i,j,ii,jj)] = 0.0;
        }
      }
      ut[S(i,j)*(neqn+1)+cs] = 0.0;
    }
  }

  printf("internal cells zeroed \n");
  for(i=1; i<nx-1; i++) /* inside boundary */
  {
    for(j=1; j<ny-1; j++)
    {
      ii = S(i,j); /* equation */
      /* for(each term */

      /* d^2 u(x,y)/dx^2 */
      mpf_nuderivd(digits, 2,nx,i,xg,cxx);
      for(k=0; k<nx; k++)
      {
        if(k==0 || k==nx-1)
	{
          ut[ii*(neqn+1)+cs] = ut[ii*(neqn+1)+cs] - cxx[k]*U(xg[k],yg[j]);
	}	 
        else
	{
          kj = S(k,j);
          ut[ii*(neqn+1)+kj] = ut[ii*(neqn+1)+kj] + cxx[k];
	}
      } /* k */

      /* d^2 u(x,y)/dy^2 */
      mpf_nuderivd(digits, 2,ny,j,yg,cyy);
      for(k=0; k<ny; k++)
      {
        if(k==0 || k==ny-1)
	{
          ut[ii*(neqn+1)+cs] = ut[ii*(neqn+1)+cs] - cyy[k]*U(xg[i],yg[k]);
        }
        else
	{
          ik = S(i,k);
          ut[ii*(neqn+1)+ik] = ut[ii*(neqn+1)+ik] + cyy[k];
	}
      } /* k */

      /* RHS constant */
      ut[ii*(neqn+1)+cs] = ut[ii*(neqn+1)+cs] + c(xg[i],yg[j]);
      /* end terms */
    } /* j */
  } /* i */
} /* end init_matrix */

static void print_B()
{
  int i, j, ii, jj;

  printf("matrix B includes RHS\n");
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<ny-1; j++)
    {
      for(ii=1; ii<nx-1; ii++)
      {
        for(jj=1; jj<ny-1; jj++)
	{
          printf("i=%2d, j=%2d, ii=%2d, jj=%2d, ut=%g\n",
		 i, j, ii, jj ,ut[SS(i,j,ii,jj)]);
	} /* jj */
      } /* ii */
      printf("RHS=                         %g\n", ut[S(i,j)*(neqn+1)+cs]);
    } /* j */
  } /* i */
  printf("\n");
  printf("\n");
} /*  end print_B */

static void check_rows() /* only for(checking */
{	/* row is in ut solution coordinates */
  int i, j, ii, jj, row, col;
  double sum;

  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<ny-1; j++)
    {
      row = S(i,j);
      sum = 0.0;
      for(ii=1; ii<nx-1; ii++)
      {
        for(jj=1; jj<ny-1; jj++)
	{
          col = S(ii,jj);
          sum = sum + ut[row*(neqn+1)+col]*U(xg[ii],yg[jj]); /* u is check solution */
	} /* ii */
      } /* jj */
      sum = sum - ut[row*(neqn+1)+cs];
      if(abs(sum)>0.001)
      {
        printf("row i=%2d, j=%2d, is bad err=%g\n", i, j, sum);
      }
    } /* j */
  } /*i */
} /* end check_rows */

static void check_soln()  /* check when solution unknown */
{
  double smaxerr, err, x, y, uxx, uyy;
  int i, j, ii, jj, k;

  /*  ux(x,y) + uy(x,y) = c(x,y) */
  printf("check_soln against PDE\n");
  smaxerr = 0.0;

  for(i=1; i<nx-1; i++)  /* through dof equations */
  {
    mpf_nuderivd(digits, 2,nx,i,xg,cxx);
    x = xg[i];
    for(j=1; j<ny-1; j++)
    {
      mpf_nuderivd(digits, 2,ny,j,yg,cyy);
      y = yg[j];
      err = -c(x,y); /* add other side of equation terms */
      uxx = 0.0;
      for(k=0; k<nx; k++) /* compute numerical derivative at (i,j) */
      {
        if(k==0 || k==nx-1)
	{
          uxx = uxx + cxx[k]*U(xg[k],yg[j]); /* boundary */
	}
        else
	{
          uxx = uxx + cxx[k]*us[S(k,j)];
        }
      }
      uyy = 0.0;
      for(k=0; k<ny; k++)
      {
        if(k==0 || k==ny-1)
	{
          uyy = uyy + cyy[k]*U(xg[i],yg[k]); /* boundary */
	}
        else
	{
          uyy = uyy + cyy[k]*us[S(i,k)];
	}
      }
      err = err + uxx + uyy;
      if(abs(err) > smaxerr) smaxerr = abs(err);
    } /* ii y */
  } /* i x */

  printf("check soln against PDE max error=%g\n", smaxerr);
} /* end Check_Soln */

/* write_soln2  for gnuplot */
static void  write_soln2(double Ux[], char file_name[])
{
  FILE * Fout;
  int i, j;

  Fout = fopen(file_name, "w");
  if(Fout==NULL)
  {
    printf("ERROR unable to open %s for writing \n", file_name);
    exit(1);
  }
  printf("writing %s\n", file_name);
  for(i=0; i<nx; i++)
  {
    for(j=0; j<ny; j++)
    {
      if(i==0 || i==nx-1 || j==0 || j==ny-1)
      {
        fprintf(Fout,"%f %f %f \n",
                xg[i], yg[j], U(xg[i],yg[j])); /* boundary */
      }
      else
      { 
        fprintf(Fout,"%f %f %f \n",
                xg[i],yg[j],Ux[S(i,j)]); /* computed solution */
      }
    } /* j */
    fprintf(Fout,"\n");
  } /* i */
  fclose(Fout);
  printf("finished writing %s\n", file_name);
}  /* end write_soln2 */

int main(int argc, char * argv[])
{
  int i, j;

  printf("pde_laplace.c running\n");
  printf("Laplace differential equation to solve\n");
  printf("d^2u/dx^2 + d^2u/dy^2 = c(x,y) = 0\n");
  printf("uniform grid on rectangle 0,1 to 0,1\n");
  printf("known Solution, for testing method\n");
  printf("u(x,y) = cos(kx)*cosh(kx)\n");
  printf("\n");
  hx = (xmax-xmin)/(double)(nx-1);
  hy = (ymax-ymin)/(double)(ny-1);

  for(i=0; i<nx; i++)
  {
    xg[i] = xmin+(double)i*hx;
    printf("xg(%2d)=%f\n", i, xg[i]);
  }
  for(j=0; j<ny; j++)
  {
    yg[j] = ymin+(double)j*hy;
    printf("yg(%2d)=%f\n", j, yg[j]);
  }
  printf("xmin=%f, xmax=%f, hx=%f, nx=%2d\n", xmin, xmax, hx, nx);
  printf("ymin=%f, ymax=%f, hy=%f, ny=%2d\n", ymin, ymax, hy, ny);
  printf("u(0.5,0.5)=%f\n", U(0.5,0.5));
  printf("k=%f\n", k);
  printf("\n");

  init_matrix();
  printf("matrix initialized\n");
  /* print_B(); */
  check_rows();
  simeqb(neqn, ut, us);
  check_soln();
  printf("\n");
  print();
  write_soln2(us, "pde_laplace_c.dat");
  printf("finished pde_laplace.c\n");
  return 0;
} /* end pde_laplace.c */