// pde22_s2_eq.c   uniform grid boundary value PDE
//              known solution for(testing method
//              U(x,y) = sin^2(Pi*x)*sin^2(Pi*y)
//
//              differential equation to solve                 
//              d^2u/dx^2 + d^2u/dy^2 = Uxx(x,y) + Uyy(x,y) = c(x,y) given
//              Uxx(x,y) =  2*Pi^2*cos(2*Pi*x)*sin^2(Pi*y)
//              Uyy(x,y) = -2*Pi^2*sin^2(Pi*x)*cos(2*Pi*y)
//              U(x,y)=0 on boundary omega x=0, x=1, y=0, y=1
//
// 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 <time.h>
#include <math.h>
#include "nuderiv.h"
#include "simeq.h"

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

#define nx 13  // including boundary at 0 and nx-1, 0..ny-1 
#define ny 13  // including boundary at 0 and ny-1, 0..nx-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;

#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 sin(Pi*x)*sin(Pi*x)*sin(Pi*y)*sin(Pi*y);
} // end U 

static double Uxx(double x, double y)
{               // solution for(boundary and optional check 
  return 2.0*Pi*Pi*cos(2.0*Pi*x)*sin(Pi*y)*sin(Pi*y);
} // end U 

static double Uyy(double x, double y)
{               // solution for(boundary and optional check 
  return 2.0*Pi*Pi*sin(Pi*x)*sin(Pi*x)*cos(2.0*Pi*y);
} // end U 

static double c(double x, double y)
{               // RHS of differential equation to solve 
  return Uxx(x,y) + Uyy(x,y);
} //  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 
      nuderiv(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 
      nuderiv(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 
  {
    nuderiv(2,nx,i,xg,cxx);
    x = xg[i];
    for(j=1; j<ny-1; j++)
    {
      nuderiv(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;
  double start, now;

  printf("pde22_s2_eq.c running\n");
  printf("differential equation to solve\n");
  printf("d^2u/dx^2 + d^2u/dy^2 = c(x,y)\n");
  printf("Uxx(x,y) = 2*Pi^2*cos(2*Pi*x)*sin^2(Pi*y)\n");
  printf("Uyy(x,y) = -2*Pi^2*sin^2(Pi*x)*cos(2*Pi*y)\n");
  printf("c(x,y) = Uxx(x,y) + Uyy(x,y);\n");
  printf("uniform grid on rectangle 0,1 to 0,1\n");
  printf("known Solution, for testing method\n");
  printf("u(x,y) = sin^2(Pi*x)*sin^2(Pi*y)\n");
  printf("\n");
  start = (double)clock()/(double)CLOCKS_PER_SEC;
  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("c(0.5,0.5)=%f\n", c(0.5,0.5));
  printf("\n");

  init_matrix();
  printf("matrix initialized\n");
  // print_B(); 
  check_rows();
  simeqb(neqn, ut, us);
  check_soln();
  printf("\n");
  print();
  write_soln2(us, "pde22_s2_eq_c.dat");
  now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("finished pde22_s2_eq.c in %f seconds\n", now-start);
  return 0;
} // end pde22_s2_eq.c