// fem22_la.c  Galerkin FEM
// solve uxx(x,y)+ uyy(x,y) = F(x,y)
// F(x,y) = 2*Pi^2*cos(2*Pi*x)*sin^2(Pi*y) +
//          2*Pi^2*sin^2(Pi*x)*cos(2*Pi*y)
// boundary conditions computed using u(x,y)=0
// analytic solution may be given by u(x) = sin^2(Pi*x)*Sin^2(Pi*y)
// 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))
#define Pi 3.1415926535897932384626433832795028841971
#define nx 11
#define ny 11

static double kg[nx*nx*ny*ny]; // (i*nx+j)*ny*ny + (ii*ny+jj)  
static double xg[nx];
static double yg[ny];
static double fg[nx*ny];
static double ug[nx*ny];
static double Ua[nx*ny];
static double xmin = 0.0; // problem parameters 
static double xmax = 1.0;
static double ymin = 0.0;
static double ymax = 1.0;
static int    debug = 0; // 1 for some, 2 for much
static FILE* Fout;

double F(double x, double y) // forcing function, F(x,y)  
{
  return 2.0*Pi*Pi*cos(2.0*Pi*x)*sin(Pi*y)*sin(Pi*y) +
    2.0*Pi*Pi*sin(Pi*x)*sin(Pi*x)*cos(2.0*Pi*y);
}

double uana(double x, double y) // analytic solution and boundaries 
{
  return sin(Pi*x)*sin(Pi*x)*sin(Pi*y)*sin(Pi*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 ( phipp(x,j,nx-1,xg)* phi(y,jj,ny-1,yg) +
           phi(x,j,nx-1,xg)  * phipp(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("fem22_la.c running \n");
  printf("Given uxx + uyy = 2Pi^2(cos(2Pi*x)sin^2(Pi*y)+sin^2(Pi*x)cos(2Pi*y))\n");
  printf("xmin<=x<=xmax ymin<=y<=ymax Boundaries \n");
  printf("Analytic solution for checking u(x,y) = sin^2(Pi*x)sin^2(Pi*y) \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);
      if(debug>0)
        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);
    if(debug>0)
      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);
    if(debug>0)
      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);
    if(debug>0)
      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]);
	    }
          }
          if(debug>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(debug>1)
        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 if debug > 1 \n");
  printf("f computed forcing function, see above if debug > 1 \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");

  // write result fem22_la_c.dat for gnuplot
  Fout = fopen("fem22_la_c.dat","w");
  for(i=0; i<nx; i++)
  {
    for(ii=0; ii<ny; ii++)
    {
      fprintf(Fout, "%f %f %f \n", xg[i], yg[ii], Ua[i*ny+ii]);
    }
    fprintf(Fout, "\n");
  }
  fclose(Fout);

  return 0;
} // end main 

// end fem22_la.c