// pde44h_eq.java  homogenous Biharmonic 4th order, 4 dimensions
// 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 u(x,y,z,t) = 0
//
// boundary conditions computed using u(x,y,z,t)
// analytic solution is u(x,y,z,t) = sim(x,y,z,t)
//
// replace continuous derivatives with discrete derivatives
// solve for Uijkl numerical approximation U(x_i,y_j,z_k,t_l)
// A * U = F, solve simultaneous equations for U
//
// fourth order numerical difference order=4, npoints=5, term=0
// d^4u/dx^4 = (1/(b4[0]*hx^4))*(a4[0][0]*u(x,y,z,t) + a4[0][1]*u(x+hx,y,z,t) +
// a4[0][2]*u(x+2hx,y,z,t) + a4[0][3]*u(x+3hx,y,z,t) +a4[0][4]*u(x+4hx,y,z,t))
//
// d^4u/dx^4 using subscripts i,j,k,l for i a minimum subscript of x
// d^4u/dx^4 =(1/(b4[0]*hx^4))*(a4[0][0]*u(i,j,k,l) + a4[0][1]*u(i+1,j,k,l) +
// a4[0][2]*u(i+2,j,k,l) + a4[0][3]*u(i+3,j,k,l) +a4[0][4]*u(i+4,j,k,l))
//
// d^4u/dy^4 using subscripts i,j,k,l for j a minimum subscript of y, etc 
// d^4u/dy^4 =(1/(b4[0]*hy^4))*(a4[0][0]*u(i,j,k,l) + a4[0][1]*u(i,j+1,k,l) +
// a4[0][2]*u(i,j+2,k,l) + a4[0][3]*u(i,j+3,k,l) +a4[0][4]*u(i,j+4,k,l))
//
// d^4u/dx^4 using subscripts i,j,k,l for i-2 and i+2 allowed subscripts of x
// d^4u/dx^4 =(1/(b4[2]*hx^4))*(a4[2][0]*u(i-2,j,k,l) + a4[2][1]*u(i-1,j,k,l) +
// a4[2][2]*u(i,j,k,l) + a4[2][3]*u(i+1,j,k,l) +a4[2][4]*u(i+2,j,k,l))
//
// d^4u/dx^4 using subscripts i,j,k,l for i a maximum subscript of x
// d^4u/dx^4 =(1/(b4[4]*hx^4))*(a4[4][0]*u(i-4,j,k,l) + a4[4][1]*u(i-3,j,k,l) +
// a4[4][2]*u(i-2,j,k,l) + a4[4][3]*u(i-1,j,k,l) +a4[4][4]*u(i,j,k,l))
//
// subscripts in code use x[i], y[ii], z[iii], t[iiii]
// derivatives are computed with b,hx,a included as in
// cx[j] for first derivative of x at point j
// cxxxx[j] for fourth derivative of x at point j
// cyy[jj]*czz[jjj] for d^4u/dy^2dz^2 at jj, jjj
 
// The subscript versions of the derivatives are substituted into the
// equation to be solved. The resulting equation is analytically solved
// for u(i,j,k,l) = other u terms with coefficients and + f(xi,yj,zk,tl)
// The boundaries xmin..xmax, ymin..ymax, zmin..zmax, tmin..tmax are known
// The number of grid points in each dimension is known nx, ny, nz, nt
// hx=(xmax-xmin)/(nx-1), hy=(ymax-ymin)/(ny-1),
// hz=(zmax-zmin)/(nz-1), ht=(tmax-tmin)/(tx-1), 
// thus xi=xmin+i*hx, yj=ymin+j*hy, zk=zmin+k*hz, tl=tmin+l*ht
// then the linear system of equations is built,
// nxyzt=(nx-2)*(ny-2)*(nz-2)*(nt-2) equations because the boundary value
//                                   is known on each end.
// The matrix A is nxyzt rows by nxyzt columns (nxyzt+1) columns including F
// The forcing function vector F is nxyzt elements adjunct to A
// The solution vector U is nxyzt elements.
// The building of the A matix requires 4 nested loops for rows
// i  in 1..nx-1, j  in 1..ny-1,  k  in 1..nz-1, l  in 1..nt-1 for rows of A
// then each term in the differential equation fills in that row
// and cs = (nx-2)*(ny-2)*(nz-2)*(nt-2) the RHS, fills the last column of A
// subscripting functions are used to compute linear subscripts of A, U
// row=(i-1)*(ny-2)*(nz-2)*(nt-2)+(ii-1)*(nz-2)*(nt-2)+(iii-1)*(nt-2)+(iiii-1)
// row*(nxyzt+1)+col for A
//
// care must be taken to move the negative of boundary elements to
// the right hand size of the equation to be solved. These are subtracted
// from the computed value of the forcing function, cs, for that row.
// tests are j==0 || j==nx-1   ...  jjjj==0 || jjjj==nt-1 for boundaries
 
