/* fem_check33_la.c  Galerkin FEM
 * solve uxxx(x,y,z) + 2 uyyy(x,y,z) + 3 uzzz(x,y,z) +
 *     4 uxxy(x,y,z) + 5 uxxz(x,y,z) + 6 uyyx(x,y,z) +
 *     7 uyyz(x,y,z) + 8 uzzx(x,y,z) + 9 uzzy(x,y,z) +
 *     10 uxx(x,y,z) + 11 uyy(x,y,z) + 12 uzz(x,y,z) +
 *     13 uxy(x,y,z) + 14 uxz(x,y,z) + 15 uyz(x,y,z) +
 *     16  ux(x,y,z) + 17  uy(x,y,z) + 18  uz(x,y,z) +
 *     19   u(x,y,z) = F(x,y,z)
 * F(x,y,z) =
 *    406.0 + 208.0*x*y*z*z + 240.0*x*x*y*z + 88.0*x*x*z*z +
 *    136.0*x*x*y*z*z + 128.0*x*y*y*z*z + 217.0*x + 302.0*y +
 *    80.0*y*y*z*z + 144.0*x*x*y*y*z + 96.0*x*x*y*y +
 *    224.0*x*y*y*z + 76.0*x*x*y*y*z*z + 144.0*x*x*y +
 *    80.0*y*y*z + 64.0*y*z*z + 133.0*y*z + 128.0*x*y*y +
 *    112.0*x*x*z + 96.0*x*z*z + 114.0*x*z + 95.0*x*y + 431.0*z +
 *    19.0*x*x*x*x + 38.0*y*y*y*y + 57.0*z*z*z*z + 64.0*x*x*x +
 *    136.0*y*y*y + 216.0*z*z*z + 120.0*x*x + 264.0*y*y + 432.0*z*z; 
 *
 * boundary conditions computed using u(x,y,z)
 * analytic solution may be given by u(x,y,z) = x^4 + 2 y^4 + 3 z^4 +
 *                   4 x^2 y^2 z^2 + 5 x y + 6 x z + 7 y z + 8
 * gauss-legendre integration used
 * solve for Uijk numerical approximation U(x_i,y_j,z_k) of u(x_i,y_j,z_k)
 * 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

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

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

static double kg[nx*ny*nz*nx*ny*nz];
         /* (i*ny*nz+ii*nz+iii)*ny*ny*nz + (j*ny*nz+jj*nz+jjj)  */
