/* pde_bihar2d_eq.c  discretization solution biharmonic 4th order
 *                   non uniform grid                                                                   
 * solve  uxxxx(x,y) + 2*uxxyy(x,y) + uyyyy(x,y) + 
 *        2*uxx(x,y) + 2*uyy(x,y)  + 2*u(x,y) = f(x,y)
 *
 *  f(x,y) := 2.4 + 0.08*x*x + 0.32*y*y + 1.6*x*y + 0.02*x*x*x*x +
 *            0.04*y*y*y*y - 0.08*x*x*y*y + 0.2*x*x*x*y + 0.4*y*y*y*x
 *
 *  in the rectangle  xmin< x <xmax and  ymin< y <ymax  ,
 *
 *     subject to boundary conditions of the first
 *     kind, Dirichlet
 *
 *  Dirichlet boundary conditions are computed using:
 *  ub(x,y) := 0.01*x*x*x*x + 0.02*y*y*y*y - 0.04*x*x*y*y +
 *             0.1*x*x*x*y + 0.2*y*y*y*x - x*y +1
 *
 *  uses   simeq.c  nuderiv.c 
 */

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include "simeq.h"
#include "nuderiv.h"

#undef  abs
#define abs(a) ((a)<0.0?(-(a)):(a))
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))

#define nx  8
#define ny  8

static double xmin = 0.0; /* problem parameters */
static double xmax = 4.0; 
static double ymin = 0.0;
static double ymax = 4.0;

/* global i, j, nx, ii, jj, ny, xg, yg needed by functions*/
static int i, j, ii, jj;
static double xg[nx]; /* X grig */
static double yg[ny]; /* Y grid */
static double hx, hy;
static double ug[nx*ny];        /* solution */
static double kg[nx*ny][nx*ny]; /* simultaneous equations */
static double fg[nx*ny];        /* kg * ug = fg */ 
static double Ua[nx*ny];        /* analytic solution for error check */
static double maxerr;


/* function prototypes */
double F(double x, double y); /* RHS */
double ub(double x, double y);
void discretize();
void dset(int k, int kk);
void write_soln();
void check_soln();

