/* fem_bihar3d_la.c  Galerkin FEM biharmonic 4th order
 *                                                                   
 * solve  uxxxx(x,y,z) + uyyyy(x,y,z) + uzzzz(x,y,z) + 2*uxxyy(x,y,z) + 
 *        2*uxxzz(x,y,z) + 2*uyyzz(x,y,z) = f(x,y,z)
 *
 *        f(x,y,z) := 0.24  ,
 *
 *        in the cube   xmin<x<xmax and ymin<y<ymax and zmin<z<xmax  ,
 *
 *  Dirichlet boundary conditions are computed using:
 *  ub(x,y,z) := 0.01*x*x*x*x + 0.02*y*y*y*y + 0.03*z*z*z*z - 0.04*x*x*y*y -
 *               0.05*y*y*z*z - 0.06*x*x*z*z + 0.1*x*x*x*y + 
 *               0.3*y*y*y*z + 0.2*z*z*z*x - x*y - y*z - x*z + 1
 *
 *  uses laphi.c  simeq.c  deriv.c  gaulegf.c 
 */

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include "simeq.h"
#include "laphi.h"

#undef  abs
#define abs(a) ((a)<0.0?(-(a)):(a))
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))

#define nx  6
#define ny  6
#define nz  6
#define npx 6
#define npy 6
#define npz 6

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;

/* i, j, nx, ii, jj, ny, iii, jjj, xg, yg,zg */
/* needed by galk or galf, available at call */
static int i, j, ii, jj, iii, jjj;
static double xg[nx]; /* X grig */
static double yg[ny]; /* Y grid */
static double zg[ny]; /* z grid */
static double hx, hy, hz;
static double ug[nx*ny*nz]; /* solution */
static double maxerr;


/* function prototypes */
int s(int i, int j, int k); // subscript for kg and Ua
double F(double x, double y, double z); /* RHS */
double ub(double x, double y, double z);
double galk(double x, double y, double z); /* Galerkin k function for this problem */
double galf(double x, double y, double z); /* Galerkin f function for this problem */
void write_soln();
void check_soln();

