/* fem_bihar2dps_la.c  Galerkin FEM biharmonic 4th order
 *                                                                   
 * solve  uxxxx(x,y) + uyyyy(x,y) + 2*uxxyy(x,y) = f(x,y)
 *
 *  f(x,y) := 4.0*sin(x)*sin(y)
 *
 *  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) := sin(x)*sin(y)
 *
 *  high order shape function used
 *  hight order Gauss-Legendre quadrature used
 *  small degree of freedom
 *  uses laphi.c  simeq.c  deriv.c  gaulegf.c 
 */

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <time.h>
#include "simeq.h"
#include "laphi.h"
#include "deriv.h"
#include "gaulegf.h"

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

#define nx  9
#define ny  9
#define npx 9
#define npy 9

static double xmin = 0.0; /* problem parameters */
static double xmax = 3.14159; 
static double ymin = 0.0;
static double ymax = 3.14159;
static double start, now, prev; /* timing */

/* i, j, nx, ii, jj, ny, xg, yg needed by galk or galf, available at call*/
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 maxerr;


/* function prototypes */
double F(double x, double y); /* RHS */
double ub(double x, double y);
double galk(double x, double y); /* Galerkin k function for this problem */
double galf(double x, double y); /* Galerkin f function for this problem */
void write_soln();
void check_soln();

int main(int argc, char * argv[])
{
    int nxy = nx*ny;
    double x, y;
    /* npx, npy  Gauss Legendre integration order */
    double xx[npx];
    double wx[npx];
    double yy[npy];
    double wy[npy];
    double val, err, avgerr;
    int px, py;
    double kg[nx*ny][nx*ny];
    double fg[nx*ny];
    double Ua[nx*ny];
    int stat;

    printf("fem_bihar2dps_la.c fourth order biharmonic PDE\n");
    start = (double)time(NULL);
    prev = start;
    printf("solve  uxxxx(x,y) + 2*uxxyy(x,y) + uyyyy(x,y) = f(x,y)\n");
    printf("f(x,y) := 4*sin(x)*sin(y) \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) := sin(x)*sin(y) \n");
    printf(" High order shape functions used. \n");
    printf(" High order Gauss-Legendre quadrature used. \n");
    printf("%d  Degrees of Freedom \n");
    printf("xmin=%7.4f, xmax=%7.4f, ymin=%7.4f, ymax=%7.4f\n",
           xmin, xmax, ymin, ymax);
    printf("nx=%d, ny=%d\n",nx, ny);
    printf("x grid and analytic solution at ymax \n");
    hx = (xmax-xmin)/(double)(nx-1);
    for(i=0; i<nx; i++)
    {
      xg[i] = xmin+i*hx;
      printf("i=%d, xg[i]=%f, ub(xg[i],ymax)=%f\n",
             i, xg[i], ub(xg[i],ymax));
    }  
    printf(" \n");
    printf("y grid and analytic solution at xmax \n");
    hy = (ymax-ymin)/(double)(ny-1);
    for(ii=0; ii<ny; ii++)
    {
      yg[ii] = ymin+ii*hy;
      printf("ii=%d, yg[ii]=%f, ub(xmax,yg[ii])=%f \n",
             ii, yg[ii], ub(xmax, yg[ii]));
    }  
    printf(" \n");
    /* put solution in Ua vector */ 
    for(i=0; i<nx; i++)
    {
      x = xmin+(double)i*hx;
      for(ii=0; ii<ny; ii++)
      {
	y = ymin+(double)ii*hy;
        Ua[i*ny+ii] = ub(x,y);
      }
    }
    printf(" \n");

    /* initialize k and f */ 
    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;
	  }
        }
      }
    }
    for(i=0; i<nx; i++)
    {
      for(ii=0; ii<ny; ii++)
      {
        fg[i*ny+ii] = 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);
    }

    printf("calling gauleg xmin=%f, xmax=%f, npx=%d \n",
           xmin, xmax, npx);
    gaulegf(xmin, xmax, xx, wx, npx); /* xx and wx set up for integrations */ 
    printf("xx[1]=%f, xx[2]=%f, wx[1]=%f, wx[2]=%f \n",
           xx[1],xx[2],wx[1],wx[2]);
    printf("calling gauleg ymin=%f, ymax=%f, npy=%d \n",
           ymin,ymax,npy);
    gaulegf(ymin, ymax, yy, wy, npy); /* yy and wy set up for integrations */
    printf("yy[1]=%f, yy[2]=%f, wy[1]=%f, wy[2]=%f \n",
           yy[1],yy[2],wy[1],wy[2]);
    i=1; ii=1;  j=1; jj=1;
    val = galk(xx[2],yy[2]);
    printf("galk(xx[2],yy[2])=%g\n",val);
    val = galf(xx[2],yy[2]);
    printf("galf(xx[2],yy[2])=%g\n",val);

    now = (double)time(NULL);
    fprintf(stderr, "time to initialize = %f seconds \n", now-prev);
    printf("time to initialize = %f seconds \n", now-prev);
    prev = now;

    /* compute entries in stiffness matrix */ 
    printf("compute stiffness matrix \n");
    for(i=1; i<nx-1; i++)
    {
      for(j=0; j<nx; j++)
      {
        for(ii=1; ii<ny-1; ii++)
        {
          for(jj=0; jj<ny; jj++)
          {
            /* compute integral for k(i,j) */  
            val = 0.0;
            for(px=1; px<=npx; px++)
            {
              for(py=1; py<=npy; py++)
              {
                val = val + wx[px]*wy[py]*galk(xx[px],yy[py]);
	      }
            }
            kg[(i*ny+ii)][(j*ny+jj)] = val;
          } /* end jj loop */ 
        } /* end ii loop */
      } /* end j loop */
    } /* end i loop */

    /* compute integral for f(i) */  
    for(i=1; i<nx-1; i++)
    {
      for(ii=1; ii<ny-1; ii++)
      {
        val = 0.0;
        for(px=1; px<=npx; px++)
        {
          for(py=1; py<=npy; py++)
          {
            val = val + wx[px]*wy[py]*galf(xx[px],yy[py]);
          }
        }
        fg[i*ny+ii] = val;
      }
    }
    printf("k computed stiffness matrix, see above \n ");
    printf("f computed forcing function, see above \n ");
    now = (double)time(NULL);
    fprintf(stderr, "time for stiffness and forcing = %f seconds \n", now-prev);
    printf("time for stiffness and forcing = %f seconds \n", now-prev);
    prev = now;
    simeq(nx*ny, (double *)kg, fg, ug);

    now = (double)time(NULL);
    fprintf(stderr, "time for solve = %f seconds \n", now-prev);
    printf("time for solve = %f seconds \n", now-prev);
    prev = now;

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

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

    now = (double)time(NULL);
    fprintf(stderr, "time for solve = %f seconds \n", now-prev);
    printf("time for write and check solution = %f seconds \n", now-prev);
    prev = now;

    avgerr = 0.0;
    maxerr = 0.0;
    printf("ug computed Galerkin, 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, npx=%d, npy=%d \n",
           xmax,ymax,nx,ny,npx,npy);
    printf("                  maxerr=%g, avgerr=%g \n",
           maxerr,avgerr/(double)(nx*ny));
    printf("\n");

  now = (double)time(NULL);
  fprintf(stderr, "total CPU time = %f seconds \n", now-start);
  printf("total CPU time = %f seconds \n", now-start);
  printf("calling gnuplot \n");
  fflush(stdout);

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