int main(int argc, char * argv[])
{
    int nxy = nx*ny;
    double x, y;
    double val, err, avgerr;
    int stat;
    double B[nx*ny*(nx*ny+1)]; /* for simeqb */

    printf("pde_bihar2d_eq.c running.\n");
    printf("solve  uxxxx(x,y) + 2*uxxyy(x,y) + uyyyy(x,y) + \n");
    printf("       2*uxx(x,y) + 2*uyy(x,y)  + 2*u(x,y) = f(x,y)\n");
    printf(" \n");
    printf("f(x,y) := 2.4 + 0.08*x*x + 0.32*y*y + 1.6*x*y + 0.02*x*x*x*x +\n");
    printf("          0.04*y*y*y*y - 0.08*x*x*y*y + 0.2*x*x*x*y + 0.4*y*y*y*x\n");
    printf(" \n");
    printf("in the rectangle   xmin< x <xmax  and  ymin< y <ymax,\n");
    printf(" \n");
    printf("   subject to boundary conditions of the first kind, Dirichlet\n");
    printf("   ub(x,y) := 0.01*x*x*x*x + 0.02*y*y*y*y - 0.04*x*x*y*y +\n");
    printf("              0.1*x*x*x*y + 0.2*y*y*y*x - x*y + 1\n");
    printf("xmin=%7.3f, xmax=%7.3f, ymin=%7.3f, ymax=%7.3f\n",
           xmin, xmax, ymin, ymax);
    printf("nx=%d, ny=%d\n",nx, ny);
    printf("x grid and analytic solution at ymin \n");
    hx = (xmax-xmin)/(double)(nx-1);
    for(i=0; i<nx; i++) /* does not need to be uniform */
    {
      xg[i] = xmin+i*hx;
      printf("i=%d, xg[i]=%f, ub(xg[i],ymin)=%f\n",
             i, xg[i], ub(xg[i],ymin));
    }  
    printf(" \n");
    printf("y grid and analytic solution at xmin \n");
    hy = (ymax-ymin)/(double)(ny-1); /* does not need to be uniform */
    for(ii=0; ii<ny; ii++)
    {
      yg[ii] = ymin+ii*hy;
      printf("ii=%d, yg[ii]=%f, ub(xmin,yg[ii])=%f \n",
             ii, yg[ii], ub(xmin, yg[ii]));
    }  
    printf(" \n");
    /* put solution in Ua vector, not usually known */ 
    for(i=0; i<nx; i++)
    {
      x = xg[i];
      for(ii=0; ii<ny; ii++)
      {
	y = yg[ii];
        Ua[i*ny+ii] = ub(x,y);
      }
    }
    printf(" \n");

    /* initialize k  */ 
    for(i=0; i<nx; i++)
    {
      for(j=0; j<nx; j++)
      {
        for(ii=0; ii<ny; ii++)
        {
          for(jj=0; jj<ny; jj++)
          {
            kg[(i*ny+ii)][(j*ny+jj)] = 0.0;
	  }
        }
      }
    }

    /* set boundary conditions, four sides */ 
    for(i=0; i<nx; i++)
    {
      x = xg[i];
      ii=0;
      y = yg[ii];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = ub(x,y);
      ii=ny-1;
      y = yg[ii];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = ub(x,y);
    }
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      i=0;
      x = xg[i];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = ub(x,y);
      i=nx-1;
      x = xg[i];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = ub(x,y);
    }

    for(i=1; i<nx-1; i++) /* set remaining RHS */
    {
      x = xg[i];
      for(ii=1; ii<ny-1; ii++)
      {
        y = yg[ii];
        fg[i*ny+ii] = F(x,y);
      }
    }

    /* compute entries in stiffness, equation matrix */ 
    printf("compute stiffness, equation matrix \n");
    for(i=1; i<nx-1; i++)
    {
      for(ii=1; ii<ny-1; ii++)
      {
        /* compute discrete derivatives around (i,ii) DOF */  
        discretize(); /* everything global */
                      /* sets kg[(i*ny+ii)][(j*ny+jj)...] */ 
      } /* end ii loop */
    } /* end i loop */

    printf("k computed discrete matrix \n ");
    printf("f computed forcing function, RHS \n ");

    for(i=0; i<nxy; i++) /* move kg to B */
      for(ii=0; ii<nxy; ii++)
        B[i*(nxy+1)+ii] = kg[i][ii];
    for(i=0; i<nxy; i++) /* move fg to B */
      B[i*(nxy+1)+nxy] = fg[i];

    simeqb(nxy, B, ug); /* B is kg fg */

    write_soln(); /* to file x y soln */

    check_soln(); /* independent of nowing analytic solution */

    avgerr = 0.0;
    maxerr = 0.0;
    printf("ug computed by discretization, Ua analytic, error \n");
    for(i=0; i<nx; i++)
    {
      for(ii=0; ii<ny; ii++)
      {
        err = abs(ug[i*ny+ii]-Ua[i*ny+ii]);
        if(err>maxerr) maxerr = err;
        avgerr = avgerr + err;
        printf("ug[%d,%d]=%f, Ua=%f, err=%g \n",
               i,ii,ug[i*ny+ii],Ua[i*ny+ii],ug[i*ny+ii]-Ua[i*ny+ii]);
      }
    }
    printf("xmax=%f, ymax=%f, nx=%d, ny=%d \n",
           xmax,ymax,nx,ny);
    printf("                  maxerr=%g, avgerr=%g \n",
           maxerr,avgerr/(double)(nx*ny));
    printf("\n");

  printf("calling gnuplot \n");
  fflush(stdout);

  stat = execl("/usr/bin/gnuplot","-persist","pde_bihar2d_eq_c.plot", NULL);
  /* only returns upon error, no return on sucess */
  printf("returned from gnuplot with stat= %d \n", stat);
  return 0;
} /* end pde_bihar2d_eq main */

/* PDE problem definition functions  f and ub */

double F(double x, double y)
{
  return 2.4 + 0.08*x*x + 0.32*y*y + 1.6*x*y + 0.02*x*x*x*x +
         0.04*y*y*y*y - 0.08*x*x*y*y + 0.2*x*x*x*y + 0.4*y*y*y*x;
}