int main(int argc, char * argv[])
{
    int nxyz = nx*ny*nz;
    double x, y, z;
    /* npx, npy, npz  Gauss Legendre integration order */
    double xx[npx];
    double wx[npx];
    double yy[npy];
    double wy[npy];
    double zz[npy];
    double wz[npy];
    double val, err, avgerr;
    int px, py, pz;
    double kg[nx*ny*nz][nx*ny*nz];
    double fg[nx*ny*nz];
    double Ua[nx*ny*nz];
    int stat;

    printf("fem_bihar3d_la.c running.\n");
    printf("solve  uxxxx(x,y,z) + uyyyy(x,y,z) + uzzzz(x,y,z) + 2*uxxyy(x,y,z) + \n");
    printf("       2*uxxzz(x,y,z) + 2*uyyzz(x,y,z) = f(x,y,z)\n");
    printf(" \n");
    printf("       f(x,y,z) := 0.24  ,\n");
    printf(" \n");
    printf("      in the cube   xmin<x<xmax and ymin<y<ymax and zmin<z<xmax  ,\n");
    printf(" \n");
    printf("Dirichlet boundary conditions are computed using:\n");
    printf("ub(x,y,z) := 0.01*x*x*x*x + 0.02*y*y*y*y + 0.03*z*z*z*z - 0.04*x*x*y*y -\n");
    printf("             0.05*y*y*z*z - 0.06*x*x*z*z + 0.1*x*x*x*y + \n");
    printf("             0.3*y*y*y*z + 0.2*z*z*z*x - x*y - y*z - x*z + 1\n");
    printf(" \n");
    printf("xmin=%6.3f, xmax=%6.3f, ymin=%6.3f, ymax=%6.3f, zmin=%6.3f, zmax=%6.3f \n",
           xmin, xmax, ymin, ymax, zmin, zmax);
    printf("nx=%d, ny=%d, nz=%d \n",nx, ny, nz);
    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;
      printf("i=%d, xg[i]=%f, ub(xg[i],ymin,zmin)=%f\n",
             i, xg[i], ub(xg[i],ymin,zmin));
    }  
    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;
      printf("ii=%d, yg[ii]=%f, ub(xmin,yg[ii],zmin)=%f \n",
             ii, yg[ii], ub(xmin, yg[ii],zmin));
    }  
    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;
      printf("iii=%d, zg[iii]=%f, ub(xmin,ymin,zg[iii])=%f \n",
             iii, zg[iii], ub(xmin, ymin, zg[iii]));
    }  
    printf(" \n");
    /* put solution in Ua vector */ 
    for(i=0; i<nx; i++)
    {
      x = xmin+(double)i*hx;
      for(ii=0; ii<ny; ii++)
      {
	y = ymin+(double)ii*hy;
        for(iii=0; iii<nz; iii++)
        {
	  z = zmin+(double)iii*hz;
          Ua[s(i,ii,iii)] = ub(x,y,z);
	}
      }
    }
    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[s(i,ii,iii)][s(j,jj,jjj)] = 0.0;
	      }
	    }
	  }
        }
      }
    }
    for(i=0; i<nx; i++)
    {
      for(ii=0; ii<ny; ii++)
      {
        for(iii=0; iii<nz; iii++)
        {
          fg[i*ny+ii] = 0.0;
	}
      }
    }
    /* set boundary conditions, six sides */ 
    for(i=0; i<nx; i++) /* x,y at zmin */
    {
      x = xg[i];
      for(ii=0; ii<ny; ii++)
      {
        y = yg[ii];
        iii = 0;
        z = zg[iii];
        kg[s(i,ii,iii)][s(i,ii,iii)] = 1.0;
        fg[s(i,ii,iii)] = ub(x,y,z);
      }
    } /* end x,y at zmin */
    for(i=0; i<nx; i++) /* x,y at zmax */
    {
      x = xg[i];
      for(ii=0; ii<ny; ii++)
      {
        y = yg[ii];
        iii = nz-1;
        z = zg[iii];
        kg[s(i,ii,iii)][s(i,ii,iii)] = 1.0;
        fg[s(i,ii,iii)] = ub(x,y,z);
      }
    } /* end x,y at zmax */
    for(i=0; i<nx; i++) /* x,z at ymin */
    {
      x = xg[i];
      ii = 0;
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        kg[s(i,ii,iii)][s(i,ii,iii)] = 1.0;
        fg[s(i,ii,iii)] = ub(x,y,z);
      }
    } /* end x,z at ymin */
    for(i=0; i<nx; i++) /* x,z at ymax */
    {
      x = xg[i];
      ii = ny-1;
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        kg[s(i,ii,iii)][s(i,ii,iii)] = 1.0;
        fg[s(i,ii,iii)] = ub(x,y,z);
      }
    } /* end x,z at ymax */
    i = 0;
    x = xg[i];
    for(ii=0; ii<ny; ii++) /* y,z at xmin */
    {
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        kg[s(i,ii,iii)][s(i,ii,iii)] = 1.0;
        fg[s(i,ii,iii)] = ub(x,y,z);
      }
    } /* end y,z at xmin */
    i = nx-1;
    x = xg[i];
    for(ii=0; ii<ny; ii++) /* y,z at xmax */
    {
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        kg[s(i,ii,iii)][s(i,ii,iii)] = 1.0;
        fg[s(i,ii,iii)] = ub(x,y,z);
      }
    } /* end y,z at xmax */


    printf("calling gauleg xmin=%f, xmax=%f, npx=%d \n",
           xmin, xmax, npx);
    gaulegf(xmin, xmax, xx, wx, npx); /* xx and wx set up for integrations */ 
    printf("xx[1]=%f, xx[2]=%f, wx[1]=%f, wx[2]=%f \n",
           xx[1],xx[2],wx[1],wx[2]);
    printf("calling gauleg ymin=%f, ymax=%f, npy=%d \n",
           ymin,ymax,npy);
    gaulegf(ymin, ymax, yy, wy, npy); /* yy and wy set up for integrations */
    printf("yy[1]=%f, yy[2]=%f, wy[1]=%f, wy[2]=%f \n",
           yy[1],yy[2],wy[1],wy[2]);
    printf("calling gauleg zmin=%f, zmax=%f, npz=%d \n",
           zmin,zmax,npz);
    gaulegf(zmin, zmax, zz, wz, npz); /* zz and wz set up for integrations */
    printf("zz[1]=%f, zz[2]=%f, wz[1]=%f, wz[2]=%f \n",
           zz[1],zz[2],wz[1],wz[2]);
    i=1; ii=1;  j=1; jj=1; iii=1; jjj=1;
    val = galk(xx[1],yy[1],zz[1]);
    printf("galk(xx[1],yy[1],zz[1])=%g\n",val);
    val = galf(xx[2],yy[2],zz[2]);
    printf("galf(xx[2],yy[2],zz[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++)
          {
            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]);
		    }
	          }
                }
                kg[s(i,ii,iii)][s(j,jj,jjj)] = val;
              } /* end jjj loop */ 
            } /* end iii loop */
          } /* end jj loop */ 
        } /* end ii loop */
      } /* end j loop */
    } /* end i loop */

    /* compute integral for f(i,j,k) */  
    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]);
	      }
            }
          }
          fg[s(i,ii,iii)] = val;
	}
      }
    }
    printf("k computed stiffness matrix, see above \n ");
    printf("f computed forcing function, see above \n ");
    simeq(nx*ny*nz, kg, fg, ug);

    write_soln(); /* to file x y soln */

    check_soln(); /* independent of nowing analytic solution */

    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<ny; 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]=%f, Ua=%f, err=%g \n",
		 i,ii,iii,ug[s(i,ii,iii)],Ua[s(i,ii,iii)],
                 ug[s(i,ii,iii)]-Ua[s(i,ii,iii)]);
	}
      }
    }
    printf("xmax=%f, ymax=%f, zmax=%f, nx=%d, ny=%d, nz=%d, npx=%d, npy=%d, npz=%d \n",
           xmax,ymax,zmax,nx,ny,nz,npx,npy,npz);
    printf("                  maxerr=%g, avgerr=%g \n",
           maxerr,avgerr/(double)(nx*ny*nz));
    printf("\n");

  printf("calling gnuplot \n");
  fflush(stdout);

  stat = execl("/usr/bin/gnuplot","-persist","fem_bihar3d_la_c.plot", NULL);
  /* only returns upon error, no return on sucess */
  printf("returned from gnuplot with stat= %d \n", stat);
  return 0;
} /* end fem_bihar3d_la main */

