/* fem_check4th_la.c  Galerkin FEM
 * solve x^3 u'''' + x^2 u''' + x u'' + u' + u = x^4 + 64x^3 - x^2 - 4x + 1 
 * u(0)=1, u(1)=1  0 < x < 1 
 * analytic solution u(x) = x^4 - x^2 + 1
 * investigate gauss-legendre and adaptive integration
 * solve for Uj numerical approximation U(x_j) of u(x_j)
 * 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  */
static double kd[400]; /* adaptive integration */
static double xg[20];
static double f[20];
static double u[20];
static double fd[20];
static double ud[20];
static double Ua[20];


double F(double x) /* forcing function  */
{
  return x*x*x*x + 64.0*x*x*x - x*x - 4.0*x + 1.0;
}

double uana(double x) /* analytic solution and boundaries */
{
  return x*x*x*x - x*x + 1.0;
}

/* i, j, n, needed by galk or galf, available at call */
int i, j, n;
double galk(double x) /* Galerkin k stiffness function for this problem */
{
  return (x*x*x*phipppp(x,j,n-1,xg) +
          x*x*phippp(x,j,n-1,xg)     +
          x*phipp(x,j,n-1,xg)        +
          phip(x,j,n-1,xg)           +
          phi(x,j,n-1,xg)            )
         *phi(x,i,n-1,xg);
} /* end galk used by aquad3 and other integration */

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

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-5;
  int p, np;

  printf("fem_check4th_la.c running \n");
  printf("Given x^3 u'''' + x^2 u''' + x u'' + u' + u = \n");
  printf("x^4 + 64 x^3 - x^2 - 4 x + 1 \n");
  printf("0<x<1  Boundary u(0) = 1, u(1) = 1 \n");
  printf("Analytic solution u(x) = x^4 - x^2 + 1 \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;
      kd[i*n+j] = 0.0;
    }
    f[i] = 0.0;
    fd[i] = 0.0;
  }
  /* set boundary conditions */
  k[0*n+0] = 1.0;
  f[0] = uana(xmin);        /* u(x=0.0) */
  k[(n-1)*n+(n-1)] = 1.0;
  f[n-1] = uana(xmax);      /* u(x=1.0) */
  kd[0*n+0] = 1.0;
  fd[0] = uana(xmin);       /* u(x=0.0) */
  kd[(n-1)*n+(n-1)] = 1.0;
  fd[n-1] = uana(xmax);     /* u(x=1.0) */

  a = xmin; /* integrate over complete domain */
  b = xmax;
  np = 48;
  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, xmin, xmax, eps);
      printf("adaptive integration=%g, at i=%d, j=%d \n", val, i, j);
      kd[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;
  }

  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("\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");  

  return 0;
} /* end main */

/* end fem_check4th_la.c */