// fem_check47h_la.c  Galerkin FEM
//                    Biharmonic 4th order PDE in seven dimensions
// solve Uxxxx(x,y,z,t,u,v,w) + Uyyyy(x,y,z,t,u,v,w) + Uzzzz(x,y,z,t,u,v,w) + 
//       Utttt(x,y,z,t,u,v,w) + Uuuuu(x,y,z,t,u,v,w) + Uvvvv(x,y,z,t,u,v,w) + 
//       Uwwww(x,y,z,t,u,v,w) + 
//       2 Uxx(x,y,z,t,u,v,w) + 2 Uyy(x,y,z,t,u,v,w) + 2 Uzz(x,y,z,t,u,v,w) +
//       2 Utt(x,y,z,t,u,v,w) + 2 Uuu(x,y,z,t,u,v,w) + 2 Uvv(x,y,z,t,u,v,w) + 
//       2 Uww(x,y,z,t,u,v,w) + 
//       8 U(x,y,z,t,u,v,w)   = F(x,y,z,t,u,v,w)
//
//       F(x,y,z,t,u,v,w) = sin(x+y+z+t+u+v+w)
//
// boundary conditions computed using uana(x,y,z,t,u,v,w)
// analytic solution is U(x,y,z,t,u,v,w) = sin(x+y+z+t+u+v+w)
//
// subscripts in code use x[i], y[i2], z[i3], t[i4], u[i5], v[i6], w[i7]
//
// using Lagrange Phi functions and Gauss-Legendre integration
// solve for U numerical approximation U(x_i,y_i2,z_i3,t_i4,u_i5,v_i6,w_i7)
// kg * uu = 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  3
#define ny  3
#define nz  3
#define nt  3
#define nu  3
#define nv  3
#define nw  3
#define nxyztuvw nx*ny*nz*nt*nu*nv*nt

int s(int i, int i2, int i3, int i4, int i5, int i6, int i7)
{
  // row index to kg, index to uu and fg 
  return i*ny*nz*nt*nu*nv*nw + i2*nz*nt*nu*nv*nw +
         i3*nt*nu*nv*nw + i4*nu*nv*nw + i5*nv*nw + i6*nw + i7;
} // end s

int ss(int i, int i2, int i3, int i4, int i5, int i6, int i7,
       int j, int j2, int j3, int j4, int j5, int j6, int j7)
{
  // row,column index into kg
  return s(i,i2,i3,i4,i5,i6,i7)*nxyztuvw +
         s(j,j2,j3,j4,j5,j6,j7);
} // end ss 

static double kg[nxyztuvw*nxyztuvw];
static double xg[nx];
static double yg[ny];
static double zg[nz];
static double tg[nt];
static double ug[nu];
static double vg[nv];
static double wg[nw];
static double fg[nxyztuvw];
static double uu[nxyztuvw];
static double Ua[nxyztuvw];
static double xmin = 0.0; // problem parameters
static double xmax = 0.5;
static double ymin = 0.0;
static double ymax = 0.5;
static double zmin = 0.0;
static double zmax = 0.5;
static double tmin = 0.0;
static double tmax = 0.5;
static double umin = 0.0;
static double umax = 0.5;
static double vmin = 0.0;
static double vmax = 0.5;
static double wmin = 0.0;
static double wmax = 0.5;


static double F(double x, double y, double z, double t,
                double u, double v, double w)
             // forcing function, F(x,y,z,t,u,v,w) 
{
  return sin(x+y+z+t+u+v+w);
} // end F

static double uana(double x, double y, double z, double t,
                   double u, double v, double w)
                                        //analytic solution and boundaries
{
  return sin(x+y+z+t+u+v+w);
}

