// 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
 
// needs laphi.java
// needs simeq.java
// needs gaulegf.java

class fem_check44_la
{
  laphi L = new laphi();
  final int nx = 5;
  final int ny = 5;
  final int nz = 5;
  final int nt = 5;
  final int nxyzt = nx*ny*nz*nt; // full matrix for integration
                                 // includes boundary

  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)*nxyzt+
           (j*ny*nz*nt+jj*nz*nt+jjj*nt+jjjj);
  } // end ss 

  double kg[] = new double[nxyzt*nxyzt]; // index by ss(...)
  double xg[] = new double[nx]; // index by i
  double yg[] = new double[ny]; // index by ii
  double zg[] = new double[nz]; // index by iii
  double tg[] = new double[nt]; // index by iiii
  double fg[] = new double[nxyzt]; // RHS index by s(i,ii,iii,iiii)
  double ug[] = new double[nxyzt]; // unknown u(x,y,z,t)
  double Ua[] = new double[nxyzt]; // analytic solution for checking
  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 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,
              int i, int j, int ii, int jj, int iii, int jjj,
              int iiii, int jjjj)
  {                         // 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,
              int i, int ii, int iii, int iiii)
                                    // 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 

  fem_check44_la()
  { 
    double x, y, z, t, hx, hy, hz, ht;
    double xx[] = new double[100]; // for Gauss-Legendre 
    double wx[] = new double[100];
    double yy[] = new double[100];
    double wy[] = new double[100];
    double zz[] = new double[100];
    double wz[] = new double[100];
    double tt[] = new double[100];
    double wt[] = new double[100];
    double val, err, avgerr, maxerr;
    int npx, npy, npz, npt;
    double time_start;
    double time_now;
    double time_last;
    double valU, valUa;

    System.out.println("fem_check44_la.java 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);
    time_start  = System.currentTimeMillis()/1000.0;
    System.out.println("time start = "+time_start+" seconds ");
    time_last = time_start;

    System.out.println("x grid");
    hx = (xmax-xmin)/(double)(nx-1);
    for(int i=0; i<nx; i++)
    {
      xg[i] = xmin+(double)i*hx;
      System.out.println("xg["+i+"]="+xg[i]);
    } // i   
    System.out.println(" ");
    System.out.println("y grid");
    hy = (ymax-ymin)/(double)(ny-1);
    for(int ii=0; ii<ny; ii++)
    {
      yg[ii] = ymin+(double)ii*hy;
      System.out.println("yg["+ii+"="+yg[ii]);
    } // ii 
    System.out.println(" ");
    System.out.println("z grid");
    hz = (zmax-zmin)/(double)(nz-1);
    for(int iii=0; iii<nz; iii++)
    {
      zg[iii] = zmin+(double)iii*hz;
      System.out.println("zg["+iii+"]="+zg[iii]);
    } // iii 
    System.out.println(" ");
    System.out.println("t grid");
    ht = (tmax-tmin)/(double)(nt-1);
    for(int iiii=0; iiii<nt; iiii++)
    {
      tg[iiii] = tmin+(double)iiii*ht;
      System.out.println("tg["+iiii+"]="+tg[iiii]);
    } // iiii 
    System.out.println(" ");

    System.out.println("put solution in Ua vector ");
    // put solution in Ua vector 
    for(int i=0; i<nx; i++)
    {
      x = xg[i];
      for(int ii=0; ii<ny; ii++)
      {
        y = yg[ii];
        for(int iii=0; iii<nz; iii++)
        {
          z = zg[iii];
          for(int iiii=0; iiii<nt; iiii++)
          {
            t = tg[iiii];
            Ua[s(i,ii,iii,iiii)] = uana(x,y,z,t);
            System.out.println("solution at i="+i+", ii="+ii+", iii="+iii+", iiii="+
                               iiii+", "+Ua[s(i,ii,iii,iiii)]); 
          } // iiii 
        } // iii 
      } // ii 
    } // i 
    time_now = System.currentTimeMillis()/1000.0;
    System.out.println("time now = "+time_now+" , time for previous section = "+
                       (time_now-time_last)+" seconds");
    time_last = time_now;
    System.out.println(" ");

    System.out.println("initialize k and f ");
    // initialize k and f 
    for(int i=0; i<nx; i++)
    {
      for(int j=0; j<nx; j++)
      {
        for(int ii=0; ii<ny; ii++)
        {
          for(int jj=0; jj<ny; jj++)
          {
            for(int iii=0; iii<nz; iii++)
            {
              for(int jjj=0; jjj<nz; jjj++)
              {
                for(int iiii=0; iiii<nt; iiii++)
                {
                  for(int jjjj=0; jjjj<nt; jjjj++)
                  {
                    kg[ss(i,ii,iii,iiii,j,jj,jjj,jjjj)] = 0.0;
	          } // jjjj 
	        } // iiii 
	      } // jjj 
            } // iii 
          } // jj 
        } // ii 
      } // j 
    } // i 
    for(int i=0; i<nx; i++)
    {
      for(int ii=0; ii<ny; ii++)
      {
        for(int iii=0; iii<nz; iii++)
        {
          for(int iiii=0; iiii<nt; iiii++)
          {
            fg[s(i,ii,iii,iiii)] = 0.0;
          } // iiii 
        } // iii 
      } // ii 
    } // i 
    time_now = System.currentTimeMillis()/1000.0;
    System.out.println("time now = "+time_now+", time for previous section = "+
		       (time_now-time_last)+" seconds");
    time_last = time_now;
    System.out.println(" ");

    System.out.println("set boundary conditions, all sides ");
    // set boundary conditions, all sides 
    for(int i=0; i<nx; i++)
    {
      x = xg[i];
      for(int ii=0; ii<ny; ii++)
      {
        y = yg[ii];
        for(int iii=0; iii<nz; iii++)
        {
          z = zg[iii];
          for(int iiii=0; iiii<nt; iiii++)
          {
            t = tg[iiii];
            kg[ss(i,ii,iii,iiii,i,ii,iii,iiii)] = 1.0;
            val = uana(x,y,z,t);
            fg[s(i,ii,iii,iiii)] = val;
          } // iiii 
        } // iii 
      } // ii 
    } // i  // end boundry overkill 


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

    val = galk(xx[2],yy[2],zz[2],tt[2],1,1,1,1,1,1,1,1);
    System.out.println("galk(xx[2],yy[2],zz[2],tt[2])="+val);
    val = galf(xx[2],yy[2],zz[2],tt[2],1,1,1,1);
    System.out.println("galf(xx[2],yy[2],zz[2],tt[2])="+val);

    time_now = System.currentTimeMillis()/1000.0;
    System.out.println("time now = "+time_now+" time for previous section = "+
                       (time_now-time_last)+" seconds");
    time_last = time_now;
    System.out.println(" ");


    // compute entries in stiffness matrix 
    System.out.println("compute stiffness matrix ");
    for(int i=1; i<nx-1; i++)
    {
      for(int j=0; j<nx; j++)
      {
        for(int ii=1; ii<ny-1; ii++)
        {
          for(int jj=0; jj<ny; jj++)
          {
            for(int iii=1; iii<nz-1; iii++)
            {
              for(int jjj=0; jjj<nz; jjj++)
              {
                for(int iiii=1; iiii<nt-1; iiii++)
                {
                  for(int jjjj=0; jjjj<nt; jjjj++)
                  {
                    // compute integral for k(x,y,z,t)  
                    val = 0.0;
                    for(int px=1; px<=npx; px++)
                    {
                      for(int py=1; py<=npy; py++)
                      {
                        for(int pz=1; pz<=npz; pz++)
                        {
                          for(int pt=1; pt<=npt; pt++)
                          {
                            val = val + wx[px]*wy[py]*wz[pz]*wt[pt]*
			          galk(xx[px],yy[py],zz[pz],tt[pt],
				  i,j,ii,jj,iii,jjj,iiii,jjjj);
                          } // 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=%g, at i=%d,j=%d,ii=%d,jj=%d,iii=%d,jjj=%d,
      } // end j loop 
    } // end i loop 
    time_now = System.currentTimeMillis()/1000.0;
    System.out.println("time now = "+time_now+" time for previous section = "+
  		       (time_now-time_last)+" seconds");
    time_last = time_now;
    System.out.println(" ");

    System.out.println("compute integral for f(i) ");
    // compute integral for f(i)  
    for(int i=1; i<nx-1; i++)
    {
      for(int ii=1; ii<ny-1; ii++)
      {
        for(int iii=1; iii<nz-1; iii++)
        {
          for(int iiii=1; iiii<nt-1; iiii++)
          {
            val = 0.0;
            for(int px=1; px<=npx; px++)
            {
              for(int py=1; py<=npy; py++)
              {
                for(int pz=1; pz<=npz; pz++)
                {
                  for(int pt=1; pt<=npt; pt++)
                  {
                    val = val + wx[px]*wy[py]*wz[pz]*wt[pt]*
		          galf(xx[px],yy[py],zz[pz],tt[pt],i,ii,iii,iiii);
  	          } // pt 
	        } // pz 
              } // py 
            } // px 
            fg[s(i,ii,iii,iiii)] = val;
          } // iiii 
        } // iii 
      } // ii 
      // System.out.println("Legendre integration=%g, f at i=%d, ii=%d, iii=%d,
    } // i 
    System.out.println("k computed stiffness matrix, see above ");
    System.out.println("f computed forcing function, see above ");
    time_now = System.currentTimeMillis()/1000.0;
    System.out.println("time now = "+time_now+" time for previous section = "+
		       (time_now-time_last)+" seconds");
    time_last = time_now;
    System.out.println(" ");

    System.out.println("solve "+nxyzt+" equations in "+nxyzt+" unknowns ");

    new simeq(nxyzt, kg, fg, ug);

    System.out.println("equations solved ");
    time_now = System.currentTimeMillis()/1000.0;
    System.out.println("time now = "+time_now+" time for previous section = "+
		       (time_now-time_last)+" seconds");
    time_last = time_now;
    System.out.println(" ");


    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("ug computed Galerkin, Ua analytic, error ");
    for(int i=0; i<nx; i++)
    {
      for(int ii=0; ii<ny; ii++)
      {
        for(int iii=0; iii<nz; iii++)
        {
          for(int iiii=0; iiii<nt; iiii++)
          {
	    valU = ug[s(i,ii,iii,iiii)];
            valUa =Ua [s(i,ii,iii,iiii)];
	    err = Math.abs(valU - valUa);
            if(err>maxerr) maxerr = err;
            avgerr = avgerr + err;
            System.out.println("ug["+i+","+ii+","+iii+","+iiii+"]="+valU+", Ua="+
                               valUa+", err="+err);
          } // iiii 
        } // iii 
      } // ii 
    } // i 
    time_now = System.currentTimeMillis()/1000.0;
    System.out.println("total CPU time  = "+(time_now-time_start)+" seconds");
    System.out.println(" ");

    System.out.println("                              maxerr="+maxerr+
                       ", avgerr="+(avgerr/(double)(nx*ny*nz*nt)));
    System.out.println(" ");
    System.out.println("end fem_check44_la.java");
  } // fem_check44_la

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