/* fem_check44_la.c  Galerkin FEM
 * solve uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) +
 *     4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) +
 *     7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) +
 *     10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) =
 *     F(x,y,z,t)
 *     F(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +
 *                  180 z^2*t + 200 y^2z*t + 44 t^3 +
 *                  110 y^2z^2t + 66 xy + 12t^4 +
 *                  24 z^4 + 36 y^4 + 48 x^4 +
 *                  60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t 
 *
 * boundary conditions computed using u(x,y,z,t)
 * analytic solution is u(x,y,z,t) = t^4 + 2 z^4 + 3 y^4 + 4 x^4 +
 *                                   5 y^2 z^2 t^2 + 6 x y t + 
 *                                   7 x z + 8 t + 9
 * gauss-legendre integration used
 * solve for Uijkl numerical approximation U(x_i,y_j,z_k,t_l)
 * kg * ug = fg, solve simultaneous equations for U
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.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 nx  5
#define ny  5
#define nz  5
#define nt  5

int s(int i, int ii, int iii, int iiii)
{
  return i*ny*nz*nt+ii*nz*nt+iii*nt+iiii;
} /* end s */

int ss(int i, int ii, int iii, int iiii, int j, int jj, int jjj, int jjjj)
{
  return (i*ny*nz*nt+ii*nz*nt+iii*nt+iiii)*nx*ny*nz*nt+
         (j*ny*nz*nt+jj*nz*nt+jjj*nt+jjjj);
} /* end ss */

static double kg[nx*ny*nz*nt*nx*ny*nz*nt];
              /* (i*ny*nz*nt+ii*nz*nt+iii*nt+iiii)*ny*ny*nz*nt+
                 (j*ny*nz*nt+jj*nz*nt+jjj*nt+jjjj) */
static double xg[nx];
static double yg[ny];
static double zg[nz];
static double tg[nt];
static double fg[nx*ny*nz*nt];
static double ug[nx*ny*nz*nt];
static double Ua[nx*ny*nz*nt];
static double xmin = 0.0; /* problem parameters */
static double xmax = 1.0;
static double ymin = 0.0;
static double ymax = 1.0;
static double zmin = 0.0;
static double zmax = 1.0;
static double tmin = 0.0;
static double tmax = 1.0;


static double F(double x, double y, double z, double t)
                                            /* forcing function, F(x,y,z,t) */
{
  return 706.0 + 120.0*t*t + 140.0*y*y*z + 768.0*x + 180.0*z*z*t +
         200.0*y*y*z*t + 44.0*t*t*t + 110.0*y*y*z*z*t + 66.0*x*y +
         12.0*t*t*t*t + 24.0*z*z*z*z + 36.0*y*y*y*y + 48.0*x*x*x*x +
         60.0*y*y*z*z*t*t + 72.0*x*y*t + 84.0*x*z + 96.0*t; 
}

static double uana(double x, double y, double z, double t)
                                        /*analytic solution and boundaries*/
{
  return t*t*t*t + 2.0*z*z*z*z + 3.0*y*y*y*y + 4.0*x*x*x*x +
         5.0*y*y*z*z*t*t + 6.0*x*y*t + 7.0*x*z + 8.0*t + 9.0;
}

