/* fem_checke_la.c  Galerkin FEM
 *  solve  u'(x)=e^x  u(0)=1, u(1)=e  0 < x < 1 
 * initial try 10 points, then 20
 * investigate gauss-legendre, adaptive, and trapezoidal integration
 * solve for Uj
 * sum j=0,9 int x=0,1 phip(x,j)*phi(x,i) dx = int x=0,1 f(x)*phi(x,i)dx
 * k(i,j) = int x=0,1 phip(x,j)*phi(x,i) dx  (not i==0 or i==9, boundary)
 * f(i)   = int x=0,1 exp(x)*phi(x,i) dx
 * k * U = f, solve simultaneous equations for U
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "laphi.h"
#include "simeq.h"
#include "gaulegf.h"
#include "aquad3.h"

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

static double k[400]; /* i*n+j  Gauss-Legendre integration */
static double ks[400]; /* simple trapazoidal integration */
static double kd[400]; /* adaptive integration */
static double xg[20];
static double f[20];
static double u[20];
static double fs[20];
static double us[20];
static double fd[20];
static double ud[20];
static double Ua[20];


double F(double x) /* forcing function */
{
  return exp(x);
}

double uana(double x) /* analytic solution and boundary */
{
  return exp(x);
}

/* needed by galk or galk, available at call */
int i, j, n;
double galk(double x) /* Galerkin k stiffness function */
{
  return phip(x,j,n-1,xg)*phi(x,i,n-1,xg);
} /* end galk used by aquad3 */

double galf(double x) /* Galerkin f function */
{
  return F(x)*phi(x,i,n-1,xg);
} /* end galf used by aquad */