static double xg[nx];
static double yg[ny];
static double zg[nz];
static double fg[nx*ny*nz];
static double ug[nx*ny*nz];
static double Ua[nx*ny*nz];
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 F(double x, double y, double z) /* forcing function, F(x,y,z) */
{
  return 406.0 + 208.0*x*y*z*z + 240.0*x*x*y*z + 88.0*x*x*z*z +
         136.0*x*x*y*z*z + 128.0*x*y*y*z*z + 217.0*x + 302.0*y +
         80.0*y*y*z*z + 144.0*x*x*y*y*z + 96.0*x*x*y*y +
         224.0*x*y*y*z + 76.0*x*x*y*y*z*z + 144.0*x*x*y +
         80.0*y*y*z + 64.0*y*z*z + 133.0*y*z + 128.0*x*y*y +
         112.0*x*x*z + 96.0*x*z*z + 114.0*x*z + 95.0*x*y + 431.0*z +
         19.0*x*x*x*x + 38.0*y*y*y*y + 57.0*z*z*z*z + 64.0*x*x*x +
         136.0*y*y*y + 216.0*z*z*z + 120.0*x*x + 264.0*y*y + 432.0*z*z; 
}

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

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

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

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

  printf("fem_check33_la.c running \n");
  printf("Given uxxx(x,y,z) + 2 uyyy(x,y,z) + 3 uzzz(x,y,z) + \n");
  printf(" 4 uxxy(x,y,z) + 5 uxxz(x,y,z) + 6 uyyx(x,y,z) + \n");
  printf(" 7 uyyz(x,y,z) + 8 uzzx(x,y,z) + 9 uzzy(x,y,z) + \n");
  printf(" 10 uxx(x,y,z) + 11 uyy(x,y,z) + 12 uzz(x,y,z) + \n");
  printf(" 13 uxy(x,y,z) + 14 uxz(x,y,z) + 15 uyz(x,y,z) + \n");
  printf(" 16  ux(x,y,z) + 17  uy(x,y,z) + 18  uz(x,y,z) + \n");
  printf(" 19   u(x,y,z) = F(x,y,z) \n");
  printf(" F(x,y,z) = \n");
  printf(" 406.0 + 208.0*x*y*z*z + 240.0*x*x*y*z + 88.0*x*x*z*z + \n");
  printf(" 136.0*x*x*y*z*z + 128.0*x*y*y*z*z + 217.0*x + 302.0*y + \n");
  printf(" 80.0*y*y*z*z + 144.0*x*x*y*y*z + 96.0*x*x*y*y + \n");
  printf(" 224.0*x*y*y*z + 76.0*x*x*y*y*z*z + 144.0*x*x*y + \n");
  printf(" 80.0*y*y*z + 64.0*y*z*z + 133.0*y*z + 128.0*x*y*y + \n");
  printf(" 112.0*x*x*z + 96.0*x*z*z + 114.0*x*z + 95.0*x*y + 431.0*z + \n");
  printf(" 19.0*x*x*x*x + 38.0*y*y*y*y + 57.0*z*z*z*z + 64.0*x*x*x + \n");
  printf(" 136.0*y*y*y + 216.0*z*z*z + 120.0*x*x + 264.0*y*y + 432.0*z*z \n");
  printf("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries \n");
  printf("Analytic solution u(x,y,z) =  x^4 + 2 y^4 + 3 z^4 + \n");
  printf("      4 x^2 y^2 z^2 + 5 x y + 6 x z + 7 y z + 8 \n");
  printf("\n");

  printf("xmin=%g, xmax=%g, ymin=%g, ymax=%g, zmin=%g, zmax=%g \n",
          xmin, xmax, ymin, ymax, zmin, zmax);
  printf("nx=%d, ny=%d, nz=%d \n", nx, ny, nz);
  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 \n");
  hx = (xmax-xmin)/(double)(nx-1);
  for(i=0; i<nx; i++)
  {
    xg[i] = xmin+i*hx;
    Ua[i] = uana(xg[i],ymin,zmin); /* temp */
    printf("i=%d, Ua(%6.3f,%6.3f,%6.3f)=%7.4f \n",
            i, xg[i], ymin, zmin, Ua[i]);
  }  
  printf("\n");
  printf("y grid and analytic solution at xmin, zmin \n");
  hy = (ymax-ymin)/(double)(ny-1);
  for(ii=0; ii<ny; ii++)
  {
    yg[ii] = ymin+ii*hy;
    Ua[ii] = uana(xmin, yg[ii],zmin); /* temp */
    printf("ii=%d, Ua(%6.3f,%6.3f,%6.3f)=%7.4f \n",
            ii, xmin, yg[ii], zmin, Ua[ii]);
  }  
  printf("\n");
  printf("z grid and analytic solution at xmin,ymin \n");
  hz = (zmax-zmin)/(double)(nz-1);
  for(iii=0; iii<nz; iii++)
  {
    zg[iii] = zmin+iii*hz;
    Ua[iii] = uana(xmin, ymin, zg[iii]); /* temp */
    printf("iii=%d, Ua(%6.3f,%6.3f,%6.3f)=%7.4f \n",
            iii, xmin, ymin, zg[iii], Ua[iii]);
  }  
  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;
      for(iii=0; iii<nz; iii++)
      {
        z = zmin+iii*hz;
        Ua[s(i,ii,iii)] = uana(x,y,z);
        /* printf("solution at i=%d,x=%g, ii=%d,y=%g, iii=%d,z=%g is %g \n",
	          i, x, ii, y, iii, z, Ua[s(i,ii,iii)]); */
      }
    }
  }
  printf("\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");

  /* 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++)
            {
              kg[ss(i,ii,iii,j,jj,jjj)] = 0.0;
	    }
	  }
	}
      }
    }
  }
  for(i=0; i<nx; i++)
  {
    for(ii=0; ii<ny; ii++)
    {
      for(iii=0; iii<nz; iii++)
      {
        fg[s(i,ii,iii)] = 0.0;
      }
    }
  }
  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");

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

  for(ii=0; ii<ny; ii++)
  {
    y = yg[ii];
    for(i=0; i<nx; i++)
    {
      x = xg[i];
      iii = 0;
      z = zg[iii];
      kg[ss(i,ii,iii,i,ii,iii)] = 1.0;
      fg[s(i,ii,iii)] = uana(x,y,z);
      printf("boundary i=%d,x=%g, ii=%d,y=%g, iii=%d,z=%g is %g \n",
             i, x, ii, y, iii, z, fg[s(i,ii,iii)]);
    }
    for(i=0; i<nx; i++)
    {
      x = xg[i];
      iii = nz-1;
      z = zg[iii];
      kg[ss(i,ii,iii,i,ii,iii)] = 1.0;
      fg[s(i,ii,iii)] = uana(x,y,z);
      printf("boundary i=%d,x=%g, ii=%d,y=%g, iii=%d,z=%g is %g \n",
             i, x, ii, y, iii, z, fg[s(i,ii,iii)]);
    }
    for(iii=0; iii<nz; iii++)
    {
      z = zg[iii];
      i = 0;
      x = xg[i];
      kg[ss(i,ii,iii,i,ii,iii)] = 1.0;
      fg[s(i,ii,iii)] = uana(x,y,z);
      printf("boundary i=%d,x=%g, ii=%d,y=%g, iii=%d,z=%g is %g \n",
             i, x, ii, y, iii, z, fg[s(i,ii,iii)]);
    }
    for(iii=0; iii<nz; iii++)
    {
      z = zg[iii];
      i = nx-1;
      x = xg[i];
      kg[ss(i,ii,iii,i,ii,iii)] = 1.0;
      fg[s(i,ii,iii)] = uana(x,y,z);
      printf("boundary i=%d,x=%g, ii=%d,y=%g, iii=%d,z=%g is %g \n",
             i, x, ii, y, iii, z, fg[s(i,ii,iii)]);
    }
  } /* end ny */

  for(iii=0; iii<nz; iii++)
  {
    z = zg[iii];
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      i = 0;
      x = xg[i];
      kg[ss(i,ii,iii,i,ii,iii)] = 1.0;
      fg[s(i,ii,iii)] = uana(x,y,z);
      printf("boundary i=%d,x=%g, ii=%d,y=%g, iii=%d,z=%g is %g \n",
             i, x, ii, y, iii, z, fg[s(i,ii,iii)]);
    }
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      i = nx-1;
      x = xg[i];
      kg[ss(i,ii,iii,i,ii,iii)] = 1.0;
      fg[s(i,ii,iii)] = uana(x,y,z);
      printf("boundary i=%d,x=%g, ii=%d,y=%g, iii=%d,z=%g is %g \n",
             i, x, ii, y, iii, z, fg[s(i,ii,iii)]);
    }
    for(i=0; i<nx; i++)
    {
      x = xg[i];
      ii = 0;
      y = yg[ii];
      kg[ss(i,ii,iii,i,ii,iii)] = 1.0;
      fg[s(i,ii,iii)] = uana(x,y,z);
      printf("boundary i=%d,x=%g, ii=%d,y=%g, iii=%d,z=%g is %g \n",
             i, x, ii, y, iii, z, fg[s(i,ii,iii)]);
    }
    for(i=0; i<nx; i++)
    {
      x = xg[i];
      ii = ny-1;
      y = yg[ii];
      kg[ss(i,ii,iii,i,ii,iii)] = 1.0;
      fg[s(i,ii,iii)] = uana(x,y,z);
      printf("boundary i=%d,x=%g, ii=%d,y=%g, iii=%d,z=%g is %g \n",
             i, x, ii, y, iii, z, fg[s(i,ii,iii)]);
    }
  } /* end nz */


  npx = 12;
  npy = 12;
  npz = 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 */

  /* debug */
  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]); 
  i=1; ii=1; iii=1; j=1; jj=1; jjj=1;
  val = galk(xx[2],yy[2],zz[2]);
  printf("galk(xx[2],yy[2],zz[2])=%g, at i=%d,j=%d,ii=%d,jj=%d,iii=%d,jjj=%d \n",
          val, i, j, ii, jj, iii, jjj);
  val = galf(xx[2],yy[2],zz[2]);
  printf("galf=%g, at i=%d,ii=%d,iii=%d \n",
          val, i, ii, iii);

  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++)
            {
              /* compute integral for k(i,j,k)  */
              val = 0.0;
              for(px=1; px<=npx; px++)
              {
                for(py=1; py<=npy; py++)
                {
                  for(pz=1; pz<=npz; pz++)
                  {
                    val = val + wx[px]*wy[py]*wz[pz]*
                          galk(xx[px],yy[py],zz[pz]);
		  } /* pz */
		} /* py */
	      } /* px */
              kg[ss(i,ii,iii,j,jj,jjj)] = val;
	    } /* 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 \n",
              val, i, j, ii, jj, iii, jjj);
    } /* 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");

  /* 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++)
      {
        val = 0.0;
        for(px=1; px<=npx; px++)
        {
          for(py=1; py<=npy; py++)
          {
            for(pz=1; pz<=npz; pz++)
            {
              val = val + wx[px]*wy[py]*wz[pz]*galf(xx[px],yy[py],zz[pz]);
	    } /* pz */
          } /* py */
        } /* px */
        fg[s(i,ii,iii)] = val;
      } /* iii */
    } /* ii */
    printf("Legendre integration=%g, f at i=%d, ii=%d, iii=%d \n",
           val, i, ii, iii);
  } /* 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, nx*ny*nz);

  simeq(nx*ny*nz, 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++)
      {
        err = abs(ug[s(i,ii,iii)]-Ua[s(i,ii,iii)]);
        if(err>maxerr) maxerr = err;
        avgerr = avgerr + err;
        printf("ug[%d,%d,%d]=%8.3f, Ua=%8.3f, err=%g \n",
                i, ii, iii, ug[s(i,ii,iii)], Ua[s(i,ii,iii)], err);
      }
    }
  }
  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));
  printf("\n");
  printf("end fem_check33_la.c \n");
  return 0;
} /* end main */

/* end fem_check33_la.c */