/* i, j, ... , jjjj, needed by galk or galf,available call*/
static int i, j, ii, jj, iii, jjj, iiii, jjjj;
static double galk(double x, double y, double z, double t)
{                         /* Galerkin k stiffness function for this problem */
  return
  (    phipppp(x,j,nx-1,xg)*    phi(y,jj,ny-1,yg)*    phi(z,jjj,nz-1,zg)*
                                                      phi(t,jjjj,nt-1,tg)+
   2.0*    phi(x,j,nx-1,xg)*phipppp(y,jj,ny-1,yg)*    phi(z,jjj,nz-1,zg)*
                                                      phi(t,jjjj,nt-1,tg)+
   3.0*    phi(x,j,nx-1,xg)*    phi(y,jj,ny-1,yg)*phipppp(z,jjj,nz-1,zg)*
                                                      phi(t,jjjj,nt-1,tg)+
   4.0*    phi(x,j,nx-1,xg)*    phi(y,jj,ny-1,yg)*    phi(z,jjj,nz-1,zg)*
                                                      phipppp(t,jjjj,nt-1,tg)+
   5.0*   phip(x,j,nx-1,xg)*   phip(y,jj,ny-1,yg)*    phi(z,jjj,nz-1,zg)*
                                                     phip(t,jjjj,nt-1,tg)+
   6.0*    phi(x,j,nx-1,xg)*  phipp(y,jj,ny-1,yg)*  phipp(z,jjj,nz-1,zg)*
                                                      phi(t,jjjj,nt-1,tg)+
   7.0*    phi(x,j,nx-1,xg)*    phi(y,jj,ny-1,yg)*   phip(z,jjj,nz-1,zg)*
                                                    phipp(t,jjjj,nt-1,tg)+
   8.0* phippp(x,j,nx-1,xg)*    phi(y,jj,ny-1,yg)*    phi(z,jjj,nz-1,zg)*
                                                      phi(t,jjjj,nt-1,tg)+
   9.0*    phi(x,j,nx-1,xg)*  phipp(y,jj,ny-1,yg)*    phi(z,jjj,nz-1,zg)*
                                                     phip(t,jjjj,nt-1,tg)+
  10.0*    phi(x,j,nx-1,xg)*    phi(y,jj,ny-1,yg)*   phip(z,jjj,nz-1,zg)*
                                                     phip(t,jjjj,nt-1,tg)+
  11.0*    phi(x,j,nx-1,xg)*    phi(y,jj,ny-1,yg)*    phi(z,jjj,nz-1,zg)*
                                                     phip(t,jjjj,nt-1,tg)+
  12.0*    phi(x,j,nx-1,xg)*    phi(y,jj,ny-1,yg)*    phi(z,jjj,nz-1,zg)*
                                                      phi(t,jjjj,nt-1,tg))*
          (phi(x,i,nx-1,xg)*    phi(y,ii,ny-1,yg)*    phi(z,iii,nz-1,zg)*
                                                      phi(t,iiii,nt-1,tg));
} /* end galk used by integration */

static double galf(double x, double y, double z, double t) 
                                    /* Galerkin f function for this problem */
{
  return F(x,y,z,t)*phi(x,i,nx-1,xg)*phi(y,ii,ny-1,yg)*phi(z,iii,nz-1,zg)*
                                                       phi(t,iiii,nt-1,tg);
} /* end galf used by integration */