// uses simeq.java
// uses deriv.java
import java.text.*;
import java.io.*;

class pde44h_eq
{
  final int nx = 6; // problem parameters
  final int ny = 6;
  final int nz = 6;
  final int nt = 6;
  final int nxyzt = (nx-2)*(ny-2)*(nz-2)*(nt-2);

  double A[][] = new double[nxyzt][nxyzt+1]; // last column is RHS 
  double U[] = new double[nxyzt]; // computed solution
  double xg[] = new double[nx];
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double tg[] = new double[nt];
  double Ua[] = new double[nxyzt]; // check solution
  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 hx, hy, hz, ht;

  DecimalFormat f6 = new DecimalFormat("0.000");
  DecimalFormat fe = new DecimalFormat("0.000E00");

  pde44h_eq()
  {
    double x, y, z, t;
    int k, cs;
    // derivative point arrays, dynamic in this version 
    double cxx[] = new double[nx];
    double cyy[] = new double[ny];
    double czz[] = new double[nz];
    double ctt[] = new double[nt]; 
    double cxxxx[] = new double[nx];
    double cyyyy[] = new double[ny];
    double czzzz[] = new double[nz];
    double ctttt[] = new double[nt]; 
    double val, err, avgerr, maxerr;
    double time_start;
    double time_now;
    double time_last;

    System.out.println("pde44h_eq.java running ");
    System.out.println("Given uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) + ");
    System.out.println("    utttt(x,y,z,t) + 2 uxx(x,y,z,t) + 2 uyy(x,y,z,t) + ");
    System.out.println("    2 uzz(x,y,z,t) + 2 u(x,y,z,t) = 0");
    System.out.println("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries ");
    System.out.println("Analytic solution u(x,y,z,t) =  sin(x+y+z+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);
    time_start = (double)System.currentTimeMillis()/1000.0;
    time_last = time_start;

    System.out.println("x grid and analytic solution at ymin,zmin,tmin ");
    hx = (xmax-xmin)/(double)(nx-1);
    for(int i=0; i<nx; i++)
    {
      xg[i] = xmin+(double)i*hx;
      Ua[i] = uana(xg[i],ymin,zmin,tmin); // temp 
      System.out.println("i="+i+", Ua("+f6.format(xg[i])+","+
                         f6.format(ymin)+","+f6.format(zmin)+","+
                         f6.format(tmin)+")="+f6.format(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(int ii=0; ii<ny; ii++)
    {
      yg[ii] = ymin+(double)ii*hy;
      Ua[ii] = uana(xmin, yg[ii],zmin,tmin); // temp 
      System.out.println("ii="+ii+", Ua("+f6.format(xmin)+","+
                         f6.format(yg[ii])+","+f6.format(zmin)+","+
                         f6.format(tmin)+")="+f6.format(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(int iii=0; iii<nz; iii++)
    {
      zg[iii] = zmin+(double)iii*hz;
      Ua[iii] = uana(xmin, ymin, zg[iii],tmin); // temp 
      System.out.println("iii="+iii+", Ua("+f6.format(xmin)+","+
                         f6.format(ymin)+","+f6.format(zg[iii])+","+
                         f6.format(tmin)+")="+f6.format(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(int iiii=0; iiii<nt; iiii++)
    {
      tg[iiii] = tmin+(double)iiii*ht;
      Ua[iiii] = uana(xmin, ymin, zmin, tg[iiii]); // temp 
      System.out.println("iiii="+iiii+", Ua("+f6.format(xmin)+","+
                         f6.format(ymin)+","+f6.format(zmin)+","+
                         f6.format(tg[iiii])+")="+f6.format(Ua[iiii]));
    } // iiii 
    System.out.println(" ");

    System.out.println("check rderiv(4,nx,0,hx,cxxxx");
    new rderiv(4,nx,0,hx,cxxxx);
    for(int i=0; i<nx; i++) System.out.println("cxxxx["+i+"]="+cxxxx[i]);
    System.out.println("check rderiv(4,nx,nx-1,hx,cxxxx");
    new rderiv(4,nx,nx-1,hx,cxxxx);
    for(int i=0; i<nx; i++) System.out.println("cxxxx["+i+"]="+cxxxx[i]);
    System.out.println(" ");

    System.out.println("put solution in Ua vector ");
    // 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)] = uana(x,y,z,t);
          } // iiii 
        } // iii 
      } // ii 
    } // i 
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("time for previous section = "+(time_now-time_last)+
                       " seconds ");
    time_last = time_now;
    System.out.println(" ");

    System.out.println("initialize A ");
    // initialize A 
    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++)
          {
            for(int j=1; j<nx-1; j++)
            {
              for(int jj=1; jj<ny-1; jj++)
              {
                for(int jjj=1; jjj<nz-1; jjj++)
                {
                  for(int jjjj=1; jjjj<nt-1; jjjj++)
                  {
		    A[s(i,ii,iii,iiii)][s(j,jj,jjj,jjjj)] = 0.0;
	          } // jjjj 
	        } // jjj 
	      } // jj 
            } // j 
          } // iiii 
        } // iii 
      } // ii 
    } // i 

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

    // compute entries in A matrix, inside boundary 
    System.out.println("compute A matrix ");
    cs = (nx-2)*(ny-2)*(nz-2)*(nt-2); // constant RHS column 
    val = 0.0;
    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++)
          {
            k = s(i,ii,iii,iiii); // row index 

            // for each term in first equation 

	    // for d^4u/dx^4 
            new rderiv(4,nx,i,hx,cxxxx);
            for(int j=0; j<nx; j++)
            {
              val = cxxxx[j];
              if(j==0 || j==nx-1)
	           A[k][cs] -= val*uana(xg[j],yg[ii],zg[iii],tg[iiii]);
              else A[k][s(j,ii,iii,iiii)] += val;
	    }

	    // for d^4u/dy^4 
            new rderiv(4,ny,ii,hy,cyyyy);
            for(int jj=0; jj<ny; jj++)
            {
              val = cyyyy[jj];
              if(jj==0 || jj==ny-1)
	           A[k][cs] -=  val*uana(xg[i],yg[jj],zg[iii],tg[iiii]);
              else A[k][s(i,jj,iii,iiii)] += val;
	    }

	    // for d^4u/dz^4 
            new rderiv(4,nz,iii,hz,czzzz);
            for(int jjj=0; jjj<nz; jjj++)
            {
              val = czzzz[jjj];
              if(jjj==0 || jjj==nz-1)
	           A[k][cs] -= val*uana(xg[i],yg[ii],zg[jjj],tg[iiii]);
              else A[k][s(i,ii,jjj,iiii)] += val;
	    }

	    // for d^4u/dt^4 
            new rderiv(4,nt,iiii,ht,ctttt);
            for(int jjjj=0; jjjj<nt; jjjj++)
            {
              val = ctttt[jjjj];
              if(jjjj==0 || jjjj==nt-1)
	           A[k][cs] -= val*uana(xg[i],yg[ii],zg[iii],tg[jjjj]);
              else A[k][s(i,ii,iii,jjjj)] += val;
	    }

	    // for d^2u/dx^2 
            new rderiv(2,nx,i,hx,cxx);
            for(int j=0; j<nx; j++)
            {
              val = cxx[j];
              if(j==0 || j==nx-1)
	           A[k][cs] -= val*uana(xg[j],yg[ii],zg[iii],tg[iiii]);
              else A[k][s(j,ii,iii,iiii)] += val;
	    }

	    // for d^2u/dy^2 
            new rderiv(2,ny,ii,hy,cyy);
            for(int jj=0; jj<ny; jj++)
            {
              val = cyy[jj];
              if(jj==0 || jj==ny-1)
	           A[k][cs] -=  val*uana(xg[i],yg[jj],zg[iii],tg[iiii]);
              else A[k][s(i,jj,iii,iiii)] += val;
	    }

	    // for d^2u/dz^2 
            new rderiv(2,nz,iii,hz,czz);
            for(int jjj=0; jjj<nz; jjj++)
            {
              val = czz[jjj];
              if(jjj==0 || jjj==nz-1)
	           A[k][cs] -= val*uana(xg[i],yg[ii],zg[jjj],tg[iiii]);
              else A[k][s(i,ii,jjj,iiii)] += val;
	    }

	    // for u(x,y,z,t) 
            A[k][s(i,ii,iii,iiii)] += 2.0;

            // end terms 

	  } // end iiii loop 
        } // end iii loop 
      } // end ii loop 
    } // end i loop 

    // printA();

    System.out.println("A computed stiffness matrix ");
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("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, A, U);

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

    System.out.println("check PDE against known solution");
    check_soln(Ua);
    System.out.println("check PDE against computed solution");
    check_soln(U);
    System.out.println(" ");

    write_soln();

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("          U 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++)
          {
            err = Math.abs(U[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+"]="+
			       f6.format(U[s(i,ii,iii,iiii)])+", Ua="+
                               f6.format(Ua[s(i,ii,iii,iiii)])+", err="+
                               fe.format(err));
          } // iiii 
        } // iii 
      } // ii 
    } // i 

    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("total CPU time  = "+(time_now-time_start)+" seconds"); 
    System.out.println(" ");


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

  int s(int i, int ii, int iii, int iiii) // for x,y,z,t 
  {
    return (i-1)*(ny-2)*(nz-2)*(nt-2) + (ii-1)*(nz-2)*(nt-2) +
           (iii-1)*(nt-2) + (iiii-1);
  } // end s 

  int sk(int k, int i, int ii, int iii, int iiii) // for x,y,z,t 
  {
    return k*(nxyzt+1) + (i-1)*(ny-2)*(nz-2)*(nt-2) + (ii-1)*(nz-2)*(nt-2) +
           (iii-1)*(nt-2) + (iiii-1);
  } // end sk 

  int ss(int i, int ii, int iii, int iiii, int j, int jj, int jjj, int jjjj)
  {
    return s(i,ii,iii,iiii)*(nxyzt+1) + s(j,jj,jjj,jjjj);
  } // end ss 

  // problem RHS, 0

  // problem boundaries (and analytic solution for checking)
  double uana(double x, double y, double z, double t)
  {
      return Math.sin(x+y+z+t);
  }

  void printA()
  {
    int cs, k;
    cs = (nx-2)*(ny-2)*(nz-2)*(nt-2); // constant RHS column 
    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++)
          {
            k = s(i,ii,iii,iiii); // row index 
            for(int j=1; j<nx-1; j++)
            {
              for(int jj=1; jj<ny-1; jj++)
              {
                for(int jjj=1; jjj<nz-1; jjj++)
                {
                  for(int jjjj=1; jjjj<nt-1; jjjj++)
                  {
                    System.out.println("A[row("+i+","+ii+","+iii+","+iiii+
                                   "),col("+j+","+jj+","+jjj+","+jjjj+")] ="+
			           A[s(i,ii,iii,iiii)][s(j,jj,jjj,jjjj)]);
	          } // jjjj 
	        } // jjj 
	      } // jj 
            } // j 
            System.out.println("RHS A["+k+","+cs+"]="+A[k][cs]);
            // equivalent A[sk(k,i,ii,iii,iiii+1)] 
          } // iiii 
        } // iii 
      } // ii 
    } // i 
    System.out.println(" ");
  } // end printA 

  void check_soln(double u[])
  {
    // 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 u(x,y,z,t) = 0
    double val, err, smaxerr;

    smaxerr = 0.0;
    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;
            // uxxxx(x,y,z,t)
            val += eval_derivative(4, 0, 0, 0, i, ii, iii, iiii, u);
	    // + uyyyy(x,y,z,t)
            val += eval_derivative(0, 4, 0, 0, i, ii, iii, iiii, u);
	    // + uzzzz(x,y,z,t)
            val += eval_derivative(0, 0, 4, 0, i, ii, iii, iiii, u);
            // + utttt(x,y,z,t)
            val += eval_derivative(0, 0, 0, 4, i, ii, iii, iiii, u);
	    // + 2  uxx(x,y,z,t)
            val += 2.0*eval_derivative(2, 0, 0, 0, i, ii, iii, iiii, u);
	    // + 2 uyy(x,y,z,t)
            val += 2.0*eval_derivative(0, 2, 0, 0, i, ii, iii, iiii, u);
            // + 2 uzz(x,y,z,t)
            val += 2.0*eval_derivative(0, 0, 2, 0, i, ii, iii, iiii, u);
	    // + 2 u(x,y,z,t)
            val += 2.0*uana(xg[i],yg[ii],zg[iii],tg[iiii]);
            //        F(x,y,z,t) should be zero
	    err = Math.abs(val);
            if(err>smaxerr) smaxerr = err;
	  } // end iiii
	} // end iii
      } // end ii
    } // end i
    System.out.println("check_soln maxerr="+smaxerr);
  } // end check_soln

  void write_soln()
  {
    try
    {
      FileWriter fout = new FileWriter("pde44h_eq_java.dat");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde44h_eq.dat");
      // not including boundary
      fileout.write("4 "+(nx-2)+" "+(ny-2)+" "+(nz-2)+" "+(nt-2)+"\n");
      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++)
            {
              fileout.write(" "+xg[i]+" "+yg[ii]+" "+zg[iii]+" "+
                            tg[iiii]+" "+U[s(i,ii,iii,iiii)]+" \n");
	    } // end iiii
	  } // end iii
        } // end ii
      } // end i

      fileout.close();
      System.out.println("finished writing pde44h_eq_java.dat");
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("pde44h_eq_java.dat  not opened");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("pde44h_eq_java.java output-file error");
    } // end try
  } // end write_soln

  double eval_derivative(int xord, int yord, int zord, int tord,
                         int i, int ii, int iii, int iiii, double u[])
  {
    double cx[] = new double[nx];
    double cy[] = new double[ny];
    double cz[] = new double[nz];
    double ct[] = new double[nt];
    int nn = Math.max(nx, Math.max(ny, Math.max(nz, nt)));
    double p1[] = new double[nn];
    double p2[] = new double[nn];
    double p3[] = new double[nn];
    double discrete;
    double x, y, z, t;

    new rderiv(xord, nx, i, hx, cx);
    x = xg[i];
    new rderiv(yord, ny, ii, hy, cy);
    y = yg[ii];
    new rderiv(zord, nz, iii, hz, cz);
    z = zg[iii];
    new rderiv(tord, nt, iiii, ht, ct);
    t = tg[iiii];
    discrete = 0.0;

    // four cases of one partial x, y, z, t to any order
    if(xord!=0 && yord==0 && zord==0 && tord==0)
      for(int j=0; j<nx; j++)
	if(j==0 || j==nx-1) discrete += cx[j]*uana(xg[j],y,z,t);
        else discrete += cx[j]*u[s(j,ii,iii,iiii)];
    else if(xord==0 && yord!=0 && zord==0 && tord==0)
      for(int jj=0; jj<ny; jj++)
	if(jj==0 || jj==ny-1) discrete += cy[jj]*uana(x,yg[jj],z,t);
        else discrete += cy[jj]*u[s(i,jj,iii,iiii)];
    else if(xord==0 && yord==0 && zord!=0 && tord==0)
      for(int jjj=0; jjj<nz; jjj++)
	if(jjj==0 || jjj==nz-1) discrete += cz[jjj]*uana(x,y,zg[jjj],t);
        else discrete += cz[jjj]*u[s(i,ii,jjj,iiii)];
    else if(xord==0 && yord==0 && zord==0 && tord!=0)
      for(int jjjj=0; jjjj<nt; jjjj++)
	if(jjjj==0 || jjjj==nt-1) discrete += ct[jjjj]*uana(x,y,z,tg[jjjj]);
        else discrete += ct[jjjj]*u[s(i,ii,iii,jjjj)];
    // six cases of two partials xy, xz, xt, yz, yt, zt to any order 
    else if(xord!=0 && yord!=0 && zord==0 && tord==0)
    {
      for(int j=0; j<nx; j++)
      {
        p1[j] = 0.0;
        for(int jj=0; jj<ny; jj++)
	  if(j==0 || j==nx-1 || jj==0 || jj==ny-1)
	       p1[j] += cy[jj]*uana(xg[j],yg[jj],zg[iii],tg[iiii]);
          else p1[j] += cy[jj]*u[s(j,jj,iii,iiii)];
        discrete += cx[j]*p1[j];
      } // end j
    }
    else if(xord!=0 && yord==0 && zord!=0 && tord==0)
    {
      for(int j=0; j<nx; j++)
      {
        p1[j] = 0.0;
        for(int jjj=0; jjj<nz; jjj++)
	  if(j==0 || j==nx-1 || jjj==0 || jjj==nz-1)
	       p1[j] += ct[jjj]*uana(xg[j],yg[ii],zg[jjj],tg[iiii]);
          else p1[j] += cz[jjj]*u[s(j,ii,jjj,iiii)];
        discrete += cx[j]*p1[j];
      }
    }
    else if(xord!=0 && yord==0 && zord==0 && tord!=0)
    {
      for(int j=0; j<nx; j++)
      {
        p1[j] = 0.0;
        for(int jjjj=0; jjjj<nt; jjjj++)
	  if(j==0 || j==nx-1 || jjjj==0 || jjjj==nt-1)
	       p1[j] += ct[jjjj]*uana(xg[j],yg[ii],zg[iii],tg[jjjj]);
          else p1[j] += ct[jjjj]*u[s(j,ii,iii,jjjj)];
        discrete += cx[j]*p1[j];
      }
    }
    else if(xord==0 && yord!=0 && zord!=0 && tord==0)
    {
      for(int jj=0; jj<ny; jj++)
      {
        p1[jj] = 0.0;
        for(int jjj=0; jjj<nz; jjj++)
	  if(jj==0 || jj==ny-1 || jjj==0 || jjj==nz-1)
	       p1[jj] += cz[jjj]*uana(xg[i],yg[jj],zg[jjj],tg[iiii]);
	  else p1[jj] += cz[jjj]*u[s(i,jj,jjj,iiii)];
        discrete += cy[jj]*p1[jj];
      }
    }
    else if(xord==0 && yord!=0 && zord==0 && tord!=0)
    {
      for(int jj=0; jj<ny; jj++)
      {
        p1[jj] = 0.0;
        for(int jjjj=0; jjjj<nt; jjjj++)
	  if(jj==0 || jj==ny-1 || jjjj==0 || jjjj==nt-1)
	       p1[jj] += ct[jjjj]*uana(xg[i],yg[jj],zg[iii],tg[jjjj]);
          else p1[jj] += ct[jjjj]*u[s(i,jj,iii,jjjj)];
        discrete += cy[jj]*p1[jj];
      }
    }
    else if(xord==0 && yord==0 && zord!=0 && tord!=0)
    {
      for(int jjj=0; jjj<nz; jjj++)
      {
        p1[jjj] = 0.0;
        for(int jjjj=0; jjjj<nt; jjjj++)
	  if(jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1)
               p1[jjj] += ct[jjjj]*uana(xg[i],yg[ii],zg[jjj],tg[jjjj]);
          else p1[jjj] += ct[jjjj]*u[s(i,ii,jjj,jjjj)];
        discrete += cz[jjj]*p1[jjj];
      }
    }
    // four cases of three partials xyz, xyt, xzt, yzt to any order 
    else if(xord!=0 && yord!=0 && zord!=0 && tord==0)
    {
      for(int j=0; j<nx; j++)
      {
        p1[j] = 0.0;
        for(int jj=0; jj<ny; jj++)
	{
          p2[jj] = 0.0;
          for(int jjj=0; jjj<nz; jjj++)
	    if(j==0 || j==nx-1 ||
               jj==0 || jj==ny-1 ||
               jjj==0 || jjj==nz-1)
                 p2[jj] += cz[jjj]*uana(xg[j],yg[jj],zg[jjj],tg[iiii]);
            else p2[jj] += cz[jjj]*u[s(j,jj,jjj,iiii)];
            p1[j] += cy[jj]*p2[jj];
	  }
        discrete += cx[j]*p1[j];
      }
    }
    else if(xord!=0 && yord!=0 && zord==0 && tord!=0)
    {
      for(int j=0; j<nx; j++)
      {
        p1[j] = 0.0;
        for(int jj=0; jj<ny; jj++)
	{
          p2[jj] = 0.0;
          for(int jjjj=0; jjjj<nt; jjjj++)
	    if(j==0 || j==nx-1 ||
               jj==0 || jj==ny-1 ||
               jjjj==0 || jjjj==nt-1)
                 p2[jj] += ct[jjjj]*uana(xg[j],yg[jj],zg[iii],tg[jjjj]);
            else p2[jj] += ct[jjjj]*u[s(j,jj,iii,jjjj)];
          p1[j] += cy[jj]*p2[jj];
        }
        discrete += cx[j]*p1[j];
      }
    }
    else if(xord!=0 && yord==0 && zord!=0 && tord!=0)
    {
      for(int j=0; j<nx; j++)
      {
        p1[j] = 0.0;
        for(int jjj=0; jjj<nz; jjj++)
	{
          p2[jjj] = 0.0;
          for(int jjjj=0; jjjj<nt; jjjj++)
	    if(j==0 || j==nx-1 ||
               jjj==0 || jjj==nz-1 ||
               jjjj==0 || jjjj==nt-1)
                 p2[jjj] += ct[jjjj]*uana(xg[j],yg[ii],zg[jjj],tg[jjjj]);
            else p2[jjj] += ct[jjjj]*u[s(j,ii,jjj,jjjj)];
          p1[j] += cz[jjj]*p2[jjj];
	}
        discrete += cx[j]*p1[j];
      }
    }
    else if(xord==0 && yord!=0 && zord!=0 && tord!=0)
    {
      for(int jj=0; jj<ny; jj++)
      {
        p1[jj] = 0.0;
        for(int jjj=0; jjj<nz; jjj++)
	{
          p2[jjj] = 0.0;
          for(int jjjj=0; jjjj<nt; jjjj++)
	    if(jj==0 || jj==ny-1 ||
               jjj==0 || jjj==nz-1 ||
               jjjj==0 || jjjj==nt-1)
                 p2[jjj] += ct[jjjj]*uana(xg[i],yg[jj],zg[jjj],tg[jjjj]);
            else p2[jjj] += ct[jjjj]*u[s(i,jj,jjj,jjjj)];
          p1[jj] += cz[jjj]*p2[jjj];
	} // end jjj
        discrete += cy[jj]*p1[jj];
      } // end jj
    }
    // one case of four partials
    else
    {
      for(int j=0; j<nx; j++)
      {
        p1[j] = 0.0;
        for(int jj=0; jj<ny; jj++)
	{
          p2[jj] = 0.0;
          for(int jjj=0; jjj<nz; jjj++)
	  {
            p3[jjj] = 0.0;
            for(int jjjj=0; jjjj<nt; jjjj++)
		if(j==0 || j==nx-1 ||
                   jj==0 || jj==ny-1 ||
                   jjj==0 || jjj==nz-1 ||
                   jjjj==0 || jjjj==nt-1)
                   p3[jjj] += ct[jjjj]*uana(xg[j],yg[jj],zg[jjj],tg[jjjj]);
	      else p3[jjj] += ct[jjjj]*u[s(j,jj,jjj,jjjj)];
            p2[jj] += cz[jjj]*p3[jjj];
	  } // end jjj
          p1[j] += cy[jj]*p2[jj];
	} // end jj
        discrete += cx[j]*p1[j];
      } // end j
    } // end if
    return discrete;
  } // end eval_deriv

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




