/* fem_check13_la.c  solve linear second order equation
 * x^2 u'' + 2 x u' + 3u =3 -(100/3)*x +108*x^2 -80*x^3
 * with boundary  u(0)=1, u(1)=1  0 < x < 1 
 * try 10 points
 * simple integration
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "laphi.h"
#include "simeq.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 xg[20];
static double u[20];
static double f[20];
static double Ua[20];

/*           forcing function */
#define F(x) (3.0 -(100.0/3.0)*x +108.0*x*x -80.0*x*x*x)

/*           exact solution for checking */
#define U(x) (1.0 -(20.0/3.0)*x +12.0*x*x -(16.0/3.0)*x*x*x)

/*           Galerkin phi substituted for u */
#define gal(x) (x*x*phipp(x,j,n-1,xg)+2.0*x*phip(x,j,n-1,xg)+3.0*phi(x,j,n-1,xg))

int main(int argc, char *argv[])
{
  int i, j, n;
  double xmin = 0.0;
  double xmax = 1.0;
  double x, h, hi;
  double val, err, avgerr, maxerr;
  int p, np;

  printf("fem_check13_la.c running \n");
  printf("Given x^2 u''(x) + 2x u'(x) + 3u(x) = \n");
  printf("  3 -(100/3)*x +108*x^2 -80*x^3 \n");
  printf("  0<x<1  Boundary u(0) = 1, u(1) = 1 \n");
  printf("Analytic solution u(x) = 1 -(20/3)*x +12*x^2 -(16/3)*x^3 \n");

  printf("Try 10 points \n");
  n = 10;
  h = (xmax-xmin)/(double)(n-1);
  printf("F forcing function and exact solution \n");
  for(i=0; i<n; i++)
  {
    xg[i] = xmin+i*h;
    Ua[i] = U(xg[i]);
    printf("i=%d, F(%7.3f)=%7.3f, exact=%7.3f \n",
           i, xg[i], F(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;
    }
    f[i] = 0.0;
  }
  /* set boundary conditions */
  k[0*n+0] = 1.0;
  f[0] = U(xmin);
  k[(n-1)*n+(n-1)] = 1.0;
  f[n-1] = U(xmax);

  np = 20;
  hi = (xmax-xmin)/(double)(np-1); /* for simple 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=0; p<np; p++)
      {
        x = xmin + hi/2.0 + p*hi;
        val = val + gal(x)*phi(x,i,n-1,xg)*hi;
      }
      printf("val=%g, for k at i=%d, j=%d \n", val, i, j);
      k[i*n+j] = val;
    } /* end j loop */

    /* compute integral for f(i)  */
    val = 0.0;
    for(p=0; p<np; p++)
    {
      x = xmin + hi/2.0 + p*hi;
      val = val + F(x)*phi(x,i,n-1,xg)*hi;
    }
    printf("val=%g, for f at i=%d \n", val, i);
    f[i] = val;
  }

  printf("k computed stiffness matrix, see above \n");
  printf("f computed forcing function, see above \n");
  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]=%7.4f, Ua=%7.4f, err=%g \n",
            i, u[i], Ua[i], u[i]-Ua[i]);
  }
  printf("                                        maxerr=%g, avgerr=%g \n",
         maxerr, avgerr/(double)n);
  printf("\n");
  
  return 0;
} /* end main of fem_check13_la.c */