/* fem_check_abc_la.c  Galerkin FEM
 * solve a(x,y) uxx(x,y) + b(x,y) uxy(x,y) + c(x,y) uyy(x,y) +
 *       d(x,y) ux(x,y)  + e(x,y) uy(x,y)  + f(x,y) u(x,y) = c(x,y)
 * boundary conditions computed using u(x) in abc.c
 * xmin, xmax, ymin, ymax, nx, ny given in abc.h
 * analytic solution may be given by u(x) in abc.c, ifcheck>0
 * 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"
#include "abc_la.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];


double F(double x, double y) /* forcing function  */
{
  return c(x,y);
}

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

/* 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 (a1(x,y)*phipp(x,j,nx-1,xg)*  phi(y,jj,ny-1,yg) +
          b1(x,y)* phip(x,j,nx-1,xg)* phip(y,jj,ny-1,yg) +
          c1(x,y)*  phi(x,j,nx-1,xg)*phipp(y,jj,ny-1,yg) +
          d1(x,y)* phip(x,j,nx-1,xg)*  phi(y,jj,ny-1,yg) +
          e1(x,y)*  phi(x,j,nx-1,xg)* phip(y,jj,ny-1,yg) +
          f1(x,y)*  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_check_abc_la.c running \n");
  printf("Given a1 uxx + b1 uxy + c1 uyy + d1 ux + e1 uy + f1 u = \n");
  printf("c(x,y) \n");
  printf("xmin<=x<=xmax ymin<=y<=ymax Boundaries \n");
  printf("Analytic solution u(x) =  \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");
  if(ifcheck) /* 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 = 48;
  npy = 48;
  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 */

  /* 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]);
	    }
          }
          if(ifcheck>1)
            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]);
        }
      }
      if(ifcheck>1)
        printf("Legendre integration=%g, f at i=%d, ii=%d \n", val, i, ii);
      fg[i*ny+ii] = val;
    }
  }

  simeq(nx*ny, kg, fg, ug);

  if(ifcheck)
  {
    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));
  }
  else /* no ifcheck */
  {
    printf("solution ug computed by Galerkin method \n");
    for(i=0; i<nx; i++)
    {
      for(ii=0; ii<ny; ii++)
      {
        printf("ug[%d,%d]=%6.3f \n",
                i, ii, ug[i*ny+ii]);
      }
    }
  } /* ifcheck */
  printf("end fem_check_abc_la.c \n");

  return 0;
} /* end main */

/* end fem_check_abc_la.c */