// pde_bihar44t_eq.java  discretization solution biharmonic 4th order 4D 
//
// solve  uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + 
//        uzzzz(x,y,z,t) + utttt(x,y,z,t) + 
//        2*uxx(x,y,z,t) + 2*uyy(x,y,z,t) + 2*uzz(x,y,z,t) + 
//        2*utt(x,y,z,t) +  2*u(x,y,z,t) = f(x,y,z,t)
//
// f(x,y,z,t) := 72*t+24*x^2*t+48*(y-z)^2*t+2*(x^4+(y-z)^4)*t
//
//  in the hyper cube  xmin< x <xmax and  ymin< y <ymax  ,
//                     zmin< z <zmax and  tmin< t <tmax 
//
//     subject to boundary conditions of the first
//     kind, Dirichlet
//
//  Dirichlet boundary conditions are computed using:
//  ub(x,y,z,t) := (x^4+(y-z)^4)*t
//
//  uses   simeq.c  nderiv.c 
 
import java.io.*;
//                uses rderiv.java and simeq.java

class pde_bihar44t_eq
{
  //  nx, ny, ny, nz, xg, yg, zg, tg needed by many, available at call 
  int nx, ny, nz, nt;
  int i4;
  double xg[];   // X grig
  double yg[];   // Y grid
  double zg[];   // Z grig
  double tg[];   // T grid
  double ug[];   // solution
  double kg[][]; // matrix
  double fg[];   // RHS
  double maxerr;
  double cxxxx[][];
  double cyyyy[][];
  double czzzz[][];
  double ctttt[][];
  double cxx[][];
  double cyy[][];
  double czz[][];
  double ctt[][];
  // cxxyy etc. built dynamically
  boolean debug = true;
  double hx, hy, hz, ht;
  int nxyzt;

  pde_bihar44t_eq(){}; // construct the class once, run many times

  void run(double xmin, double xmax, double ymin, double ymax,
           double zmin, double zmax, double tmin, double tmax,
           int anx, int any, int anz, int ant)
  {
    nx = anx; // needed global
    ny = any;
    nz = anz;
    nt = ant;
    nxyzt = (nx-2)*(ny-2)*(nz-2)*(nt-2);
    xg = new double[nx];
    yg = new double[ny];
    zg = new double[nz];
    tg = new double[nt];
    double x, y, z, t;
    double val, err, avgerr;

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

    cxxxx = new double[nx][nx];
    cyyyy = new double[ny][ny];
    czzzz = new double[nz][nz];
    ctttt = new double[nt][nt];
    cxx = new double[nx][nx];
    cyy = new double[ny][ny];
    czz = new double[nz][nz];
    ctt = new double[nt][nt];

    System.out.println("pde_bihar44t_eq.java running.");
    System.out.println("solve  uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) + utttt(x,y,z,t) + ");
    System.out.println("       2*uxx(x,y,z,t) + 2*uyy(x,y,z,t)  + 2*uzz(x,y,z,t)  + ");
    System.out.println("       2*utt(x,y,z,t) + 2*u(x,y,z,t) = f(x,y,z,t)");
    System.out.println(" ");
    System.out.println("f(x,y,z,t) := 72*t+24*x^2*t+48*(y-z)^2*t+2*(x^4+(y-z)^4)*t");
    System.out.println("in the hyper cube   xmin<x<xmax  and  ymin<y<ymax,");
    System.out.println("                    zmin<z<zmax  and  tmin<t<tmax,");
    System.out.println(" ");
    System.out.println("   subject to boundary conditions of the first");
    System.out.println("   kind, Dirichlet.");
    System.out.println("   ub(x,y,z,t) := (x^4+(y-z)^4)*t");
    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);

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

