/* fem_check22_la.c  Galerkin FEM
 * solve uxx(x,y) + 2 uxy(x,y) + 3 uyy(x,y) +
 *       4 ux(x,y)  + 5 uy(x,y)  + 6 u(x,y) = c(x,y)
 * boundary conditions computed using u(x)
 * analytic solution may be given by u(x) = x^3 + y^3 + xy + 1
 * gauss-legendre integration used
 * 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"

#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 kg[160000]; /* (i*nx+j)*ny*ny + (ii*ny+jj)  */
static double xg[400];
static double yg[400];
static double fg[400];
static double ug[400];
static double Ua[400];
static double xmin = 0.0; /* problem parameters */
static double xmax = 1.0;
static double ymin = 0.0;
static double ymax = 1.0;
static int nx = 4;
static int ny = 4;


double F(double x, double y) /* forcing function, c(x,y)  */
{
  return 6.0*x*x*x+6.0*y*y*y+12.0*x*x+15.0*y*y+6.0*x*y+11.0*x+22.0*y+8.0;
}

double uana(double x, double y) /* analytic solution and boundaries */
{
  return x*x*x+y*y*y+x*y+1;
}

/* i, j, nx, ii, jj, ny, needed by galk or galf, available at call */
int i, j, ii, jj;
double galk(double x, double y)
{                         /* Galerkin k stiffness function for this problem */
  return (      phipp(x,j,nx-1,xg)*  phi(y,jj,ny-1,yg) +
           2.0 * phip(x,j,nx-1,xg)* phip(y,jj,ny-1,yg) +
           3.0 *  phi(x,j,nx-1,xg)*phipp(y,jj,ny-1,yg) +
           4.0 * phip(x,j,nx-1,xg)*  phi(y,jj,ny-1,yg) +
           5.0 *  phi(x,j,nx-1,xg)* phip(y,jj,ny-1,yg) +
           6.0 *  phi(x,j,nx-1,xg)*  phi(y,jj,ny-1,yg) )*
                  phi(x,i,nx-1,xg)*  phi(y,ii,ny-1,yg);
} /* end galk used by aquad3 and other 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 aquad and other integration */

int main(int argc, char *argv[])
{
  /* i, j and n above, used by aquad3 */
  double x, y, hx, hy;
  double a, b, sum;
  double xx[100], wx[100]; /* for Gauss-Legendre */
  double yy[100], wy[100]; /* for Gauss-Legendre */
  double val, err, avgerr, maxerr;
  int px, npx, py, npy;

  printf("fem_check22_la.c running \n");
  printf("Given uxx + 2 uxy + 3 uyy + 4 ux + 5 uy + 6 u = \n");
  printf("6x^3+6y^3+12x^2+15y^2+6xy+11x+22y+8 \n");
  printf("xmin<=x<=xmax ymin<=y<=ymax Boundaries \n");
  printf("Analytic solution u(x,y) =  x^3 + y^3 +xy + 1 \n");

  printf("xmin=%g, xmax=%g, ymin=%g, ymax=%g \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++)
  {
    xg[i] = xmin+i*hx;
    Ua[i] = uana(xg[i],ymin); /* temp */
    printf("i=%d, Ua(%6.3f)=%6.3f \n", i, xg[i], Ua[i]);
  }  
  printf("\n");
  printf("y grid and analytic solution at xmin\n");
  hy = (ymax-ymin)/(double)(ny-1);
  for(ii=0; ii<ny; ii++)
  {
    yg[ii] = ymin+ii*hy;
    Ua[ii] = uana(xmin, yg[ii]); /* temp */
    printf("ii=%d, Ua(%6.3f)=%6.3f \n", ii, yg[ii], Ua[ii]);
  }  
  printf("\n");
  /* put solution in Ua vector */
  for(i=0; i<nx; i++)
  {
    x = xmin+i*hx;
    for(ii=0; ii<ny; ii++)
    {
      y = ymin+ii*hy;
      Ua[i*ny+ii] = uana(x,y);
      printf("solution at i=%d,x=%g, ii=%d,y=%g is %g \n",
             i, x, ii, y, Ua[i*ny+ii]);
    }
  }
  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)*nx*ny+(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)*nx*ny+(i*ny+ii)] = 1.0;
    fg[i*ny+ii] = uana(x,y);
    printf("boundary i=%d,x=%g, ii=%d,y=%g is %g \n",
           i, x, ii, y, fg[i*ny+ii]);
    ii=ny-1;
    y = yg[ii];
    kg[(i*ny+ii)*nx*ny+(i*ny+ii)] = 1.0;
    fg[i*ny+ii] = uana(x,y);
    printf("boundary i=%d,x=%g, ii=%d,y=%g is %g \n",
           i, x, ii, y, fg[i*ny+ii]);
  }
  for(ii=0; ii<ny; ii++)
  {
    y = yg[ii];
    i=0;
    x = xg[i];
    kg[(i*ny+ii)*nx*ny+(i*ny+ii)] = 1.0;
    fg[i*ny+ii] = uana(x,y);
    printf("boundary i=%d,x=%g, ii=%d,y=%g is %g \n",
           i, x, ii, y, fg[i*ny+ii]);
    i=nx-1;
    x = xg[i];
    kg[(i*ny+ii)*nx*ny+(i*ny+ii)] = 1.0;
    fg[i*ny+ii] = uana(x,y);
    printf("boundary i=%d,x=%g, ii=%d,y=%g is %g \n",
           i, x, ii, y, fg[i*ny+ii]);
  }


  npx = 12;
  npy = 12;
  printf("calling gauleg xmin=%g, xmax=%g, npx=%d \n", xmin, xmax, npx);
  gaulegf(xmin, xmax, xx, wx, npx); /* xx and wx set up for integrations */
  printf("calling gauleg ymin=%g, ymax=%g, npy=%d \n", ymin, ymax, npy);
  gaulegf(ymin, ymax, yy, wy, npy); /* yy and wy set up for integrations */
  printf("xx[1]=%g, xx[2]=%g, wx[1]=%g, wx[2]=%g \n",xx[1],xx[2],wx[1],wx[2]);
  printf("yy[1]=%g, yy[2]=%g, wy[1]=%g, wy[2]=%g \n",yy[1],yy[2],wy[1],wy[2]);
  i=1; j=1; ii=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);

  /* 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]);
	    }
          }
          printf("Legendre integration=%g, at i=%d, j=%d, ii=%d, jj=%d \n",
                 val, i, j, ii, jj);
          kg[(i*ny+ii)*nx*ny+(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]);
        }
      }
      printf("Legendre integration=%g, f at i=%d, ii=%d \n", val, i, ii);
      fg[i*ny+ii] = val;
    }
  }
  printf("k computed stiffness matrix, see above \n");
  printf("f computed forcing function, see above \n");
  simeq(nx*ny, kg, fg, ug);

  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]=%6.3f, Ua=%6.3f, err=%g \n",
              i, ii, ug[i*ny+ii], Ua[i*ny+ii], ug[i*ny+ii]-Ua[i*ny+ii]);
    }
  }
  printf("                                        maxerr=%g, avgerr=%g \n",
         maxerr, avgerr/(double)(nx*ny));
  printf("\n");

  return 0;
} /* end main */

/* end fem_check22_la.c */