// fem_check44_la.java  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

class fem_check44_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 nt = 5;
  final int nxyzt = nx*ny*nz*nt;
  double xg[] = new double[nx];
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double tg[] = new double[nt];
    int i, j, ii, jj, iii, jjj, iiii, jjjj;
  laphi L = new laphi();

  fem_check44_la()
  {
      double x, y, z, t, hx, hy, hz, ht;
    final int npx = 12; // Gauss Legendre integration order
    final int npy = 12;
    final int npz = 12;
    final int npt = 12;
    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 tt[] = new double[npt+1];
    double wt[] = new double[npt+1]; // for Gauss-Legendre 
    double val, err, avgerr, maxerr;
    int px, py, pz, pt;

    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 tmin = 0.0;
    double tmax = 1.0;

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

    double t_start, t_init, t_kf, t_simeq, t_end;

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

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

    t_start = System.currentTimeMillis();
    System.out.println("x grid and analytic solution at ymin,zmin,tmin ");
    hx = (xmax-xmin)/(double)(nx-1);
    for(i=0; i<nx; i++)
    {
      xg[i] = xmin+i*hx;
      Ua[i] = uana(xg[i], ymin, zmin, tmin); // temp 
      System.out.println("i="+i+", Ua("+xg[i]+")="+Ua[i]);
    } // i  
    System.out.println(" ");
    System.out.println("y grid and analytic solution at xmin,zmin,tmin ");
    hy = (ymax-ymin)/(double)(ny-1);
    for(ii=0; ii<ny; ii++)
    {
      yg[ii] = ymin+ii*hy;
      Ua[ii] = uana(xmin, yg[ii], zmin, tmin); // temp 
      System.out.println("ii="+ii+", Ua("+yg[ii]+")="+Ua[ii]);
    } // ii
    System.out.println(" ");
    System.out.println("z grid and analytic solution at xmin,ymin,tmin ");
    hz = (zmax-zmin)/(double)(nz-1);
    for(iii=0; iii<nz; iii++)
    {
      zg[iii] = zmin+iii*hz;
      Ua[iii] = uana(xmin, ymin, zg[iii], tmin); // temp 
      System.out.println("iii="+iii+", Ua("+zg[iii]+")="+Ua[iii]);
    } // iii
    System.out.println(" ");
    System.out.println("t grid and analytic solution at xmin,ymin,zmin ");
    ht = (tmax-tmin)/(double)(nt-1);
    for(iiii=0; iiii<nt; iiii++)
    {
      tg[iiii] = tmin+iiii*ht;
      Ua[iiii] = uana(xmin, ymin, zmin, tg[iiii]); // temp 
      System.out.println("iiii="+iiii+", Ua("+tg[iiii]+")="+Ua[iiii]);
    } // iiii
    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;
          for(iiii=0; iiii<nt; iiii++)
          {
	    t = tmin+(double)iiii*ht;
            Ua[s(i,ii,iii,iiii)] = uana(x,y,z,t);
	  }
	}
      }
    }
    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++)
              {
                for(iiii=0; iiii<nt; iiii++)
                {
                  for(jjjj=0; jjjj<nt; jjjj++)
                  {
		    kg[ss(i,ii,iii,iiii,j,jj,jjj,jjjj)] = 0.0;
		  }
		}
	      }
	    }
	  }
        }
      }
    }

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

    // set boundary conditions, overkill faces 
    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);
	  } // iiii
        } // iii
      } // ii
    } // i

    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 
    System.out.println("calling gauleg tmin="+tmin+", tmax="+tmax+
                       ", npt="+npt);
    new gaulegf(tmin, tmax, tt, wt, npt); // tt and wt 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]);
    System.out.println("tt[1]="+tt[1]+", tt[2]="+tt[2]+
                       ", wt[1]="+wt[1]+", wt[2]="+wt[2]);
    i=1; ii=1; iii=1; iiii=1; j=1; jj=1; jjj=1; jjjj=1; 
    val = galk(xx[2],yy[2],zz[2],tt[2]);
    System.out.println("galk(xx[2],yy[2],zz[2],tt[2])="+val);
    val = galf(xx[2],yy[2],zz[2],tt[2]);
    System.out.println("galf(xx[2],yy[2],zz[2],tt[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++)
              {
                for(iiii=1; iiii<nt-1; iiii++)
                {
                  for(jjjj=0; jjjj<nt; jjjj++)
                  {
                    // compute integral for k(i,j,k,l)  
                    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 
        System.out.println("Legendre integration="+val+", at i="+i+
       	                   ", j="+j+", ii="+ii+", jj="+jj);
        System.out.println("        , iii"+iii+", jjj"+jjj+
                           ", iiii"+iiii+", jjjj"+jjjj);
      } // 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++)
	{
	  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;
	  } // end iiii loop
        } // end iii loop
      } // end ii loop
      System.out.println("Legendre integration="+val+", f at i="+i+
                         ", ii="+ii+", iii="+iii+", iiii="+iiii);
    } // 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 "+nxyzt+" equations in "+nxyzt+" unknowns");
    new simeq(nxyzt, 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++)
        {
          for(iiii=0; iiii<nt; iiii++)
          {
	      err = Math.abs(ug[s(i,ii,iii,iiii)]-Ua[s(i,ii,iii,iiii)]);
            if(err>maxerr) maxerr = err;
            avgerr = avgerr + err;
            System.out.println("ug["+i+","+ii+","+iii+","+iiii+"]="+
                                ug[s(i,ii,iii,iiii)]+
			 ", Ua="+Ua[s(i,ii,iii,iiii)]+
			 ", err="+(ug[s(i,ii,iii,iiii)]-Ua[s(i,ii,iii,iiii)]));
	  }
	}
      }
    }
    System.out.println("                         maxerr="+maxerr+", avgerr="+
                       avgerr/(double)(nxyzt));
    t_end = System.currentTimeMillis();
    System.out.println("total took "+((t_end-t_start)/1000.0)+" seconds");

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

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

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

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

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

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

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