    // put solution in Ua vector 
    for(int i=1; i<nx-1; i++)
    {
      x = xg[i];
      for(int ii=1; ii<ny-1; ii++)
      {
	y = yg[ii];
        for(int iii=1; iii<nz-1; iii++)
        {
          z = zg[iii];
          for(int iiii=1; iiii<nt-1; iiii++)
          {
	    t = tg[iiii];
            Ua[s(i,ii,iii,iiii)] = ub(x,y,z,t);
	  }
	}
      }
    }
    if(debug) System.out.println(" ");

    // initialize k and f 
    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<nx-1; j++)
      {
        for(int ii=1; ii<ny-1; ii++)
        {
          for(int jj=1; jj<ny-1; jj++)
          {
            for(int iii=1; iii<nz-1; iii++)
            {
              for(int jjj=1; jjj<nz-1; jjj++)
              {
                for(int iiii=1; iiii<nt-1; iiii++)
                {
                  for(int jjjj=1; jjjj<nt-1; jjjj++)
                  {
	            kg[s(i,ii,iii,iiii)][s(j,jj,jjj,jjjj)] = 0.0;
		  }
		}
	      }
	    }
	  }
        }
      }
    }
    for(int i=1; i<nx-1; i++)
    {
      x = xg[i];
      for(int ii=1; ii<ny-1; ii++)
      {
        y = yg[ii];
        for(int iii=1; iii<nz-1; iii++)
        {
          z = zg[iii];
          for(int iiii=1; iiii<nt-1; iiii++)
          {
            t = tg[iiii];
	    fg[s(i,ii,iii,iiii)] = F(x,y,z,t);
	  }
	}
      }
    }

    dsetup(); // setup discrete coefficients 

    // compute entries in kg (stiffness matrix) discretization matrix 
    System.out.println("compute discretization matrix  ");
    for(int i=1; i<nx-1; i++) // not boundary rows
    {
      for(int ii=1; ii<ny-1; ii++) // not boundary rows
      {
        for(int iii=1; iii<nz-1; iii++) // not boundary rows
        {
          for(int iiii=1; iiii<nt-1; iiii++) // not boundary rows
          {
            // compute discrete derivatives around (i,ii) DOF  
	    discretize(i, ii, iii, iiii);
                      // sets kg[s(i,ii,iii,iiii)][s(j,jj,jjj,jjjj)...]
	  } // end iiii loop
	} // end iii loop
      } // end ii loop 
    } // end i loop 

    System.out.println("k computed discrete matrix");
    System.out.println("f computed RHS");

    // write_eqn(); // all
    write_eqnij(); // i,j,val non zero

    new simeq(kg, fg, ug);

    write_soln(); // to file x y soln

    check_soln(); // independent of nowing analytic solution

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("ug computed, Ua analytic, error  ");
    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++)
           {
	     i4 = s(i,ii,iii,iiii);
             err = Math.abs(ug[i4]-Ua[i4]);
             if(err>maxerr) maxerr = err;
             avgerr = avgerr + err;
             System.out.println("ug["+i+","+ii+","+iii+","+iiii+"]="+
                               ug[i4]+", Ua="+Ua[i4]+", err="+(ug[i4]-Ua[i4]));
	   }
	 }
      }
    }
    System.out.println("xmin="+xmin+", xmax="+xmax+", hx="+hx+", nx="+nx);
    System.out.println("ymin="+ymin+", ymax="+ymax+", hy="+hy+", ny="+ny);
    System.out.println("zmin="+zmin+", zmax="+zmax+", hz="+hz+", nz="+nz);
    System.out.println("tmin="+tmin+", tmax="+tmax+", ht="+ht+", nt="+nt);
    System.out.println("                  maxerr="+maxerr+", avgerr="+
                       avgerr/(double)(nxyzt));
    System.out.println("");
  } // end pde_bihar44t_eq constructor

  int s(int i, int ii, int iii, int iiii) // index 1..nx-1  etc
  {
      return (i-1)*(ny-2)*(nz-2)*(nt-2) + (ii-1)*(nz-2)*(nt-2) +
             (iii-1)*(nt-2) + (iiii-1);
  } // end s

  // PDE problem definition functions  f and ub

  double F(double x, double y, double z, double t)
  {
      return 72.0*t+24.0*x*x*t+48.0*(y-z)*(y-z)*t+2.0*(x*x*x*x+
	     (y-z)*(y-z)*(y-z)*(y-z))*t;
  }

  double ub(double x, double y, double z, double t)
  {
      return (x*x*x*x+(y-z)*(y-z)*(y-z)*(y-z))*t;
  }

  void discretize(int i, int ii, int iii, int iiii) 
                     // sets kg[s(i,ii,iii,iiii)][s(j,jj,jjj, jjjj)...]
                     // this code unique to this specific PDE 
  {
    int i4;
    double val;

    i4 = s(i,ii,iii,iiii); // equation index 
    // for each term in each equation 

    // for uxxxx  d^4u/dx^4 
    for(int j=0; j<nx; j++)
    {
      val = cxxxx[i][j];
      if(j==0 || j==nx-1)
        fg[i4] -= val*ub(xg[j],yg[ii],zg[iii],tg[iiii]);
      else kg[i4][s(j,ii,iii,iiii)] += val;
    }

    // for uyyyy  d^4u/dy^4 
    for(int jj=0; jj<ny; jj++)
    {
      val = cyyyy[ii][jj];
      if(jj==0 || jj==ny-1)
        fg[i4] -=  val*ub(xg[i],yg[jj],zg[iii],tg[iiii]);
      else kg[i4][s(i,jj,iii,iiii)] += val;
    }

    // for uzzzz  d^4u/dz^4 
    for(int jjj=0; jjj<nz; jjj++)
    {
      val = czzzz[iii][jjj];
      if(jjj==0 || jjj==nz-1)
        fg[i4] -=  val*ub(xg[i],yg[ii],zg[jjj],tg[iiii]);
      else kg[i4][s(i,ii,jjj,iiii)] += val;
    }

    // for utttt  d^4u/dt^4 
    for(int jjjj=0; jjjj<nt; jjjj++)
    {
      val = ctttt[iiii][jjjj];
      if(jjjj==0 || jjjj==nt-1)
        fg[i4] -=  val*ub(xg[i],yg[ii],zg[iii],tg[jjjj]);
      else kg[i4][s(i,ii,iii,jjjj)] += val;
    }

    // for uxx  d^2u/dx^2 
    for(int j=0; j<nx; j++)
    {
      val = 2.0*cxx[i][j];
      if(j==0 || j==nx-1)
 	fg[i4] -= val*ub(xg[j],yg[ii],zg[iii],tg[iiii]);
      else kg[i4][s(j,ii,iii,iiii)] += val;
    } // end j 

    // for uyy  d^2u/dy^2 
    for(int jj=0; jj<ny; jj++)
    {
      val = 2.0*cyy[ii][jj];
      if(jj==0 || jj==ny-1)
	fg[i4] -= val*ub(xg[i],yg[jj],zg[iii],tg[iiii]);
      else kg[i4][s(i,jj,iii,iiii)] += val;
    } // end jj 

    // for uzz  d^2u/dz^2 
    for(int jjj=0; jjj<nz; jjj++)
    {
      val = 2.0*czz[iii][jjj];
      if(jjj==0 || jjj==nz-1)
	fg[i4] -= val*ub(xg[i],yg[ii],zg[jjj],tg[iiii]);
      else kg[i4][s(i,ii,jjj,iiii)] += val;
    } // end jjj 

    // for utt  d^2u/dt^2 
    for(int jjjj=0; jjjj<nt; jjjj++)
    {
      val = 2.0*ctt[iiii][jjjj];
      if(jjjj==0 || jjjj==nt-1)
	fg[i4] -= val*ub(xg[i],yg[ii],zg[iii],tg[jjjj]);
      else kg[i4][s(i,ii,iii,jjjj)] += val;
    } // end jjjj 

    kg[i4][i4] += 2.0; // U  
    if(debug) System.out.println("discretize equation i4="+i4+", RHS="+fg[i4]);
  } // end discretize 

  void dsetup()
  {
    // each such as cxxxx[i][j] is j_th coefficient at DOF i 
    for(int i=0; i<nx; i++) rderiv(4, nx, i, hx, cxxxx[i]);
    if(debug)
    {
      System.out.println("debug print of cxxxx[i][j] ");
      for(int i=0; i<nx; i++)
      {
        System.out.println("i="+i);
        for(int j=0; j<nx; j++) System.out.println(i+", "+j+", "+cxxxx[i][j]);
        System.out.println(" ");
      }
    }
    for(int i=0; i<nx; i++) rderiv(2, nx, i, hx, cxx[i]);
    if(debug)
    {
      System.out.println("debug print of cxx[i][j] ");
      for(int i=0; i<nx; i++)
      {
        System.out.println("i="+i);
        for(int j=0; j<nx; j++) System.out.println(i+", "+j+", "+cxx[i][j]);
        System.out.println(" ");
      }
    }
    for(int ii=0; ii<ny; ii++) rderiv(4, ny, ii, hy, cyyyy[ii]);
    for(int ii=0; ii<ny; ii++) rderiv(2, ny, ii, hy, cyy[ii]);
    for(int iii=0; iii<nz; iii++) rderiv(4, nz, iii, hz, czzzz[iii]);
    for(int iii=0; iii<nz; iii++) rderiv(2, nz, iii, hz, czz[iii]);
    for(int iiii=0; iiii<nt; iiii++) rderiv(4, nt, iiii, ht, ctttt[iiii]);
    for(int iiii=0; iiii<nt; iiii++) rderiv(2, nt, iiii, ht, ctt[iiii]);
  } // end dsetup

  void write_soln()
  {
    double x, y, z, t;
    int dim = 4;
    try
    {
      FileWriter fout = new FileWriter("pde_bihar44t_eq.dat");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde_bihar44t_eq.dat");
      fileout.write(dim+" "+nx+" "+ny+" "+nz+" "+nt+" \n");

      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];
	      if(i==0 || i==nx-1 || ii==0 || ii==ny-1 ||
		 iii==0 || iii==nz-1 || iiii==0 || iiii==nt-1)
		   fileout.write(" "+x+" "+y+" "+z+" "+t+" "+
				 ub(x,y,z,t)+" \n");
	      else
		fileout.write(" "+x+" "+y+" "+z+" "+t+" "+
			      ug[s(i,ii,iii,iiii)]+" \n");
	    } // end t
 	  } // end z
	} // end y
      } // end x
      fileout.close();
      System.out.println("finished writing pde_bihar44t_eq.dat");
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("pde_bihar44t_eq.dat  not opened");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("pde_bihar44t_eq.java output-file error");
    } // end try
  } // end write_soln

  void write_eqn()
  {
    try
    {
      FileWriter fout = new FileWriter("pde_bihar44t_eq.eqn");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde_bihar44t_eq.dat");
      fileout.write(nxyzt+" "+(nxyzt+1)+"\n");

      for(int i=0; i<nxyzt; i++)
      {
	for(int j=0; j<nxyzt; j++)
	{
	  fileout.write(kg[i][j]+" \n");
	} // end j
	fileout.write(fg[i]+" \n");
      } // end i
      fileout.close();
      System.out.println("finished writing pde_bihar44t_eq.dat");
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("pde_bihar44t_eq.eqn  not opened");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("pde_bihar44t_eq.java output-file error");
    } // end try
  } // end write_eqn

  void write_eqnij()
  {
    try
    {
      FileWriter fout = new FileWriter("pde_bihar44t_eqij.eqn");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde_bihar44t_eqij.dat");
      fileout.write(nxyzt+" "+(nxyzt+1)+"\n");

      for(int i=0; i<nxyzt; i++)
      {
	for(int j=0; j<nxyzt; j++)
	{
	  if(kg[i][j]!=0.0)
	     fileout.write(i+" "+j+" "+kg[i][j]+" \n");
	} // end j
	fileout.write(i+" "+nxyzt+" "+fg[i]+" \n");
      } // end i
      fileout.close();
      System.out.println("finished writing pde_bihar44t_eqij.dat");
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("pde_bihar44t_eqij.eqn  not opened");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("pde_bihar44t_eq.java output-file error");
    } // end try
  } // end write_eqnij

  void check_soln()
                  // for all interior points
  {
    double x, y, z, t, err, maxerr, avgerr, rmserr;
    double sumerr, sumsqerr, r_deriv, a_deriv, eval_deriv;
    double t_deriv = 0.0;
    int cnt = 0;
    int i4;
    boolean check = false; // set for cases below

    System.out.println("check_soln values ");
    sumerr = 0.0;
    sumsqerr = 0.0;
    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // full equation 
    {
      x = xg[i];
      for(int ii=1; ii<ny-1; ii++)
      {
        y = yg[ii];
        for(int iii=1; iii<nz-1; iii++)
        {
          z = zg[iii];
          for(int iiii=1; iiii<nt-1; iiii++)
          {
            t = tg[iiii];
            i4 = s(i,ii,iii,iiii);

            if(debug && i==nx-2 && ii==ny-2 && iii==nz-2 && iiii==nt-2)
               check = true;
            if(check)
	    {
              System.out.println("checking on check against PDE i4="+i4);
              System.out.println("xg[i]="+xg[i]+", yg[ii]="+yg[ii]+
                   ", zg[iii]="+zg[iii]+", tg[iiii]="+tg[iiii]);
	    }

            // PDE 2 u(x,y,z,t) term 
            eval_deriv = 2.0*ug[i4];
            if(check) System.out.println("2.0*ug ="+(2.0*ug[i4])+", err="+
					 (eval_deriv-2.0*ub(x,y,z,t)));
            // t_deriv test derivative based on computed solution
            // a_deriv test derivative based on exact solution

            // PDE uxxxx term
            t_deriv = 0.0;           
            for(int j=1; j<nx-1; j++)
            {
	      t_deriv += cxxxx[i][j]*ug[s(j,ii,iii,iiii)];
              if(check) System.out.println("cxxxx["+i+"]["+j+"]="+
                                           cxxxx[i][j]+", ug[s("+j+
                                  ",ii,iii,iiii)]="+ug[s(j,ii,iii,iiii)]);
            } 
            t_deriv += cxxxx[i][0]*ub(xg[0],y,z,t);
            if(check) System.out.println("cxxxx["+i+"][0]="+
	              cxxxx[i][0]+", ub(xg[0],y,z,t)="+ub(xg[0],y,z,t));
            t_deriv += cxxxx[i][nx-1]*ub(xg[nx-1],y,z,t);
            if(check) System.out.println("cxxxx["+i+"][nx-1]="+
	      cxxxx[i][nx-1]+", ub(xg[nx-1],y,z,t)="+ub(xg[nx-1],y,z,t));
            eval_deriv += t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int j=1; j<nx-1; j++)
              {
		a_deriv += cxxxx[i][j]*ub(xg[j],y,z,t);
              } 
              a_deriv += cxxxx[i][0]*ub(xg[0],y,z,t);
              a_deriv += cxxxx[i][nx-1]*ub(xg[nx-1],y,z,t);
	      r_deriv = 24.0*t;
	      err = r_deriv - t_deriv;
              System.out.println("uxxxx="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("uxxxx="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end uxxxx check

            // PDE uyyyy term 
            t_deriv = 0.0;
            for(int jj=1; jj<ny-1; jj++)
            {
	      t_deriv += cyyyy[ii][jj]*ug[s(i,jj,iii,iiii)];
              if(check) System.out.println("cyyyy["+ii+"]["+jj+"]="+
                                           cyyyy[ii][jj]+", ug[s(i,"+jj+
                                  ",iii,iiii)]="+ug[s(i,jj,iii,iiii)]);
            }
            t_deriv += cyyyy[ii][0]*ub(x,yg[0],z,t);
            if(check) System.out.println("cyyyy["+ii+"][0]="+
	     cyyyy[ii][0]+", ub(x,yg[0],z,t)="+ub(x,yg[0],z,t));
            t_deriv += cyyyy[ii][ny-1]*ub(x,yg[ny-1],z,t);
            if(check) System.out.println("cyyyy["+ii+"][ny-1]="+
	     cyyyy[ii][ny-1]+", ub(x,yg[ny-1],z,t)="+ub(x,yg[ny-1],z,t));
            eval_deriv += t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int jj=1; jj<ny-1; jj++)
              {
		a_deriv += cyyyy[ii][jj]*ub(x,yg[jj],z,t);
              } 
              a_deriv += cyyyy[ii][0]*ub(x,yg[0],z,t);
              a_deriv += cyyyy[ii][ny-1]*ub(x,yg[ny-1],z,t);
	      r_deriv = 24.0*t;
	      err = r_deriv - t_deriv;
              System.out.println("uyyyy="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("uyyyy="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end uyyyy check

            // PDE uzzzz term 
            t_deriv = 0.0;           
            for(int jjj=1; jjj<nz-1; jjj++)
            {
	      t_deriv += czzzz[iii][jjj]*ug[s(i,ii,jjj,iiii)];
              if(check) System.out.println("czzzz["+iii+"]["+jjj+"]="+
                                           czzzz[iii][jjj]+", ug[s(i,ii,"+jjj+
                                  ",iiii)]="+ug[s(i,ii,jjj,iiii)]);
            } 
            t_deriv += czzzz[iii][0]*ub(x,y,zg[0],t);
            if(check) System.out.println("czzzz["+iii+"][0]="+
	     czzzz[iii][0]+", ub(x,y,zg[0],t)="+ub(x,y,zg[0],t));
            t_deriv += czzzz[iii][nz-1]*ub(x,y,zg[nz-1],t);
            if(check) System.out.println("czzzz["+iii+"][nz-1]="+
	     czzzz[iii][nz-1]+", ub(x,y,zg[nz-1],t)="+ub(x,y,zg[nz-1],t));
            eval_deriv += t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int jjj=1; jjj<nz-1; jjj++)
              {
		a_deriv += czzzz[iii][jjj]*ub(x,y,zg[jjj],t);
              } 
              a_deriv += czzzz[iii][0]*ub(x,y,zg[0],t);
              a_deriv += czzzz[iii][nz-1]*ub(x,y,zg[nz-1],t);
	      r_deriv = 24.0*t;
	      err = r_deriv - t_deriv;
              System.out.println("uzzzz="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("uzzzz="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end uzzzz check

            // PDE utttt term 
            t_deriv = 0.0;           
            for(int jjjj=1; jjjj<nt-1; jjjj++)
            {
	      t_deriv += ctttt[iiii][jjjj]*ug[s(i,ii,iii,jjjj)];
              if(check) System.out.println("ctttt["+iiii+"]["+jjjj+"]="+
                                       ctttt[iiii][jjjj]+", ug[s(i,ii,iii,"+
                                  jjjj+")]="+ug[s(i,ii,iii,jjjj)]);
            }
            t_deriv += ctttt[iiii][0]*ub(x,y,z,tg[0]);
            if(check) System.out.println("ctttt["+iiii+"][0]="+
			ctttt[iiii][0]+", ub(x,y,z,tg[0])="+ub(x,y,z,tg[0]));
            t_deriv += ctttt[iiii][nt-1]*ub(x,y,z,tg[nt-1]);
            if(check) System.out.println("ctttt["+iiii+"][nt-1]="+
	      ctttt[iiii][nt-1]+", ub(x,y,z,tg[nz-1])="+ub(x,y,z,tg[nt-1]));
            eval_deriv += t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int jjjj=1; jjjj<nt-1; jjjj++)
              {
		a_deriv += ctttt[iiii][jjjj]*ub(x,y,z,tg[jjjj]);
              } 
              a_deriv += ctttt[iiii][0]*ub(x,y,z,tg[0]);
              a_deriv += ctttt[iiii][nt-1]*ub(x,y,z,tg[nt-1]);
	      r_deriv = 0;
	      err = r_deriv - t_deriv;
              System.out.println("utttt="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("utttt="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end utttt check

            // PDE  2*uxx term 
            t_deriv = 0.0;
            for(int j=1; j<nx-1; j++)
            {
	      t_deriv += cxx[i][j]*ug[s(j,ii,iii,iiii)];
              if(check) System.out.println("cxx["+i+"]["+j+"]="+
                                           cxx[i][j]+", ug[s("+j+
                                  ",ii,iii,iiii)]="+ug[s(j,ii,iii,iiii)]);
            } // end j 
            t_deriv += cxx[i][0]*ub(xg[0],y,z,t);
            if(check) System.out.println("cxx["+i+"][0]="+
	              cxx[i][0]+", ub(xg[0],y,z,t)="+ub(xg[0],y,z,t));
            t_deriv += cxx[i][nx-1]*ub(xg[nx-1],y,z,t);
            if(check) System.out.println("cxx["+i+"][nx-1]="+
	      cxx[i][nx-1]+", ub(xg[nx-1],y,z,t)="+ub(xg[nx-1],y,z,t));
            eval_deriv += 2.0*t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int j=1; j<nx-1; j++)
              {
		a_deriv += cxx[i][j]*ub(xg[j],y,z,t);
              } 
              a_deriv += cxx[i][0]*ub(xg[0],y,z,t);
              a_deriv += cxx[i][nx-1]*ub(xg[nx-1],y,z,t);
	      r_deriv = 12.0*x*x*t;
	      err = r_deriv - t_deriv;
              System.out.println("uxx="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("uxx="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end uxx check

            // PDE  2*uyy term 
            t_deriv = 0.0;           
            for(int jj=1; jj<ny-1; jj++)
            {
              t_deriv += cyy[ii][jj]*ug[s(i,jj,iii,iiii)];
            } // end jj 
            t_deriv += cyy[ii][0]*ub(x,yg[0],z,t);
            t_deriv += cyy[ii][ny-1]*ub(x,yg[ny-1],z,t);
            eval_deriv += 2.0*t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int jj=1; jj<ny-1; jj++)
              {
	        a_deriv += cyy[ii][jj]*ub(x,yg[jj],z,t);
              } 
              a_deriv += cyy[ii][0]*ub(x,yg[0],z,t);
              a_deriv += cyy[ii][ny-1]*ub(x,yg[ny-1],z,t);
	      r_deriv = 12.0*(y-z)*(y-z)*t;
	      err = r_deriv - t_deriv;
              System.out.println("uyy="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("uyy="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end uyy check

            // PDE  2*uzz term 
            t_deriv = 0.0;           
            for(int jjj=1; jjj<nz-1; jjj++)
            {
              t_deriv += czz[iii][jjj]*ug[s(i,ii,jjj,iiii)];
            } // end jjj 
            t_deriv += czz[iii][0]*ub(x,y,zg[0],t);
            t_deriv += czz[iii][nz-1]*ub(x,y,zg[nz-1],t);
            eval_deriv += 2.0*t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int jjj=1; jjj<nz-1; jjj++)
              {
		  a_deriv += czz[iii][jjj]*ub(x,y,zg[jjj],t);
              } 
              a_deriv += czz[iii][0]*ub(x,y,zg[0],t);
              a_deriv += czz[iii][nz-1]*ub(x,y,zg[nz-1],t);
	      r_deriv = 12.0*(y-z)*(y-z)*t;
	      err = r_deriv - t_deriv;
              System.out.println("uzz="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("uzz="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end uzz check

            // PDE  2*utt term 
            t_deriv = 0.0;           
            for(int jjjj=1; jjjj<nt-1; jjjj++)
            {
              t_deriv += ctt[iiii][jjjj]*ug[s(i,ii,iii,jjjj)];
            } // end jjjj 
            t_deriv += ctt[iiii][0]*ub(x,y,z,tg[0]);
            t_deriv += ctt[iiii][nt-1]*ub(x,y,z,tg[nt-1]);
            eval_deriv += 2.0*t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int jjjj=1; jjjj<nt-1; jjjj++)
              {
		  a_deriv += ctt[iiii][jjjj]*ub(x,y,z,tg[jjjj]);
              } 
              a_deriv += ctt[iiii][0]*ub(x,y,z,tg[0]);
              a_deriv += ctt[iiii][nt-1]*ub(x,y,z,tg[nt-1]);
	      r_deriv = 0.0;
	      err = r_deriv - t_deriv;
              System.out.println("utt="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("utt="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end utt check

            err = eval_deriv-F(x,y,z,t);
            if(check)
	    {
	      System.out.println("eval_deriv="+eval_deriv+", F(x,y,z,t)="+F(x,y,z,t));
	      System.out.println("err="+err);
	    }
            err = Math.abs(err);
            sumerr += err;
            sumsqerr += err*err;
            maxerr = Math.max(err,maxerr);
            cnt++;
	  } // end t 
        } // end z  
      } // end y 
    } // end x 
    rmserr = Math.sqrt(sumsqerr/(double)cnt);
    avgerr = sumerr/(double)cnt;
    System.out.println("maxerr="+maxerr+", rmserr="+rmserr+", avgerr="+avgerr);  
  } // end check_soln 

  void rderiv(int order, int npts, int point, double h, double c[])
  {
    nderiv D = new nderiv();
    int a[] = new int[npts]; // basic discretization from nderiv
    int b[] = new int[1];    // out parameter, denominator

    D.deriv(order,npts,point,a,b);
    for(int i=0; i<npts; i++)
    {
      c[i] = (double)a[i]/((double)b[0]);
      for(int j=0; j<order; j++) c[i] = c[i]/h;
    }
  } // end rderiv

  public static void main (String[] args)
  {
    pde_bihar44t_eq F2 = new pde_bihar44t_eq();
  
    F2.run(-1.0, 1.0, -1.0, 1.0, -0.2, 0.2, 0.0, 2.0, 5, 5, 5, 5);
    F2.run(-1.0, 1.0, -1.0, 1.0, -0.2, 0.2, 0.0, 2.0, 7, 7, 7, 7);
          // xmin, xmax, ymin, ymax, zmin, zmax, tmin, tmax,
          // nx, ny, nz, nt  
    F2.run(-1.0, 1.0, -1.0, 1.0, -0.2, 0.2, 0.0, 2.0, 7, 5, 5, 5);
    F2.run(-1.0, 1.0, -1.0, 1.0, -0.2, 0.2, 0.0, 2.0, 5, 7, 5, 5);
    F2.run(-1.0, 1.0, -1.0, 1.0, -0.2, 0.2, 0.0, 2.0, 5, 5, 7, 5);
    F2.run(-1.0, 1.0, -1.0, 1.0, -0.2, 0.2, 0.0, 2.0, 5, 5, 5, 7);
    try
    {
      Runtime runtime = Runtime.getRuntime();
      Process proc = runtime.exec("java -cp . plot4d < pde_bihar44t_eq.dat ");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("pde_bihar44t_eq.java runtime error");
    } // end try
  }
} // end class pde_bihar44t_eq