// fem_check33_la.java  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)
// 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

class fem_check33_la
{
  // i, j, nx, ii, jj, ny, xg, yg needed by galk or galf, available at call 
  final int nx = 5;
  final int ny = 5;
  final int nz = 5;
  final int nxyz = nx*ny*nz;
  double xg[] = new double[nx];
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  int i, j, ii, jj, iii, jjj;
  laphi L = new laphi();

  fem_check33_la()
  {
      double x, y, z, hx, hy, hz;
    final int npx = 3; // Gauss Legendre integration order
    final int npy = 3;
    final int npz = 3;
    double xx[] = new double[npx+1];
    double wx[] = new double[npx+1]; // for Gauss-Legendre 
    double yy[] = new double[npy+1];
    double wy[] = new double[npy+1]; // for Gauss-Legendre 
    double zz[] = new double[npz+1];
    double wz[] = new double[npz+1]; // for Gauss-Legendre 
    double val, err, avgerr, maxerr;
    int px, py, pz;

    double xmin = 0.0; // problem parameters 
    double xmax = 1.0;
    double ymin = 0.0;
    double ymax = 1.0;
    double zmin = 0.0;
    double zmax = 1.0;

    double kg[] = new double[nxyz*nxyz];  
    double fg[] = new double[nxyz];
    double ug[] = new double[nxyz];
    double Ua[] = new double[nxyz];

    double t_start, t_init, t_kf, t_simeq, t_end;

    System.out.println("fem_check33_la.c running");
    System.out.println("Given uxxx(x,y,z) + 2 uyyy(x,y,z) + 3 uzzz(x,y,z) +");
    System.out.println(" 4 uxxy(x,y,z) + 5 uxxz(x,y,z) + 6 uyyx(x,y,z) +");
    System.out.println(" 7 uyyz(x,y,z) + 8 uzzx(x,y,z) + 9 uzzy(x,y,z) +");
    System.out.println(" 10 uxx(x,y,z) + 11 uyy(x,y,z) + 12 uzz(x,y,z) +");
    System.out.println(" 13 uxy(x,y,z) + 14 uxz(x,y,z) + 15 uyz(x,y,z) +");
    System.out.println(" 16  ux(x,y,z) + 17  uy(x,y,z) + 18  uz(x,y,z) +");
    System.out.println(" 19   u(x,y,z) = F(x,y,z)");
    System.out.println("xmin<=x<=xmax ymin<=y<=ymax Boundaries  ");
    System.out.println("Analytic solution u(x,y,z) =  x^4 + 2 y^4 + 3 z^4 +");
    System.out.println("      4 x^2 y^2 z^2 + 5 x y + 6 x z + 7 y z + 8");
    System.out.println(" ");

    System.out.println("xmin="+xmin+", xmax="+xmax+", ymin="+ymin+", ymax="+ymax);
    System.out.println("zmin="+zmin+", zmax="+zmax);
    System.out.println("nx="+nx+", ny="+ny+", nz="+nz);

    t_start = System.currentTimeMillis();
    System.out.println("x grid and analytic solution at ymin,zmin  ");
    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 
      System.out.println("i="+i+", Ua("+xg[i]+")="+Ua[i]);
    }  
    System.out.println(" ");
    System.out.println("y grid and analytic solution at xmin,zmin ");
    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 
      System.out.println("ii="+ii+", Ua("+yg[ii]+")="+Ua[ii]);
    }  
    System.out.println(" ");
    System.out.println("z grid and analytic solution at xmin,ymin ");
    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 
      System.out.println("iii="+iii+", Ua("+zg[iii]+")="+Ua[iii]);
    }  
    System.out.println(" ");

    // 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)] = uana(x,y,z);
          System.out.println("solution at i="+i+",x="+x+", ii="+ii+",y="+y);
          System.out.println(", iii="+iii+",z="+z+", is"+Ua[s(i,ii,iii)]);
	}
      }
    }
    System.out.println(" ");

    // 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;
	}
      }
    }

    // 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);
	System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y);
	System.out.println(", iii="+iii+",z="+z+", is "+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);
	System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y);
	System.out.println(", iii="+iii+",z="+z+", is "+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);
	System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y);
	System.out.println(", iii="+iii+",z="+z+", is "+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);
	System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y);
	System.out.println(", iii="+iii+",z="+z+", is "+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);
	System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y);
	System.out.println(", iii="+iii+",z="+z+", is "+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);
	System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y);
	System.out.println(", iii="+iii+",z="+z+", is "+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);
	System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y);
	System.out.println(", iii="+iii+",z="+z+", is "+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);
	System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y);
	System.out.println(", iii="+iii+",z="+z+", is "+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);
	System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y);
	System.out.println(", iii="+iii+",z="+z+", is "+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);
	System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y);
	System.out.println(", iii="+iii+",z="+z+", is "+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);
	System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y);
	System.out.println(", iii="+iii+",z="+z+", is "+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);
	System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y);
	System.out.println(", iii="+iii+",z="+z+", is "+fg[s(i,ii,iii)]);
      }
    } // end nz

    System.out.println("calling gauleg xmin="+xmin+", xmax="+xmax+
                       ", npx="+npx);
    new gaulegf(xmin, xmax, xx, wx, npx); // xx and wx set up for integrations 
    System.out.println("calling gauleg ymin="+ymin+", ymax="+ymax+
                       ", npy="+npy);
    new gaulegf(ymin, ymax, yy, wy, npy); // yy and wy set up for integrations 
    System.out.println("calling gauleg zmin="+zmin+", zmax="+zmax+
                       ", npz="+npz);
    new gaulegf(zmin, zmax, zz, wz, npz); // zz and wz set up for integrations 

    // debug
    System.out.println("xx[1]="+xx[1]+", xx[2]="+xx[2]+
                       ", wx[1]="+wx[1]+", wx[2]="+wx[2]);
    System.out.println("yy[1]="+yy[1]+", yy[2]="+yy[2]+
                       ", wy[1]="+wy[1]+", wy[2]="+wy[2]);
    System.out.println("zz[1]="+zz[1]+", zz[2]="+zz[2]+
                       ", wz[1]="+wz[1]+", wz[2]="+wz[2]);
    i=1; ii=1; iii=1; j=1; jj=1; jjj=1;
    val = galk(xx[2],yy[2],zz[2]);
    System.out.println("galk(xx[2],yy[2],zz[2])="+val);
    val = galf(xx[2],yy[2],zz[2]);
    System.out.println("galf(xx[2],yy[2],zz[2])="+val);

    t_init = System.currentTimeMillis();
    System.out.println("initialization took "+((t_init-t_start)/1000.0)+" seconds");
    // compute entries in stiffness matrix 
    System.out.println("compute stiffness matrix  ");
    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
                System.out.println("Legendre integration="+val+", at i="+i+
                     ", j="+j+", ii="+ii+", jj="+jj+", iii"+iii+", jjj"+jjj);
                kg[ss(i,ii,iii,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)  
    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
          System.out.println("Legendre integration="+val+", f at i="+i+
                           ", ii="+ii+", iii="+iii);
          fg[s(i,ii,iii)] = val;
        } // end iii loop
      } // end ii loop
    } // end i loop
    System.out.println("k computed stiffness matrix, see above  ");
    System.out.println("f computed forcing function, see above  ");
    t_kf = System.currentTimeMillis();
    System.out.println("kf build took "+((t_kf-t_init)/1000.0)+" seconds");

    System.out.println("solve "+nxyz+" equations in "+nxyz+" unknowns");
    new simeq(nxyz, kg, fg, ug);
    System.out.println("equations solved");
    t_simeq = System.currentTimeMillis();
    System.out.println("simeq took "+((t_simeq-t_kf)/1000.0)+" seconds");

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("ug computed Galerkin, Ua analytic, error  ");
    for(i=0; i<nx; i++)
    {
      for(ii=0; ii<ny; ii++)
      {
        for(iii=0; iii<nz; iii++)
        {
	  err = Math.abs(ug[s(i,ii,iii)]-Ua[s(i,ii,iii)]);
          if(err>maxerr) maxerr = err;
          avgerr = avgerr + err;
          System.out.println("ug["+i+","+ii+","+iii+"]="+ug[s(i,ii,iii)]+
			     ", Ua="+Ua[s(i,ii,iii)]+
			     ", err="+(ug[s(i,ii,iii)]-Ua[s(i,ii,iii)]));
	}
      }
    }
    System.out.println("                         maxerr="+maxerr+", avgerr="+
                       avgerr/(double)(nxyz));
    t_end = System.currentTimeMillis();
    System.out.println("total took "+((t_end-t_start)/1000.0)+" seconds");

    System.out.println(" ");
  } // end fem_check33_la constructor

  // utility functions to compute indices
  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

  // PDE problem definition functions
  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;
  }

  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;
  }

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

  public static void main (String[] args)
  {
    new fem_check33_la();
  }
} // end class fem_check33_la