int main(int argc, char *argv[])
{
  /* i, j and n above, used by aquad3 */
  double xmin = 0.0;
  double xmax = 1.0;
  double x, h;
  double a, b, sum;
  double xx[100], ww[100]; /* for Gauss-Legendre */
  double val, err, avgerr, maxerr;
  double eps = 1.0e-7;
  int p, np;

  printf("fem_checke_la.c running \n");
  printf("Given du/dx = exp(x)  0<x<1  Boundary u(0) = 1, u(1) = e \n");
  printf("Analytic solution u(x) = exp(x) \n");

  printf("Try 10 points \n");
  n = 10;
  h = (xmax-xmin)/(double)(n-1);
  printf("x grid and analytic solution \n");
  for(i=0; i<n; i++)
  {
    xg[i] = xmin+i*h;
    Ua[i] = uana(xg[i]);
    printf("i=%d, Ua(%6.3f)=%6.3f \n", i, xg[i], Ua[i]);
  }  
  printf("\n");
  
  /* initialize k and f */
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
    {
      k[i*n+j] = 0.0;
      ks[i*n+j] = 0.0;
      kd[i*n+j] = 0.0;
    }
    f[i] = 0.0;
    fs[i] = 0.0;
    fd[i] = 0.0;
  }
  /* set boundary conditions */
  k[0*n+0] = 1.0;
  f[0] = uana(0.0);        /* u(x=0.0) = 1.0 */
  k[(n-1)*n+(n-1)] = 1.0;
  f[n-1] = uana(1.0);      /* u(x=1.0) = exp(1) */
  ks[0*n+0] = 1.0;
  fs[0] = uana(0.0);       /* u(x=0.0) = 1.0 */
  ks[(n-1)*n+(n-1)] = 1.0;
  fs[n-1] = uana(1.0);     /* u(x=1.0) = exp(1) */
  kd[0*n+0] = 1.0;
  fd[0] = uana(0.0);       /* u(x=0.0) = 1.0 */
  kd[(n-1)*n+(n-1)] = 1.0;
  fd[n-1] = uana(1.0);     /* u(x=1.0) = exp(1) */

  a = xmin; /* integrate over complete domain */
  b = xmax;
  np = 98;
  printf("calling gauleg a=%g, b=%g, np=%d \n", a, b, np);
  gaulegf(a, b, xx, ww, np); /* xx and ww set up for integrations */
  h = (xmax-xmin)/(double)np; /* for trapezoidal integration */

  /* compute entries in stiffness matrix */
  printf("compute stiffness matrix \n");
  for(i=1; i<n-1; i++)
  {
    for(j=0; j<n; j++)
    {
      /* compute integral for k(i,j)  */
      val = 0.0;
      for(p=1; p<=np; p++)
      {
        val = val + ww[p]*galk(xx[p]);
      }
      printf("Legendre integration=%g, at i=%d, j=%d \n", val, i, j);
      k[i*n+j] = val;

      val = aquad3(galk, xmin, xmax, eps);
      printf("adaptive integration=%g, at i=%d, j=%d \n", val, i, j);
      kd[i*n+j] = val;

      val = 0.0;
      for(p=1; p<np; p++)
      {
        val = val + galk(xmin+p*h)*h;
      }
      val = val + (galk(xmin)+ galk(xmax))*h/2.0;      
      printf("simple   integration=%g, at i=%d, j=%d \n", val, i, j);
      ks[i*n+j] = val;
    } /* end j loop */

    /* compute integral for f(i)  */
    val = 0.0;
    for(p=1; p<=np; p++)
    {
      val = val + ww[p]*galf(xx[p]);
    }
    printf("Legendre integration=%g, f at i=%d \n", val, i);
    f[i] = val;

    val = aquad3(galf, xmin, xmax, eps);
    printf("adaptive integration=%g, f at i=%d \n", val, i);
    fd[i] = val;

    val = 0.0;
    for(p=1; p<np; p++)
    {
      val = val + galf(xmin+p*h)*h;
    }
    val = val + (galf(xmin) + galf(xmax))*h/2.0;    
    printf("simple   integration=%g, f at i=%d \n", val, i);
    fs[i] = val;
  }

  printf("k computed stiffness matrix \n");
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
      printf("%6.3f ", k[i*n+j]);
    printf("\n");
  }
  printf("f computed forcing function \n");
  for(i=0; i<n; i++)
  {
    printf("%6.3f ", f[i]);
  }
  printf("\n");
  simeq(n, k, f, u);

  avgerr = 0.0;
  maxerr = 0.0;
  printf("u computed Galerkin, Ua analytic, error \n");
  for(i=0; i<n; i++)
  {
    err = abs(u[i]-Ua[i]);
    if(err>maxerr) maxerr = err;
    avgerr = avgerr + err;
    printf("u[%d]=%6.3f, Ua=%6.3f, err=%g \n",
            i, u[i], Ua[i], u[i]-Ua[i]);
  }
  printf("                                        maxerr=%g, avgerr=%g \n",
         maxerr, avgerr/(double)n);
  printf("\n");


  printf("kd computed adaptive stiffness matrix \n");
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
      printf("%6.3f ", kd[i*n+j]);
    printf("\n");
  }
  printf("fd computed forcing function \n");
  for(i=0; i<n; i++)
  {
    printf("%6.3f ", fd[i]);
  }
  printf("\n");
  simeq(n, kd, fd, ud);

  avgerr = 0.0;
  maxerr = 0.0;
  printf("ud computed Galerkin, Ua analytic, error \n");
  for(i=0; i<n; i++)
  {
    err = abs(ud[i]-Ua[i]);
    if(err>maxerr) maxerr = err;
    avgerr = avgerr + err;
    printf("ud[%d]=%6.3f, Ua=%6.3f, err=%g \n",
            i, ud[i], Ua[i], ud[i]-Ua[i]);
  }
  printf("                                        maxerr=%g, avgerr=%g \n",
         maxerr, avgerr/(double)n);
  printf("\n");  
  printf("\n");  


  printf("ks computed simple stiffness matrix \n");
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
      printf("%6.3f ", ks[i*n+j]);
    printf("\n");
  }
  printf("fs computed forcing function \n");
  for(i=0; i<n; i++)
  {
    printf("%6.3f ", fs[i]);
  }
  printf("\n");
  simeq(n, ks, fs, us);

  avgerr = 0.0;
  maxerr = 0.0;
  printf("us computed Galerkin, Ua analytic, error \n");
  for(i=0; i<n; i++)
  {
    err = abs(us[i]-Ua[i]);
    if(err>maxerr) maxerr = err;
    avgerr = avgerr + err;
    printf("us[%d]=%6.3f, Ua=%6.3f, err=%g \n",
            i, us[i], Ua[i], us[i]-Ua[i]);
  }
  printf("                                        maxerr=%g, avgerr=%g \n",
         maxerr, avgerr/(double)n);
  printf("\n");  
  printf("\n");  


  printf("Try 20 points \n");
  n = 20;
  h = (xmax-xmin)/(double)(n-1);
  printf("x grid and analytic solution \n");
  for(i=0; i<=n; i++) /* need one extra */
  {
    xg[i] = xmin+i*h;
    Ua[i] = uana(xg[i]);
    printf("i=%d, Ua(%6.3f)=%6.3f \n", i, xg[i], Ua[i]);
  }
  printf("\n");
  
  /* initialize k and f */
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
    {
      k[i*n+j] = 0.0;
      ks[i*n+j] = 0.0;
      kd[i*n+j] = 0.0;
    }
    f[i] = 0.0;
    fs[i] = 0.0;
    fd[i] = 0.0;
  }
  /* set boundary conditions */
  k[0*n+0] = 1.0;
  f[0] = uana(0.0);        /* u(x=0.0) = 1.0 */
  k[(n-1)*n+(n-1)] = 1.0;
  f[n-1] = uana(1.0);      /* u(x=1.0) = exp(1) */
  ks[0*n+0] = 1.0;
  fs[0] = uana(0.0);       /* u(x=0.0) = 1.0 */
  ks[(n-1)*n+(n-1)] = 1.0;
  fs[n-1] = uana(1.0);     /* u(x=1.0) = exp(1) */
  kd[0*n+0] = 1.0;
  fd[0] = uana(0.0);       /* u(x=0.0) = 1.0 */
  kd[(n-1)*n+(n-1)] = 1.0;
  fd[n-1] = uana(1.0);     /* u(x=1.0) = exp(1) */

  np = 64; /* number of points to evaluate integral */
  a = xg[0];
  b = xg[n-1];
  printf("calling gauleg a=%g, b=%g, np=%d \n", a, b, np);
  gaulegf(a, b, xx, ww, np); /* xx and ww set up for integrations */

  /* compute entries in stiffness matrix */
  printf("compute stiffness matrix \n");
  for(i=1; i<n-1; i++)
  {
    for(j=0; j<n; j++)
    {
      /* compute integral for k(i,j)  */
      val = 0.0;
      for(p=1; p<=np; p++)
      {
        val = val + ww[p]*galk(xx[p]);
      }
      printf("Legendre integration=%g, at i=%d, j=%d \n", val, i, j);
      k[i*n+j] = val;

      val = aquad3(galk, a, b, eps);
      printf("adaptive integration=%g, at i=%d, j=%d \n", val, i, j);
      kd[i*n+j] = val;

      val = 0.0;
      h = (b-a)/100.0;
      for(p=1; p<100; p++)
      {
        val = val + galk(a+p*h)*h;
      }
      val = val +(galk(a) + galk(b))*h/2.0;
      printf("simple   integration=%g, at i=%d, j=%d \n", val, i, j);
      ks[i*n+j] = val;
    } /* end j loop */

    /* compute integral for f(i)  */
    val = 0.0;
    for(p=1; p<=np; p++)
    {
      val = val + ww[p]*galf(xx[p]);
    }
    printf("Legendre integration=%g, f at i=%d \n", val, i);
    f[i] = val;

    val = aquad3(galf, a, b, eps);
    printf("adaptive integration=%g, f at i=%d \n", val, i);
    fd[i] = val;

    val = 0.0;
    h = (b-a)/100.0;
    for(p=1; p<100; p++)
    {
      val = val + galf(a+p*h)*h;
    }
    val = val +(galf(a) + galf(b))*h/2.0;
    printf("simple   integration=%g, f at i=%d \n", val, i);
    fs[i] = val;
  }

  printf("k computed stiffness matrix \n");
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
      printf("%6.3f ", k[i*n+j]);
    printf("\n");
  }
  printf("f computed forcing function \n");
  for(i=0; i<n; i++)
  {
    printf("%6.3f ", f[i]);
  }
  printf("\n");
  simeq(n, k, f, u);
  
  avgerr = 0.0;
  maxerr = 0.0;
  printf("u computed Galerkin, Ua analytic, error \n");
  for(i=0; i<n; i++)
  {
    err = abs(u[i]-Ua[i]);
    if(err>maxerr) maxerr = err;
    avgerr = avgerr + err;
    printf("u[%d]=%6.3f, Ua=%6.3f, err=%g \n",
            i, u[i], Ua[i], u[i]-Ua[i]);
  }
  printf("                                        maxerr=%g, avgerr=%g \n",
         maxerr, avgerr/(double)n);
  printf("\n");

  printf("kd computed adaptive stiffness matrix \n");
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
      printf("%6.3f ", kd[i*n+j]);
    printf("\n");
  }
  printf("fd computed forcing function \n");
  for(i=0; i<n; i++)
  {
    printf("%6.3f ", fd[i]);
  }
  printf("\n");
  simeq(n, kd, fd, ud);

  avgerr = 0.0;
  maxerr = 0.0;
  printf("ud computed Galerkin, Ua analytic, error \n");
  for(i=0; i<n; i++)
  {
    err = abs(ud[i]-Ua[i]);
    if(err>maxerr) maxerr = err;
    avgerr = avgerr + err;
    printf("ud[%d]=%6.3f, Ua=%6.3f, err=%g \n",
            i, ud[i], Ua[i], ud[i]-Ua[i]);
  }
  printf("                                        maxerr=%g, avgerr=%g \n",
         maxerr, avgerr/(double)n);
  printf("\n");  

  printf("ks computed simple stiffness matrix \n");
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
      printf("%6.3f ", ks[i*n+j]);
    printf("\n");
  }
  printf("fs computed forcing function \n");
  for(i=0; i<n; i++)
  {
    printf("%6.3f ", fs[i]);
  }
  printf("\n");
  simeq(n, ks, fs, us);

  avgerr = 0.0;
  maxerr = 0.0;
  printf("us computed Galerkin, Ua analytic, error \n");
  for(i=0; i<n; i++)
  {
    err = abs(us[i]-Ua[i]);
    if(err>maxerr) maxerr = err;
    avgerr = avgerr + err;
    printf("us[%d]=%6.3f, Ua=%6.3f, err=%g \n",
            i, us[i], Ua[i], us[i]-Ua[i]);
  }
  printf("                                        maxerr=%g, avgerr=%g \n",
         maxerr, avgerr/(double)n);
  printf("\n");  

  
  return 0;
} /* end main */

/* end fem_checke.c */