double ub(double x, double y)
{
  return 0.01*x*x*x*x + 0.02*y*y*y*y - 0.04*x*x*y*y +
         0.1*x*x*x*y + 0.2*y*y*y*x - x*y + 1.0;
}

void discretize()  /* everything global i,ii (x,y) DOF */
                   /* sets kg[(i*ny+ii)][(j*ny+jj)...] */
{
  if(i<1 || i>nx-1 || ii<1 || ii>ny-1) /* primeter is boundary */
  {
    printf("discretize problem i=%d, ii=%d\n", i, ii);
    return;
  }
  if(i==1)
  {
    if(ii==1)         dset(1, 1); /* offset is negative of x,y point */
    else if(ii==ny-2) dset(1, 3);
    else              dset(1, 2);
  }
  else if(ii==1)
  {
    if(i==1)          dset(1, 1);
    else if(i==nx-2)  dset(3, 1);
    else              dset(2, 1);
  }
  else if(i==nx-2)
  {
    if(ii==1)         dset(3, 1);
    else if(ii==ny-2) dset(3, 3);
    else              dset(3, 2);
  }
  else if(ii==ny-2)
  {
    if(i==1)          dset(1, 3);
    else if(i==nx-2)  dset(3, 3);
    else              dset(2, 3);
  }
  else                dset(2, 2); /* center case */
} /* end discretize */

void dset(int k, int kk) /* x point, y point negative for offset */
{
  double cxxxx[5];
  double cyyyy[5];
  double cxx  [5];
  double cyy  [5];
  double cxxyy[5][5];
  double pxg  [5];
  double pyg  [5];
  int m, mx, my;

  if(i-k<0 || i-k+4>nx-1 || ii-kk<0 || ii-kk+4>ny-1)
  {
    printf("dset problem k=%d, kk=%d  for  i=%d, ii=%d\n", k, kk, i, ii);
    return;
  }
  printf("dset k=%d, kk=%d, row=%d\n", k, kk, i*ny+ii);

  for(m=0; m<5; m++)
  {
    pxg[m] = xg[i-k+m];  /* offset */
    pyg[m] = yg[ii-kk+m];
  }
  nuderiv(4, 5, k, pxg, cxxxx);
  nuderiv(2, 5, k, pxg, cxx);
  nuderiv(4, 5, kk, pyg, cyyyy);
  nuderiv(2, 5, kk, pyg, cyy);
  for(mx=0; mx<5; mx++)
  {
    for(my=0; my<5; my++) cxxyy[mx][my] = cxx[mx]*cyy[my];
  } 

  for(mx=0; mx<5; mx++) /* process all 25 cells */
  {
    for(my=0; my<5; my++)
    {
      /* optional, if boundary not in A and U
       * if(my==ii && (i-k+mx==0 || i-k+mx==nx-1)) // x boundary
       * {
       *   fg[i*ny+ii] -= ub(xg[i-k+mx],yg[ii-kk+my])*cxxxx[mx];
       *   fg[i*ny+ii] -= 2.0*ub(xg[i-k+mx],yg[ii-kk+my])*cxx[mx];
       * }
       * else if(mx==i && (ii-kk+my==0 || ii-kk+my==ny-1)) // y boundary
       * {
       *   fg[i*ny+ii] -= ub(xg[i-k+mx],yg[ii-kk+my])*cyyyy[my];
       *   fg[i*ny+ii] -= 2.0*ub(xg[i-k+mx],yg[ii-kk+my])*cyy[my];
       * }
       * else if(i-k+mx==0 || i-k+mx==nx-1 || ii-kk+my==0 || ii-kk+my==ny-1)
       * {                                                     // xxyy boundary
       *   fg[i*ny+ii] -= 2.0*ub(xg[i-k+mx],yg[ii-kk+my])*cxxyy[mx][my];
       * }
       *}
       *else
       *{
       */
      if(kk==my) /* Uxxxx and Uxx */
      {
        kg[(i*ny+ii)][(i-k+mx)*ny+ii] += cxxxx[mx];
        kg[(i*ny+ii)][(i-k+mx)*ny+ii] += 2.0*cxx[mx];
      }
      if(k==mx) /* Uyyyy and Uyy */
      {
        kg[(i*ny+ii)][i*ny+(ii-kk+my)] += cyyyy[my];
        kg[(i*ny+ii)][i*ny+(ii-kk+my)] += 2.0*cyy[my];
      }
      /* Uxxyy */
      kg[(i*ny+ii)][(i-k+mx)*ny+(ii-kk+my)] += 2.0*cxxyy[mx][my];
    } /* end my */
  } /* end mx */
  kg[(i*ny+ii)][(i*ny+ii)] += 2.0; /* U */
} /* end dset */