int s(int i, int j, int k) // subscript for kg and Ua
{
  return i*ny*nz+ii*nz+iii;
}

/* PDE problem definition functions  f and ub */

double F(double x, double y, double z)
{
  return 0.24;
}

double ub(double x, double y, double z)
{
  return 0.01*x*x*x*x + 0.02*y*y*y*y + 0.03*z*z*z*z - 0.04*x*x*y*y -
         0.05*y*y*z*z - 0.06*x*x*z*z + 0.1*x*x*x*y + 0.2*z*z*z*x -
         0.03*y*y*y*z - x*y - y*z - x*z + 1.0;
}

double galk(double x, double y, double z)
{       /* Galerkin k stiffness function for this problem      */ 
        /* solve  uxxxx(x,y,z) + uyyyy(x,y,z) + uzzzz(x,y,z) +     */
        /*      2*uxxyy(x,y,z) + 2*uxxzz(x,y,z)  + 2*uyyzz(x,y,z) = f(x,y,z) */
  return (phipppp(x,j,nx-1,xg)*     phi(y,jj,ny-1,yg)*     phi(z,jjj,nz-1,zg) +
	      phi(x,j,nx-1,xg)* phipppp(y,jj,ny-1,yg)*     phi(z,jjj,nz-1,zg) +
              phi(x,j,nx-1,xg)*     phi(y,jj,ny-1,yg)* phipppp(z,jjj,nz-1,zg) +
       2.0* phipp(x,j,nx-1,xg)*   phipp(y,jj,ny-1,yg)*     phi(z,jjj,nz-1,zg) +
       2.0* phipp(x,j,nx-1,xg)*     phi(y,jj,ny-1,yg)*   phipp(z,jjj,nz-1,zg) +
         2.0* phi(x,j,nx-1,xg)*   phipp(y,jj,ny-1,yg)*   phipp(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 */ 

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 */ 

void write_soln()
{
  FILE * fout;

  fout = fopen("fem_bihar3d_la_c.dat", "w");
  if(fout==NULL)
  {
    printf("file fem_bihar3d_la_c.dat can not be opened for writing\n");
    return;
  }
  printf("writing fem_bihar3d_la_c.dat \n");
  for(i=0; i<nx; i++)
  {
    for(ii=0; ii<ny; ii++)
    {
      iii = nz-1;
      fprintf(fout," %f %f %f \n",xg[i],yg[ii],ug[s(i,ii,iii)]);
    } /* end y */
    fprintf(fout," \n");
  } /* end x */
  fflush(fout);
  fclose(fout);
  printf("finished writing fem_bihar3d_la_c.dat \n");
} /* end write_soln */


void check_soln()
{
               /* for all interior points */
  double x, y, z, err, maxerr, avgerr, rmserr;
  double sumerr, sumsqerr, eval_deriv;
  double cxx[nx][nx];
  double cxxxx[nx][nx];
  double cyy[ny][ny];
  double cyyyy[ny][ny];
  double czz[nz][nz];
  double czzzz[nz][nz];
  double tcxxyy[ny];
  double tcxxzz[nz];
  double tcyyxx[nx];
  int cnt = 0;
  int k, m;

  for(i=1; i<nx-1; i++)
  {
    rderiv(2,nx,i,hx,cxx[i]);
    rderiv(4,nx,i,hx,cxxxx[i]);
  }
  for(ii=1; ii<ny-1; ii++)
  {
    rderiv(2,ny,ii,hy,cyy[ii]);
    rderiv(4,ny,ii,hy,cyyyy[ii]);
  }
  for(iii=1; iii<nz-1; iii++)
  {
    rderiv(2,nz,iii,hz,czz[iii]);
    rderiv(4,nz,iii,hz,czzzz[iii]);
  }

  sumerr = 0.0;
  sumsqerr = 0.0;
  maxerr = 0.0;
  for(i=1; i<nx-1; i++) /* full equation */
  {
    x = xmin + (double)i*hx;
    for(j=1; j<ny-1; j++)
    {
      y = ymin + (double)j*hy;
      for(k=1; k<nz-1; k++)
      {
        z = zmin + (double)k*hz;
        eval_deriv = 0.0; /* no PDE u term */
        
        for(m=0; m<nx; m++)
        {
	  eval_deriv += cxxxx[i][m]*ug[s(m,j,k)];   /* PDE  uxxxx term */
        } 
        for(m=0; m<ny; m++)
        {
	  eval_deriv += cyyyy[j][m]*ug[s(i,m,k)];    /* PDE uyyyy term */
        }
        for(m=0; m<ny; m++)
        {
	  eval_deriv += czzzz[k][m]*ug[s(i,j,m)];    /* PDE uzzzz term */
        }
        /*  + 2*uxxyy */
        for(jj=0; jj<ny; jj++) /* walk jj in y, computing uxx */
        {
	  tcxxyy[jj] = 0.0; /* numeric  uxx part of uxxyy */
          for(m=0; m<nx; m++)
          {
	    tcxxyy[jj] += cxx[i][m]*ug[s(m,jj,k)];
	  }
        }
        for(m=0; m<ny; m++)  /* uyy part of uxxyy */
        {
          eval_deriv += 2.0*cyy[j][m]*tcxxyy[m];  /* PDE uxxyy term */
        }

        err = eval_deriv-F(x,y,z);
        err = abs(err);
        sumerr += err;
        sumsqerr += err*err;
        maxerr = max(err,maxerr);
        cnt++;
      } /* end z */  
    } /* end y */
  } /* end x */
  rmserr = sqrt(sumsqerr/(double)cnt);
  avgerr = sumerr/(double)cnt;
  printf("check_soln values \n");
  printf("maxerr=%g, rmserr=%g, avgerr=%g \n",maxerr,rmserr,avgerr);
} /* end check_soln */ 

/* end fem_bihar3d_la.c */