//  modify for reeading nx,xg, ny,yg, nz,zg, nt,tg, bxyzt, bijkl
//  int do_input() // of  .tri  .coord  .bound
//  {
//    String fname;
//    int index, last, len;
//    // read in triangles, set ntri 
//    fname = root+".tri";
//    if(debug>2) System.out.println("about to read triangles from "+fname);
//    ntri=0;
//    minvert=999999;
//    nvert=0;
//
//    try
//    {
//      BufferedReader in = new BufferedReader(new FileReader(fname));
//      String input_line;
//      input_line = in.readLine();
//      while(input_line != null)
//      {
//	 if(debug>1) System.out.println(input_line);
//         // have a line, extract 3 integers
//         len = input_line.length();
//         if(len<5)  // not enough data, ignore blank lines
//	 {
//	   input_line = in.readLine();
//           continue;
//         }
//         index = 0;
//         while(input_line.substring(index,index+1).equals(" ")){index++;}
//         last = input_line.indexOf(' ',index+1);
//         t1[ntri] = Integer.parseInt(input_line.substring(index,last)); 
//         index = last+1;
//         while(input_line.substring(index,index+1).equals(" ")){index++;}
//         last = input_line.indexOf(' ',index+1);
//         t2[ntri] = Integer.parseInt(input_line.substring(index,last)); 
//         index = last+1;
//         while(input_line.substring(index,index+1).equals(" ")){index++;}
//         last = input_line.indexOf(' ',index+1);
//         if(last<1) last = len;
//         t3[ntri] = Integer.parseInt(input_line.substring(index,last)); 
//
//         if(debug>1) System.out.println("tri "+ntri+" has vertices "+t1[ntri]+
//                                        "  "+t2[ntri]+"  "+t3[ntri]);
//         if(t1[ntri]>nvert) nvert=t1[ntri];
//         if(t2[ntri]>nvert) nvert=t2[ntri];
//         if(t3[ntri]>nvert) nvert=t3[ntri];
//         if(t1[ntri]<minvert) minvert=t1[ntri];
//         if(t2[ntri]<minvert) minvert=t2[ntri];
//         if(t3[ntri]<minvert) minvert=t3[ntri];
//         ntri++;
//	 input_line = in.readLine();
//      }
//    }
//    catch(IOException exception)
//    {
//      System.out.println("can not open for reading "+fname);
//      System.out.println(exception);
//      return 1;
//    }
//    return 0;
//  } // end do_input 
//