void write_soln()
{
  FILE * fout;

  fout = fopen("pde_bihar2d_eq_c.dat", "w");
  if(fout==NULL)
  {
    printf("file pde_bihar2d_eq_c.dat can not be opened for writing\n");
    return;
  }
  printf("writing pde_bihar2d_eq_c.dat \n");
  for(i=0; i<nx; i++)
  {
    for(ii=0; ii<ny; ii++)
    {
      fprintf(fout," %f %f %f \n",xg[i],yg[ii],ug[i*ny+ii]);
    } /* end y */
    fprintf(fout," \n");
  } /* end x */
  fflush(fout);
  fclose(fout);
  printf("finished writing pde_bihar2d_eq_c.dat \n");
} /* end write_soln */

void check_soln()
{
               /* for all interior points */
  double x, y, err, maxerr, avgerr, rmserr;
  double sumerr, sumsqerr, eval_deriv;
  double cxx[nx][nx];
  double cxxxx[nx][nx];
  double cyy[ny][ny];
  double cyyyy[ny][ny];
  double tcxxyy[ny];
  int cnt = 0;
  int k, m;

  for(i=1; i<nx-1; i++)
  {
    nuderiv(2,nx,i,xg,cxx[i]);
    nuderiv(4,nx,i,xg,cxxxx[i]);
  }
  for(ii=1; ii<ny-1; ii++)
  {
    nuderiv(2,ny,ii,yg,cyy[ii]);
    nuderiv(4,ny,ii,yg,cyyyy[ii]);
  }

  sumerr = 0.0;
  sumsqerr = 0.0;
  maxerr = 0.0;
  for(i=1; i<nx-1; i++) /* full equation */
  {
    x = xmin + (double)i*hx;
    for(j=1; j<ny-1; j++)
    {
      y = ymin + (double)j*hy;
      eval_deriv = 2.0*ug[i*ny+j]; /* PDE u term */
        
      for(k=0; k<nx; k++)
      {
	eval_deriv += 2.0*cxx[i][k]*ug[k*ny+j]; /* PDE  uxx term */
	eval_deriv += cxxxx[i][k]*ug[k*ny+j];   /* PDE  uxxxx term */
      } 
      for(k=0; k<ny; k++)
      {
	eval_deriv += 2.0*cyy[j][k]*ug[i*ny+k];  /* PDE uyy term */
	eval_deriv += cyyyy[j][k]*ug[i*ny+k];    /* PDE uyyyy term */
      }
      /*  + 2*uxxyy */
      for(jj=0; jj<ny; jj++) /* walk jj in y, computing uxx */
      {
	tcxxyy[jj] = 0.0; /* numeric  uxx part of uxxyy */
        for(k=0; k<nx; k++)
        {
	  tcxxyy[jj] += cxx[i][k]*ug[k*ny+jj];
	}
      }
      for(k=0; k<ny; k++)  /* uyy part of uxxyy */
      {
        eval_deriv += 2.0*cyy[j][k]*tcxxyy[k];  /* PDE uxxyy term */
      }

      err = eval_deriv-F(x,y);
      err = abs(err);
      sumerr += err;
      sumsqerr += err*err;
      maxerr = max(err,maxerr);
      cnt++;  
    } /* end y */
  } /* end x */
  rmserr = sqrt(sumsqerr/(double)cnt);
  avgerr = sumerr/(double)cnt;
  printf("check_soln values \n");
  printf("maxerr=%g, rmserr=%g, avgerr=%g \n",maxerr,rmserr,avgerr);
} /* end check_soln */ 

/* end pde_bihar2d_eq.c */