int main(int argc, char *argv[])
{
  double x, y, z, t, hx, hy, hz, ht;
  double xx[100], wx[100]; /* for Gauss-Legendre */
  double yy[100], wy[100];
  double zz[100], wz[100];
  double tt[100], wt[100];
  double val, err, avgerr, maxerr;
  int px, npx, py, npy, pz, npz, pt, npt;
  double time_start;
  double time_now;
  double time_last;

  printf("fem_check44_la.c running \n");
  printf("Given uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) + \n");
  printf("    4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) + \n");
  printf("    7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) + \n");
  printf("    10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) = \n");
  printf("  F(x,y,z,t) \n");
  printf("  F(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +   \n");
  printf("               180 z^2*t + 200 y^2z*t + 44 t^3 +    \n");
  printf("               110 y^2z^2t + 66 xy + 12t^4 +        \n");
  printf("               24 z^4 + 36 y^4 + 48 x^4 +           \n");
  printf("               60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t \n"); 
  printf("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries \n");
  printf("Analytic solution u(x,y,z,t) =  t^4 + 2 z^4 + 3 y^4 + \n");
  printf("     4 x^4 + 5 y^2 z^2 t^2 + 6 x y t + 7 x z + 8 t + 9 \n");
  printf("\n");

  printf("xmin=%g, xmax=%g, ymin=%g, ymax=%g \n", xmin, xmax, ymin, ymax);
  printf("zmin=%g, zmax=%g, tmin=%g, tmax=%g \n", zmin, zmax, tmin, tmax);
  printf("nx=%d, ny=%d, nz=%d, nt=%d \n", nx, ny, nz, nt);
  time_start = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time start = %f seconds \n", time_start);
  time_last = time_start;

  printf("x grid and analytic solution at ymin,zmin,tmin \n");
  hx = (xmax-xmin)/(double)(nx-1);
  for(i=0; i<nx; i++)
  {
    xg[i] = xmin+(double)i*hx;
    Ua[i] = uana(xg[i],ymin,zmin,tmin); /* temp */
    printf("i=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
            i, xg[i], ymin, zmin, tmin, Ua[i]);
  } /* i */  
  printf("\n");
  printf("y grid and analytic solution at xmin,zmin,tmin \n");
  hy = (ymax-ymin)/(double)(ny-1);
  for(ii=0; ii<ny; ii++)
  {
    yg[ii] = ymin+(double)ii*hy;
    Ua[ii] = uana(xmin, yg[ii],zmin,tmin); /* temp */
    printf("ii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
            ii, xmin, yg[ii], zmin, tmin, Ua[ii]);
  } /* ii */
  printf("\n");
  printf("z grid and analytic solution at xmin,ymin,tmin \n");
  hz = (zmax-zmin)/(double)(nz-1);
  for(iii=0; iii<nz; iii++)
  {
    zg[iii] = zmin+(double)iii*hz;
    Ua[iii] = uana(xmin, ymin, zg[iii],tmin); /* temp */
    printf("iii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
            iii, xmin, ymin, zg[iii], tmin, Ua[iii]);
  } /* iii */
  printf("\n");
  printf("\n");
  printf("t grid and analytic solution at xmin,ymin,zmin \n");
  ht = (tmax-tmin)/(double)(nt-1);
  for(iiii=0; iiii<nt; iiii++)
  {
    tg[iiii] = tmin+(double)iiii*ht;
    Ua[iiii] = uana(xmin, ymin, zmin, tg[iiii]); /* temp */
    printf("iiii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
            iiii, xmin, ymin, zmin, tg[iiii], Ua[iiii]);
  } /* iiii */
  printf("\n");

  printf("put solution in Ua vector \n");
  /* put solution in Ua vector */
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        for(iiii=0; iiii<nt; iiii++)
        {
          t = tg[iiii];
          Ua[s(i,ii,iii,iiii)] = uana(x,y,z,t);
          /* printf("solution at i=%d,x=%g, ii=%d,y=%g, iii=%d,z=%g, iiii=%d,t=%g is %g \n",
	     i, x, ii, y, iii, z, iiii, t, Ua[s(i,ii,iii,iiii)]); */
        } /* iiii */
      } /* iii */
    } /* ii */
  } /* i */
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");

  printf("initialize k and f \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++)
        {
          for(iii=0; iii<nz; iii++)
          {
            for(jjj=0; jjj<nz; jjj++)
            {
              for(iiii=0; iiii<nt; iiii++)
              {
                for(jjjj=0; jjjj<nt; jjjj++)
                {
                  kg[ss(i,ii,iii,iiii,j,jj,jjj,jjjj)] = 0.0;
	        } /* jjjj */
	      } /* iiii */
	    } /* jjj */
          } /* iii */
        } /* jj */
      } /* ii */
    } /* j */
  } /* i */
  for(i=0; i<nx; i++)
  {
    for(ii=0; ii<ny; ii++)
    {
      for(iii=0; iii<nz; iii++)
      {
        for(iiii=0; iiii<nt; iiii++)
        {
          fg[s(i,ii,iii,iiii)] = 0.0;
        } /* iiii */
      } /* iii */
    } /* ii */
  } /* i */
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");

  printf("set boundary conditions, all sides \n");
  /* set boundary conditions, all sides */
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        for(iiii=0; iiii<nt; iiii++)
        {
          t = tg[iiii];
          kg[ss(i,ii,iii,iiii,i,ii,iii,iiii)] = 1.0;
          fg[s(i,ii,iii,iiii)] = uana(x,y,z,t);
          /* printf("boundary i=%d,x=%g, ii=%d,y=%g, iii=%d,z=%g, iiii=%d,t=%g is %g \n",
             i, x, ii, y, iii, z, iiii, t, fg[s(i,ii,iii,iiii)]); */
        } /* iiii */
      } /* iii */
    } /* ii */
  } /* i */ /* end boundry overkill */


  npx = 12;
  npy = 12;
  npz = 12;
  npt = 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("calling gauleg zmin=%g, zmax=%g, npz=%d \n", zmin, zmax, npz);
  gaulegf(zmin, zmax, zz, wz, npz); /* zz and wz set up for integrations */
  printf("calling gauleg tmin=%g, tmax=%g, npt=%d \n", tmin, tmax, npt);
  gaulegf(tmin, tmax, tt, wt, npt); /* tt and wt 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]);
  printf("zz[1]=%g, zz[2]=%g, wz[1]=%g, wz[2]=%g \n",zz[1],zz[2],wz[1],wz[2]);
  printf("tt[1]=%g, tt[2]=%g, wt[1]=%g, wt[2]=%g \n",tt[1],tt[2],wt[1],wt[2]);
  i=1; j=1; ii=1; jj=1; iii=1; jjj=1; iiii=1; jjjj=1;
  val = galk(xx[2],yy[2],zz[2],tt[2]);
  printf("galk(xx[2],yy[2],zz[2],tt[2])=%g\n",val);
  val = galf(xx[2],yy[2],zz[2],tt[2]);
  printf("galf(xx[2],yy[2],zz[2],tt[2])=%g\n",val);

  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");


  /* 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++)
        {
          for(iii=1; iii<nz-1; iii++)
          {
            for(jjj=0; jjj<nz; jjj++)
            {
              for(iiii=1; iiii<nt-1; iiii++)
              {
                for(jjjj=0; jjjj<nt; jjjj++)
                {
                  /* compute integral for k(x,y,z,t)  */
                  val = 0.0;
                  for(px=1; px<=npx; px++)
                  {
                    for(py=1; py<=npy; py++)
                    {
                      for(pz=1; pz<=npz; pz++)
                      {
                        for(pt=1; pt<=npt; pt++)
                        {
                          val = val + wx[px]*wy[py]*wz[pz]*wt[pt]*
                                galk(xx[px],yy[py],zz[pz],tt[pt]);
                        } /* pt */
                      } /* pz */
	            } /* py */
	          } /* px */
                  kg[ss(i,ii,iii,iiii,j,jj,jjj,jjjj)] = val;
   	        } /* end jjjj loop */
	      } /* end iiii loop */
	    } /* end jjj loop */
	  } /* end iii loop */
        } /* end jj loop */
      } /* end ii loop */
      printf("Legendre=%g, at i=%d,j=%d,ii=%d,jj=%d,iii=%d,jjj=%d,iiii=%d, jjjj=%d \n",
	     val, i, j, ii, jj, iii, jjj, iiii, jjjj); /* just a few */
    } /* end j loop */
  } /* end i loop */
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");

  printf("compute integral for f(i) \n");
  /* compute integral for f(i)  */
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      for(iii=1; iii<nz-1; iii++)
      {
        for(iiii=1; iiii<nt-1; iiii++)
        {
          val = 0.0;
          for(px=1; px<=npx; px++)
          {
            for(py=1; py<=npy; py++)
            {
              for(pz=1; pz<=npz; pz++)
              {
                for(pt=1; pt<=npt; pt++)
                {
                  val = val + wx[px]*wy[py]*wz[pz]*wt[pt]*
                        galf(xx[px],yy[py],zz[pz],tt[pt]);
  	        } /* pt */
	      } /* pz */
            } /* py */
          } /* px */
          fg[s(i,ii,iii,iiii)] = val;
        } /* iiii */
      } /* iii */
    } /* ii */
    printf("Legendre integration=%g, f at i=%d, ii=%d, iii=%d, iiii=%d \n",
           val, i, ii, iii, iiii); /* just a few */
  } /* i */
  printf("k computed stiffness matrix, see above \n");
  printf("f computed forcing function, see above \n");
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");


  printf("solve %d equations in %d unknowns \n", nx*ny*nz*nt, nx*ny*nz*nt);

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

  printf("equations solved \n");
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");


  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++)
    {
      for(iii=0; iii<nz; iii++)
      {
        for(iiii=0; iiii<nt; iiii++)
        {
          err = abs(ug[s(i,ii,iii,iiii)]-Ua[s(i,ii,iii,iiii)]);
          if(err>maxerr) maxerr = err;
          avgerr = avgerr + err;
          printf("ug[%d,%d,%d,%d]=%8.3f, Ua=%8.3f, err=%g \n",
                  i, ii, iii, iiii, ug[s(i,ii,iii,iiii)],
                  Ua[s(i,ii,iii,iiii)], err);
        } /* iiii */
      } /* iii */
    } /* ii */
  } /* i */
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("total CPU time  = %f seconds \n", time_now-time_start);
  printf("\n");


  printf("                                        maxerr=%g, avgerr=%g \n",
         maxerr, avgerr/(double)(nx*ny*nz*nt));
  printf("\n");
  printf("end fem_check44_la.c \n");
  return 0;
} /* end main */

/* end fem_check44_la.c */