/* PDE problem definition functions  f and ub */

double F(double x, double y)
{
  return 4.0*sin(x)*sin(y);
}

double ub(double x, double y)
{
  return sin(x)*sin(y);
}

double galk(double x, double y)
{           /* Galerkin k stiffness function for this problem      */ 
            /* solve  uxxxx(x,y) + 2*uxxyy(x,y) + uyyyy(x,y)  = f(x,y) */
  return ( phipppp(x,j,nx-1,xg)*     phi(y,jj,ny-1,yg) +
        2.0* phipp(x,j,nx-1,xg)*   phipp(y,jj,ny-1,yg) +
               phi(x,j,nx-1,xg)* phipppp(y,jj,ny-1,yg) )*
               phi(x,i,nx-1,xg)*     phi(y,ii,ny-1,yg);
} /* end galk used by integration */ 

double galf(double x, double y) /* Galerkin f function for this problem */
{
  return F(x,y)*phi(x,i,nx-1,xg)*phi(y,ii,ny-1,yg);
} /* end galf used by integration */ 

void write_soln()
{
  FILE * fout;

  fout = fopen("fem_bihar2dps_la_c.dat", "w");
  if(fout==NULL)
  {
    printf("file fem_bihar2dps_la_c.dat can not be opened for writing\n");
    return;
  }
  printf("writing fem_bihar2dps_la_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 fem_bihar2dps_la_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++)
  {
    rderiv(2,nx,i,hx,cxx[i]);
    rderiv(4,nx,i,hx,cxxxx[i]);
  }
  for(ii=1; ii<ny-1; ii++)
  {
    rderiv(2,ny,ii,hy,cyy[ii]);
    rderiv(4,ny,ii,hy,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 = 0.0; /* no PDE u term */
        
      for(k=0; k<nx; k++)
      {
	eval_deriv += cxxxx[i][k]*ug[k*ny+j];   /* PDE  uxxxx term */
      } 
      for(k=0; k<ny; k++)
      {
	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 fem_bihar2dps_la.c */