// i,...i7, j, ...j7, needed by galk and  galf, available at call
static int i, j, i2, j2, i3, j3, i4, j4, i5, j5, i6, j6, i7, j7;
static double galk(double x, double y, double z, double t,
                   double u, double v, double w)
{                  // Galerkin k stiffness function for this problem 
  return
   (    phipppp(x,j,nx-1,xg)*  phi(y,j2,ny-1,yg)*     phi(z,j3,nz-1,zg)*
        phi(t,j4,nt-1,tg)*     phi(u,j5,nu-1,ug)*     phi(v,j6,nv-1,vg)*
        phi(w,j7,nw-1,wg)+
        phi(x,j,nx-1,xg)*  phipppp(y,j2,ny-1,yg)*     phi(z,j3,nz-1,zg)*
        phi(t,j4,nt-1,tg)*     phi(u,j5,nu-1,ug)*     phi(v,j6,nv-1,vg)*
        phi(w,j7,nw-1,wg)+
        phi(x,j,nx-1,xg)*      phi(y,j2,ny-1,yg)* phipppp(z,j3,nz-1,zg)*
        phi(t,j4,nt-1,tg)*     phi(u,j5,nu-1,ug)*     phi(v,j6,nv-1,vg)*
        phi(w,j7,nw-1,wg)+ 
        phi(x,j,nx-1,xg)*      phi(y,j2,ny-1,yg)*     phi(z,j3,nz-1,zg)*
        phipppp(t,j4,nt-1,tg)* phi(u,j5,nu-1,ug)*     phi(v,j6,nv-1,vg)*
        phi(w,j7,nw-1,wg)+ 
        phi(x,j,nx-1,xg)*      phi(y,j2,ny-1,yg)*     phi(z,j3,nz-1,zg)*
        phi(t,j4,nt-1,tg)* phipppp(u,j5,nu-1,ug)*     phi(v,j6,nv-1,vg)*
        phi(w,j7,nw-1,wg)+ 
        phi(x,j,nx-1,xg)*      phi(y,j2,ny-1,yg)*     phi(z,j3,nz-1,zg)*
        phi(t,j4,nt-1,tg)*     phi(u,j5,nu-1,ug)* phipppp(v,j6,nv-1,vg)*
        phi(w,j7,nw-1,wg)+ 
        phi(x,j,nx-1,xg)*      phi(y,j2,ny-1,yg)*     phi(z,j3,nz-1,zg)*
        phi(t,j4,nt-1,tg)*     phi(u,j5,nu-1,ug)*     phi(v,j6,nv-1,vg)*
        phipppp(w,j7,nw-1,wg)+ 

    2.0*phipp(x,j,nx-1,xg)*    phi(y,j2,ny-1,yg)*     phi(z,j3,nz-1,zg)*
        phi(t,j4,nt-1,tg)*     phi(u,j5,nu-1,ug)*     phi(v,j6,nv-1,vg)*
        phi(w,j7,nw-1,wg)+
    2.0*phi(x,j,nx-1,xg)*    phipp(y,j2,ny-1,yg)*     phi(z,j3,nz-1,zg)*
        phi(t,j4,nt-1,tg)*     phi(u,j5,nu-1,ug)*     phi(v,j6,nv-1,vg)*
        phi(w,j7,nw-1,wg)+
    2.0*phi(x,j,nx-1,xg)*      phi(y,j2,ny-1,yg)*   phipp(z,j3,nz-1,zg)*
        phi(t,j4,nt-1,tg)*     phi(u,j5,nu-1,ug)*     phi(v,j6,nv-1,vg)*
        phi(w,j7,nw-1,wg)+ 
    2.0*phi(x,j,nx-1,xg)*      phi(y,j2,ny-1,yg)*     phi(z,j3,nz-1,zg)*
        phipp(t,j4,nt-1,tg)*   phi(u,j5,nu-1,ug)*     phi(v,j6,nv-1,vg)*
        phi(w,j7,nw-1,wg)+ 
    2.0*phi(x,j,nx-1,xg)*      phi(y,j2,ny-1,yg)*     phi(z,j3,nz-1,zg)*
        phi(t,j4,nt-1,tg)*   phipp(u,j5,nu-1,ug)*     phi(v,j6,nv-1,vg)*
        phi(w,j7,nw-1,wg)+ 
    2.0*phi(x,j,nx-1,xg)*      phi(y,j2,ny-1,yg)*     phi(z,j3,nz-1,zg)*
        phi(t,j4,nt-1,tg)*     phi(u,j5,nu-1,ug)*   phipp(v,j6,nv-1,vg)*
        phi(w,j7,nw-1,wg)+ 
    2.0*phi(x,j,nx-1,xg)*      phi(y,j2,ny-1,yg)*     phi(z,j3,nz-1,zg)*
        phi(t,j4,nt-1,tg)*     phi(u,j5,nu-1,ug)*     phi(v,j6,nv-1,vg)*
        phipppp(w,j7,nw-1,wg)+ 

    8.0*phi(x,j,nx-1,xg)*      phi(y,j2,ny-1,yg)*      phi(z,j3,nz-1,zg)*
        phi(t,j4,nt-1,tg)*     phi(u,j5,nu-1,ug)*      phi(v,j6,nv-1,vg)* 
        phi(w,j7,nw-1,wg)                                                 )*
       (phi(x,i,nx-1,xg)*      phi(y,i2,ny-1,yg)*      phi(z,i3,nz-1,zg)*
        phi(t,i4,nt-1,tg)*     phi(u,i5,nu-1,ug)*      phi(v,i6,nv-1,vg)* 
        phi(w,i7,nw-1,wg) );
} // end galk used by integration

static double galf(double x, double y, double z, double t,
                   double u, double v, double w) 
              // Galerkin f function for this problem 
{
  return F(x,y,z,t,u,v,w)*
         phi(x,i ,nx-1,xg)*phi(y,i2,ny-1,yg)*phi(z,i3,nz-1,zg)*
         phi(t,i4,nt-1,tg)*phi(u,i5,nu-1,ug)*phi(v,i6,nv-1,vg)*
         phi(w,i7,nw-1,wg);
} // end galf used by integration

int main(int argc, char *argv[])
{
  double x, y, z, t, u, v, w, hx, hy, hz, ht, hu, hv, hw;
  double xx[100], wx[100]; // for Gauss-Legendre  coord,weight 
  double yy[100], wy[100];
  double zz[100], wz[100];
  double tt[100], wt[100];
  double uuu[100], wu[100]; // name conflict
  double vv[100], wv[100];
  double www[100], ww[100]; // name conflict
  double tmp, val, err, avgerr, maxerr;
  int px, npx, py, npy, pz, npz, pt, npt, pu, npu, pv, npv, pw, npw;
  double time_start;
  double time_now;
  double time_last;

  printf("fem_check47h_la.c running \n");
  printf("Given Uxxxx(x,y,z,t,u,v,w) + Uyyyy(x,y,z,t,u,v,w) + Uzzzz(x,y,z,t,u,v,w) + \n");
  printf("      Utttt(x,y,z,t,u,v,w) + Uuuuu(x,y,z,t,u,v,w) + Uvvvv(x,y,z,t,u,v,w) + \n");
  printf("      Uwwww(x,y,z,t,u,v,w) + \n");
  printf("      2 Uxx(x,y,z,t,u,v,w) + 2 Uyy(x,y,z,t,u,v,w) + 2 Uzz(x,y,z,t,u,v,w) + \n");
  printf("      2 Utt(x,y,z,t,u,v,w) + 2 Uuu(x,y,z,t,u,v,w) + 2 Uvv(x,y,z,t,u,v,w) + \n");
  printf("      2 Uww(x,y,z,t,u,v,w) + \n");
  printf("      8 U(x,y,z,t,u,v,w)   = F(x,y,z,t,u,v,w) \n");
  printf(" F(x,y,z,t,u,v,w) = sin(x+y+z+t+u+v+w) \n");
  printf(" \n");
  printf("Biharmonic 4th order PDE in seven dimensions \n");
  printf("Solve by  Galerkin Finite Element Method \n");
  printf("  using Lagrange Phi functions and Gauss-Legendre integration \n");
  printf("\n");
  printf("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries \n");
  printf("tmin<=t<=tmax umin<=u<=umax vmin<=v<=vmax Boundaries \n");
  printf("wmin<=w<=wmax                             Boundaries \n");
  printf("Analytic solution U(x,y,z,t,u,v,w) = sin(x+y+z+t+u+v+w) \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("umin=%g, umax=%g, vmin=%g, vmax=%g \n", umin, umax, vmin, vmax);
  printf("wmin=%g, wmax=%g \n", wmin, wmax);
  printf("nx=%d, ny=%d, nz=%d, nt=%d \n", nx, ny, nz, nt);
  printf("nu=%d, nv=%d, nw=%d \n", nu, nv, nw);

  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,umin,vmin,wmin \n");
  hx = (xmax-xmin)/(double)(nx-1);
  for(i=0; i<nx; i++)
  {
    xg[i] = xmin+(double)i*hx;
    tmp = uana(xg[i],ymin,zmin,tmin,umin,vmin,wmin);
    printf("i=%d, U(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3F,%6.3f)=%7.4f \n",
	   i, xg[i],ymin,zmin,tmin,umin,vmin,wmin,tmp);
  } // i   
  printf("\n");
  printf("y grid and analytic solution at xmax,zmin,tmin,umin,vmin,wmin \n");
  hy = (ymax-ymin)/(double)(ny-1);
  for(i2=0; i2<ny; i2++)
  {
    yg[i2] = ymin+(double)i2*hy;
    tmp = uana(xmax,yg[i2],zmin,tmin,umin,vmin,wmin);
    printf("i=%d, U(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3F,%6.3f)=%7.4f \n",
	   i2, xmax,yg[i2],zmin,tmin,umin,vmin,wmin,tmp);
  } // i2
  printf("\n");
  printf("z grid and analytic solution at xmax,ymax,tmin,umin,vmin,wmin \n");
  hz = (zmax-zmin)/(double)(nz-1);
  for(i3=0; i3<nz; i3++)
  {
    zg[i3] = zmin+(double)i3*hz;
    tmp = uana(xmax,ymax,zg[i3],tmin,umin,vmin,wmin);
    printf("i=%d, U(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3F,%6.3f)=%7.4f \n",
	   i3, xmax,ymax,zg[i3],tmin,umin,vmin,wmin,tmp);
  } // i3
  printf("\n");
  printf("t grid and analytic solution at xmax,ymax,zmax,umin,vmin,wmin \n");
  ht = (tmax-tmin)/(double)(nt-1);
  for(i4=0; i4<nt; i4++)
  {
    tg[i4] = tmin+(double)i4*ht;
    tmp = uana(xmax,ymax,zmax,tg[i4],umin,vmin,wmin);
    printf("i=%d, U(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3F,%6.3f)=%7.4f \n",
	   i4, xmax,ymax,zmax,tg[i4],umin,vmin,wmin,tmp);
  } // i4 
  printf("\n");
  printf("u grid and analytic solution at xmax,ymax,zmax,tmax,vmin,wmin \n");
  hu = (umax-umin)/(double)(nu-1);
  for(i5=0; i5<nu; i5++)
  {
    ug[i5] = umin+(double)i5*hu;
    tmp = uana(xmax,ymax,zmax,tmax,ug[i5],vmin,wmin);
    printf("i=%d, U(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3F,%6.3f)=%7.4f \n",
	   i5, xmax,ymax,zmax,tmax,ug[i5],vmin,wmin,tmp);
  } // i5 
  printf("\n");
  printf("v grid and analytic solution at xmax,ymax,zmax,tmax,umax,wmin \n");
  hv = (vmax-vmin)/(double)(nv-1);
  for(i6=0; i6<nv; i6++)
  {
    vg[i6] = vmin+(double)i6*hv;
    tmp = uana(xmax,ymax,zmax,tmax,umax,vg[i6],wmin);
    printf("i=%d, U(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3F,%6.3f)=%7.4f \n",
	   i6, xmax,ymax,zmax,tmax,umax,ug[i6],wmin,tmp);
  } // i6 
  printf("\n");
  printf("w grid and analytic solution at xmax,ymax,zmax,tmax,umax,vmax \n");
  hw = (wmax-wmin)/(double)(nw-1);
  for(i7=0; i7<nw; i7++)
  {
    wg[i7] = wmin+(double)i7*hw;
    tmp = uana(xmax,ymax,zmax,tmax,umax,vmax,wg[i7]);
    printf("i=%d, U(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3F,%6.3f)=%7.4f \n",
	   i7, xmax,ymax,zmax,tmax,umax,vmax,wg[i7],tmp);
  } // i7 
  printf("\n");

  printf("put solution in Ua vector \n");
  // put solution in Ua vector 
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(i2=0; i2<ny; i2++)
    {
      y = yg[i2];
      for(i3=0; i3<nz; i3++)
      {
        z = zg[i3];
        for(i4=0; i4<nt; i4++)
        {
          t = tg[i4];
          for(i5=0; i5<nu; i5++)
          {
            u = ug[i5];
            for(i6=0; i6<nv; i6++)
            {
              v = vg[i6];
              for(i7=0; i7<nw; i7++)
              {
                w = wg[i7];
         Ua[s(i,i2,i3,i4,i5,i6,i7)] = uana(x,y,z,t,u,v,w);
         // printf("solution at i=%1d,x=%g, i2=%1d,y=%g, i3=%1d,z=%g,\n
         // i4=%1d,t=%g, i5=%1d,u=%g, i6=%1d,v=%g, 
         // i7=%1d,w=%g, is %g \n",
         // i, x, i2, y, i3, z, i4, t, i5, u, i6, v, i7, w,Ua[s(i,i2,i3,i4)]); 
              } // i7 
            } // i6 
          } // i5 
        } // i4 
      } // i3 
    } // i2 
  } // 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(i2=0; i2<ny; i2++)
      {
        for(j2=0; j2<ny; j2++)
        {
          for(i3=0; i3<nz; i3++)
          {
            for(j3=0; j3<nz; j3++)
            {
              for(i4=0; i4<nt; i4++)
              {
                for(j4=0; j4<nt; j4++)
                {
                  for(i5=0; i5<nu; i5++)
                  {
                    for(j5=0; j5<nu; j5++)
                    {
                      for(i6=0; i6<nv; i6++)
                      {
                        for(j6=0; j6<nv; j6++)
                        {
                          for(i7=0; i7<nw; i7++)
                          {
                            for(j7=0; j7<nw; j7++)
                            {
	         kg[ss(i,i2,i3,i4,i5,i6,i7,j,j2,j3,j4,j5,j6,j7)] = 0.0;
                            } // i7 
                          } // i7 
                        } // j6 
                      } // i6 
                    } // j5 
                  } // i5 
                } // j4 
	      } // i4 
	    } // j3 
          } // i3 
        } // j2 
      } // i2 
    } // j 
  } // i 
  for(i=0; i<nx; i++)
  {
    for(i2=0; i2<ny; i2++)
    {
      for(i3=0; i3<nz; i3++)
      {
        for(i4=0; i4<nt; i4++)
        {
          for(i5=0; i5<nu; i5++)
          {
            for(i6=0; i6<nv; i6++)
            {
              for(i7=0; i7<nw; i7++)
              {
		fg[s(i,i2,i3,i4,i5,i6,i7)] = 0.0;
              } // i7 
            } // i6 
          } // i5 
        } // i4 
      } // i3 
    } // i2 
  } // 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 kg diagonal and fg \n");
  // set boundary conditions, all faces 
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(i2=0; i2<ny; i2++)
    {
      y = yg[i2];
      for(i3=0; i3<nz; i3++)
      {
        z = zg[i3];
        for(i4=0; i4<nt; i4++)
        {
          t = tg[i4];
          for(i5=0; i5<nu; i5++)
          {
            u = ug[i5];
            for(i6=0; i6<nv; i6++)
            {
              v = vg[i6];
              for(i7=0; i7<nw; i7++)
              {
                w = wg[i7];
          kg[ss(i,i2,i3,i4,i5,i6,i7,i,i2,i3,i4,i5,i6,i7)] = 1.0;
          val = uana(x,y,z,t,u,v,w);
          fg[s(i,i2,i3,i4,i5,i6,i7)] = 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");

  npx = 8;
  npy = 8;
  npz = 8;
  npt = 8;
  npu = 8;
  npv = 8;
  npw = 8;
  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("calling gauleg umin=%g, umax=%g, npu=%d \n", umin, umax, npu);
  gaulegf(umin, umax, uuu, wu, npu); // uuu and wu set up for integrations 
  printf("calling gauleg vmin=%g, vmax=%g, npv=%d \n", vmin, vmax, npv);
  gaulegf(vmin, vmax, vv, wv, npv); // vv and wv set up for integrations 
  printf("calling gauleg wmin=%g, wmax=%g, npw=%d \n", wmin, wmax, npw);
  gaulegf(wmin, wmax, www, ww, npw); // www and ww 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]);
  printf("uuu[1]=%g, uuu[2]=%g, wu[1]=%g, wu[2]=%g \n",
          uuu[1], uuu[2], wu[1], wu[2]);
  printf("vv[1]=%g, vv[2]=%g, wv[1]=%g, wv[2]=%g \n",vv[1],vv[2],wv[1],wv[2]);
  printf("www[1]=%g, www[2]=%g, ww[1]=%g, ww[2]=%g \n",
          www[1], www[2], ww[1], ww[2]);

  i=1; j=1; i2=1; j2=1; i3=1; j3=1; i4=1; j4=1;
  i5=1; j5=1; i6=1; j6=1; i7=1; j7=1; 
  val = galk(xx[2],yy[2],zz[2],tt[2],uuu[2],vv[2],www[2]);
  printf("galk(xx[2],yy[2],zz[2],tt[2],uuu[2],vv[2],www[2])=%g\n",val);
  val = galf(xx[2],yy[2],zz[2],tt[2],uuu[2],vv[2],www[2]);
  printf("galf(xx[2],yy[2],zz[2],tt[2],uuu[2],vv[2],www[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(i2=1; i2<ny-1; i2++)
      {
        for(j2=0; j2<ny; j2++)
        {
          for(i3=1; i3<nz-1; i3++)
          {
            for(j3=0; j3<nz; j3++)
            {
              for(i4=1; i4<nt-1; i4++)
              {
                for(j4=0; j4<nt; j4++)
                {
                  for(i5=1; i5<nu-1; i5++)
                  {
                    for(j5=0; j5<nu; j5++)
                    {
                      for(i6=1; i6<nv-1; i6++)
                      {
                        for(j6=0; j6<nv; j6++)
                        {
                          for(i7=1; i7<nw-1; i7++)
                          {
                            for(j7=0; j7<nw; j7++)
                            {
                              // compute integral for k(x,y,z,t,u,v,w)  
                              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++)
                                    {
                                      for(pu=1; pu<=npu; pu++)
                                      {
                                        for(pv=1; pv<=npv; pv++)
                                        {
                                          for(pw=1; pw<=npw; pw++)
                                          {
                val = val + wx[px]*wy[py]*wz[pz]*wt[pt]*wu[pu]*wv[pv]*ww[pw]*
		  galk(xx[px],yy[py],zz[pz],tt[pt],uuu[pu],vv[pv],www[pw]);
                                          } // pw 
                                        } // pv 
                                      } // pu 
                                    } // pt 
                                  } // pz 
	                        } // py 
	                      } // px 
		  kg[ss(i,i2,i3,i4,i5,i6,i7,j,j2,j3,j4,j5,j6,j7)] = val;
   	                    } // end j7 loop 
	                  } // end i7 loop 
  	                } // end j6 loop 
	              } // end i6 loop 
  	            } // end j5 loop 
	          } // end i5 loop 
  	        } // end j4 loop 
	      } // end i4 loop 
	    } // end j3 loop 
	  } // end i3 loop 
        } // end j2 loop 
      } // end i2 loop 
      printf("Legendre=%g, at i=%d,j=%d,i2=%d,j2=%d,i3=%d,j3=%d,i4=%d, j4=%d \n",
	     val, i, j, i2, j2, i3, j3, i4, j4); // 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(i2=1; i2<ny-1; i2++)
    {
      for(i3=1; i3<nz-1; i3++)
      {
        for(i4=1; i4<nt-1; i4++)
        {
          for(i5=1; i5<nu-1; i5++)
          {
            for(i6=1; i6<nv-1; i6++)
            {
              for(i7=1; i7<nw-1; i7++)
              {
                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++)
                      {
                        for(pu=1; pu<=npu; pu++)
                        {
                          for(pv=1; pv<=npv; pv++)
                          {
                            for(pw=1; pw<=npw; pw++)
                            {
                              val = val + 
                              wx[px]*wy[py]*wz[pz]*wt[pt]*wu[pu]*wv[pv]*ww[pw]*
				galf(xx[px],yy[py],zz[pz],tt[pt],uuu[pu],vv[pv],www[pw]);
  	                    } // pw 
  	                  } // pv 
  	                } // pu 
  	              } // pt 
	            } // pz 
                  } // py 
                } // px 
                fg[s(i,i2,i3,i4,i5,i6,i7)] = val;
              } // i7 
            } // i6 
          } // i5 
        } // i4 
      } // i3 
    } // i2 
    printf("Legendre integration=%g, f at i=%1d, i2=%1d, i3=%1d, i4=%1d, i5=%1d, i6=%1d, i7=%1d \n",
           val, i, i2, i3, i4, i5, i6, i7); // 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", nxyztuvw, nxyztuvw);

  simeq(nxyztuvw, kg, fg, uu);

  printf("%d equations solved \n", nxyztuvw);
  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("uu computed Galerkin, Ua analytic, error \n");
  for(i=0; i<nx; i++)
  {
    for(i2=0; i2<ny; i2++)
    {
      for(i3=0; i3<nz; i3++)
      {
        for(i4=0; i4<nt; i4++)
        {
          for(i5=0; i5<nu; i5++)
          {
            for(i6=0; i6<nv; i6++)
            {
              for(i7=0; i7<nw; i7++)
              {
                err = abs(uu[s(i,i2,i3,i4,i5,i6,i7)]-
                          Ua[s(i,i2,i3,i4,i5,i6,i7)]);
                if(err>maxerr) maxerr = err;
                avgerr = avgerr + err;
       printf("uu[%1d,%1d,%1d,%1d,%1d,%1d,%1d]=%8.3f, Ua=%8.3f, err=%g \n",
		  i, i2, i3, i4, i5, i6, i7, uu[s(i,i2,i3,i4,i5,i6,i7)],
		  Ua[s(i,i2,i3,i4,i5,i6,i7)], err);
              } // i7 
            } // i6
          } // i5
        } // i4 
      } // i3 
    } // i2 
  } // i 
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("total CPU time  = %f seconds \n", time_now-time_start);
  time_last = time_now;
  printf("\n");


  printf("                                        maxerr=%g, avgerr=%g \n",
         maxerr, avgerr/(double)(nxyztuvw));
  printf("\n ");
  printf("end fem_check47h_la.c \n");
  return 0;
} // end main

// end fem_check47h_la.c