// pde23js_eq.java converted from
// pde23js_eq.js  2D uniform grid boundary value PDE 
//                known solution to test method
// 
// solve equation in uana23.js
// uana23(x,y,z), f(x,y,z), coefs(coef)
//
// set up and solve system of (nx-2)*(ny-2)*(nz-2) linear equations
//
// uses nuderiv.js
// uses invert.js
// uses simeq.js
// uses uana23.js  uana()  f()  coefs()  define problem 

class pde23js_eq
{
pde23js_eq()
{
  double start, now;
  boolean debug = false;
  double xg[] = new double[30];   // data points
  double yg[] = new double[30];
  double zg[] = new double[30];
  int neqn; 
  int cs;
  double us[] = new double[900];   // solution being computed 
  double ut[][] = new double[900][901]; // [neqn][neqn+1];  equations incl RHS 
  double coef[] = new double[19];
  double ccxx, ccyy, cczz, ccxy, ccxz, ccyz, ccx, ccy, ccz, ccc;
  double xmin, xmax, hx;
  int nx;
  double ymin, ymax, hy;
  int ny;
  double zmin, zmax, hz;
  int nz;
  double error;         // sum of absolute exact-computed         
  double avgerror;      // average error                          
  double maxerror;      // maximum cell error                     

  System.out.print("pde23js_eq.java running \n");
  System.out.print("differential equation to solve, from uana23.js \n");
  System.out.print("given f(x,y,z), uana(x,y,x) \n");
  System.out.print("uniform grid on rectangle xmin,ymin,zmin to xmax,ymax,zmax\n");
  System.out.print(" \n");
  start = System.currentTimeMillis()/1000.0;
  ccxx = 1.0;
  ccyy = 2.0;
  cczz = 3.0;
  ccxy = 4.0;
  ccxz = 5.0;
  ccyz = 6.0;
  ccx  = 7.0;
  ccy  = 8.0;
  ccz  = 9.0;
  ccc  = 10.0;
  xmin = 0.0;
  xmax = 1.0;
  nx   = 5;
  ymin = 0.1;
  ymax = 1.1;
  ny   = 6;
  zmin = 0.2;
  zmax = 1.2;
  nz   = 7;
  System.out.print("ccxx="+ccxx+", ccyy="+ccyy+", cczz="+cczz+" \n");
  System.out.print("ccxy="+ccxy+", ccxz="+ccxz+", ccyz="+ccyz+" \n");
  System.out.print("ccx="+ccx+", ccy="+ccy+", ccz="+ccz+", ccc="+ccc+" \n");
  System.out.print("xmin="+xmin+", xmax="+xmax+", nx="+nx+" \n");
  System.out.print("ymin="+ymin+", ymax="+ymax+", ny="+ny+" \n");
  System.out.print("zmin="+zmin+", zmax="+zmax+", nz="+nz+" \n");

  neqn = (nx-2)*(ny-2)*(nz-2); // nx-2  internal, non boundary cells 
  cs = neqn;
  hx = (xmax-xmin)/(nx-1);
  System.out.print("hx="+hx+" \n"); 

  for(int i=0; i<nx; i++)
  {
    xg[i] = xmin+i*hx;
      System.out.print("xg["+i+"]="+xg[i]+" \n");
  }

  hy = (ymax-ymin)/(ny-1);
  System.out.print("hy="+hy+" \n"); 
  for(int j=0; j<ny; j++)
  {
    yg[j] = ymin+j*hy;
    System.out.print("yg("+j+")="+yg[j]+" \n");
  }

  hz = (zmax-zmin)/(nz-1);
  System.out.print("hz="+hz+" \n"); 
  for(int k=0; k<nz; k++)
  {
    zg[k] = zmin+k*hz;
    System.out.print("zg("+k+")="+zg[k]+" \n");
  }

  System.out.print("test uana(xmax,ymax,zmax)="+uana(xmax,ymax,zmax)+" \n");
  System.out.print("test f(xmax,ymax,zmax)="+f(xmax,ymax,zmax)+" \n");

  for(int i=0; i<neqn; i++)
  {
    for(int j=0; j<=neqn; j++)
	ut[i][j] = 0.0;
  }
  System.out.print("ut initialized \n");
  init_matrix(nx, ny, nz, ut, xg, yg, zg,
              ccxx, ccyy, cczz, ccxy, ccxz, ccyz, ccx, ccy, ccz, ccc);
  System.out.print("matrix initialized \n");
  if(debug) printB(nx, ny, nz, ut); // for debugging
  check_rows(nx, ny, nz, ut, xg, yg, zg); 

  System.out.print("running simeq \n");
  new simeq(neqn, ut, us);
  System.out.print("back from simeq \n");

  check_soln(nx, ny, nz, us, xg, yg, zg,
             ccxx, ccyy, cczz, ccxy, ccxz, ccyz, ccx, ccy, ccz, ccc);
  print(nx, ny, nz, xg, yg, zg, us);

  check(nx, ny, nz, xg, yg, zg, us);

  now = System.currentTimeMillis()/1000.0;
  System.out.println("total time "+(now-start)+" seconds");
  System.out.print("finished pde23js_eq.java \n");
} // end pde23js_eq

int S(int i, int j, int k, int ny, int nz) // zero based index includes boundary 
{                           // this is row of ut[] 
  return (i-1)*(ny-2)*(nz-2)+(j-1)*(nz-2)+(k-1); 
         // index into us and ut, RHS, no boundary 
} // end S 

void check(int nx, int ny, int nz, double xg[], double yg[], double zg[],
           double us[]) 
              // typically not known, instrumentation 
{
  double uexact, err, error, maxerror;
  
  error = 0.0;
  maxerror = 0.0;
  for(int i=1; i<nx-1; i++)
  {
    for(int j=1; j<ny-1; j++)
    { 
      for(int k=1; k<nz-1; k++)
      {
        uexact = uana(xg[i], yg[j], zg[k]);
        err = Math.abs(us[S(i,j,k,ny,nz)]-uexact);
        error = error + err;
        if(err>maxerror) maxerror=err;
      } //k
    } // j
  } // i
  System.out.print("check error="+error+", maxerror="+maxerror+" \n");
} // end check 

void check_rows(int nx, int ny, int nz, double ut[][], double xg[],
                double yg[], double zg[]) // only for(checking 
{	// row is in ut solution coordinates 
  int row, col;
  double sum, sum1;
  int cs = (nx-2)*(ny-2)*(nz-2);

  System.out.print("in check_rows \n");
  for(int i=1; i<nx-1; i++)
  {
    for(int j=1; j<ny-1; j++)
    {
      for(int k=1; k<nz-1; k++)
      {
        row = S(i,j,k,ny,nz);
        sum1 = 0.0;
        for(int ii=1; ii<nx-1; ii++)
        {
          for(int jj=1; jj<ny-1; jj++)
	  {
            for(int kk=1; kk<nz-1; kk++)
            {
		col = S(ii,jj,kk,ny,nz);
		sum1 = sum1 + ut[row][col]*uana(xg[ii],yg[jj],zg[kk]); // u is check solution 
            } // kk
          } // jj
        } // ii
        sum = sum1 - ut[row][cs];
        if(Math.abs(sum)>0.001)
        {
          System.out.print("row i="+i+", j="+j+", k="+k+" is bad err="+sum+" \n");
          System.out.print("ut[row][cs]="+ut[row][cs]+", sum1="+sum1+" \n"); 
        } // end if
      } // k
    } // j 
  } //i 
  System.out.print("finished check_rows \n");
} // end check_rows 

void print(int nx, int ny, int nz, double xg[], double yg[], double zg[],
           double us[]) // typically not known 
{
  double error = 0.0;
  double max_error = 0.0;
  double avg_error = 0.0;
  int neqn = (nx-2)*(ny-2)*(nz-2);

  System.out.print("exact solution u, computed us, error \n");
  for(int i=1; i<nx-1; i++)
  {
    for(int j=1; j<ny-1; j++)
    {
      for(int k=1; k<nz-1; k++)
      {
        error = uana(xg[i],yg[j],zg[k]) - us[S(i,j,k,ny,nz)];
        avg_error = avg_error + Math.abs(error);
        if(Math.abs(error)>max_error) max_error = Math.abs(error);

        System.out.print("xg["+i+"]="+xg[i]+", yg["+j+"]="+yg[j]+
                       ", zg["+k+"]="+zg[k]+" \n");
	System.out.print("u="+uana(xg[i],yg[j],zg[k])+", us="+us[S(i,j,k,ny,nz)]+
                       ", err="+error+" \n");
      } // k
    } // j
  } // i
  System.out.print("max_error="+max_error+", avg_error="+avg_error/neqn+" \n");
  System.out.print(" \n");
} // end print

void printB(int nx, int ny, int nz, double ut[][])
{
  System.out.print("matrix B includes RHS \n");
  int cs = (nx-2)*(ny-2)*(nz-2);
  for(int i=1; i<nx-1; i++)
  {
    for(int j=1; j<ny-1; j++)
    {
      for(int k=1; k<nz-1; k++)
      {
        for(int ii=1; ii<nx-1; ii++)
        {
          for(int jj=1; jj<ny-1; jj++)
          {
            for(int kk=1; kk<nz-1; kk++)
            {
              System.out.print("i="+i+", j="+j+", k="+k+", ii="+ii+
                             ", jj="+jj+", kk="+kk+", ut="+
	                     ut[S(i,j,k,ny,nz)][S(ii,jj,kk,ny,nz)]+" \n");
            } // kk
	  } // jj 
        } // ii 
        System.out.print("RHS=                 "+ut[S(i,j,k,ny,nz)][cs]+" \n");
      } // k
    } // j 
  } // i 
  System.out.print(" \n");
  System.out.print(" \n");
} //  end printB 

void init_matrix(int nx, int ny, int nz, double ut[][], double xg[],
                 double yg[], double zg[],
                 double ccxx, double ccyy, double cczz,
                 double ccxy, double ccxz, double ccyz,
                 double ccx, double ccy, double ccz, double ccc)
{
  double cxx[] = new double[30];  // nuderiv coefficients 
  double cyy[] = new double[30];
  double czz[] = new double[30];
  double cxy[] = new double[30];
  double cxz[] = new double[30];
  double cyz[] = new double[30];
  double cx[] = new double[30]; 
  double cy[] = new double[30];
  double cz[] = new double[30];
  int row;
  int cs = (nx-2)*(ny-2)*(nz-2);

  System.out.print("init_matrix running nx="+nx+", ny="+ny+", nz="+nz+" \n");
  System.out.print("cs="+cs+" \n");
  System.out.print("S(1,1,1,ny,nz)="+S(1,1,1,ny,nz)+" \n");
  System.out.print("S(nx-1,ny-1,nz-1,ny,nz)="+S(nx-1,ny-1,nz-1,ny,nz)+" \n");

  // cs is RHS, constant column subscript 
  // zero internal cells 
  for(int i=1; i<nx-1; i++)
  {
    for(int j=1; j<ny-1; j++)
    {
      for(int k=1; k<nz-1; k++)
      {
        row = S(i,j,k,ny,nz);
        // System.out.print("row S("+i+","+j+","+k+","+ny+","+nz+")=
        //                  "+S(i,j,k,ny,nz)+" \n");
        for(int ii=1; ii<nx-1; ii++)
        {
          for(int jj=1; jj<ny-1; jj++)
          {
            for(int kk=1; kk<nz-1; kk++)
            {
	      ut[row][S(ii,jj,kk,ny,nz)] = 0.0;
            } // kk
          } //jj
        } // ii
        ut[row][cs] = 0.0;
      } // k
    }
  }

  System.out.print("internal cells zeroed \n");
  for(int i=1; i<nx-1; i++) // inside boundary 
  {
    for(int j=1; j<ny-1; j++)
    {
      for(int k=1; k<nz-1; k++)
      {
        row = S(i,j,k,ny,nz); // equation 
        // ccxx*uxx+ccyy*uyy+cczz*uzz+ccxy*uxy+ccxz*uxz+ccyz*uyz+
        // ccx*ux+ccy*uy+ccz*uz+ccc*u=f(x,y,z)
        // for(each term 

        // ccxx*uxx  =  d^2 uana(x,y,z)/dx^2 
        new nuderiv(2,nx,i,xg,cxx);
        for(int ki=0; ki<nx; ki++)
        {
	  // if(i==1 && j==1) System.out.print("cxx["+ki+"]="+cxx[ki]+" \n");
          if(ki==0 || ki==nx-1)
          {
            ut[row][cs] = ut[row][cs] - ccxx*cxx[ki]*uana(xg[ki],yg[j],zg[k]);
          }	 
          else
          {
            ut[row][S(ki,j,k,ny,nz)] = ut[row][S(ki,j,k,ny,nz)] + ccxx*cxx[ki];
          }
        } // ki

        // ccyy*uyy  = d^2 uana(x,y,z)/dy^2 
        new nuderiv(2,ny,j,yg,cyy);
        for(int kj=0; kj<ny; kj++)
        {
	  // if(i==1 && j==1) System.out.print("cyy["+kj+"]="+cyy[kj]+" \n");
          if(kj==0 || kj==ny-1)
	  {
	    ut[row][cs] = ut[row][cs] - ccyy*cyy[kj]*uana(xg[i],yg[kj],zg[k]);
          }
          else
	  {
	    ut[row][S(i,kj,k,ny,nz)] = ut[row][S(i,kj,k,ny,nz)] + ccyy*cyy[kj];
	  }
        } // kj

        // cczz*uzz  = d^2 uana(x,y,z)/dz^2 
        new nuderiv(2,nz,k,zg,czz);
        for(int kk=0; kk<nz; kk++)
        {
	  // if(i==1 && j==1) System.out.print("cyy["+kj+"]="+cyy[kj]+" \n");
          if(kk==0 || kk==nz-1)
	  {
	    ut[row][cs] = ut[row][cs] - cczz*czz[kk]*uana(xg[i],yg[j],zg[kk]);
          }
          else
	  {
	    ut[row][S(i,j,kk,ny,nz)] = ut[row][S(i,j,kk,ny,nz)] + cczz*czz[kk];
	  }
        } // kj

        // ccxy*uxy  =  d^2 uana(x,y,z)/dxdy 
        new nuderiv(1,nx,i,xg,cx);
        new nuderiv(1,ny,j,yg,cy);
        for(int ki=0; ki<nx; ki++)
        {
          for(int kj=0; kj<ny; kj++)
          {
            if(ki==0 || ki==nx-1 || kj==0 || kj==ny-1)
            {
	      ut[row][cs] = ut[row][cs] - ccxy*cx[ki]*cy[kj]*
                                          uana(xg[ki],yg[kj],zg[k]);
            }
            else
	    {
	      ut[row][S(ki,kj,k,ny,nz)] = ut[row][S(ki,kj,k,ny,nz)] +
                                       ccxy*cx[ki]*cy[kj];
	    }
          } // kj 
        } // ki 

        // ccxz*uxz  =  d^2 uana(x,y,z)/dxdz 
        new nuderiv(1,nx,i,xg,cx);
        new nuderiv(1,nz,k,zg,cz);
        for(int ki=0; ki<nx; ki++)
        {
          for(int kk=0; kk<nz; kk++)
          {
            if(ki==0 || ki==nx-1 || kk==0 || kk==nz-1)
            {
	      ut[row][cs] = ut[row][cs] - ccxz*cx[ki]*cz[kk]*
                                          uana(xg[ki],yg[j],zg[kk]);
            }
            else
	    {
	      ut[row][S(ki,j,kk,ny,nz)] = ut[row][S(ki,j,kk,ny,nz)] +
                                       ccxz*cx[ki]*cz[kk];
	    }
          } // kk 
        } // ki 

        // ccyz*uyz  =  d^2 uana(x,y,z)/dydx 
        new nuderiv(1,ny,j,yg,cy);
        new nuderiv(1,nz,k,zg,cz);
        for(int kj=0; kj<ny; kj++)
        {
          for(int kk=0; kk<nz; kk++)
          {
            if(kj==0 || kj==ny-1 || kk==0 || kk==nz-1)
            {
	      ut[row][cs] = ut[row][cs] - ccyz*cy[kj]*cz[kk]*
                                          uana(xg[i],yg[kj],zg[kk]);
            }
            else
	    {
		ut[row][S(i,kj,kk,ny,nz)] = ut[row][S(i,kj,kk,ny,nz)] +
                                            ccyz*cy[kj]*cz[kk];
	    }
          } // kk
        } // kj 


        // ccx*ux  =  d uana(x,y,z)/dx 
        new nuderiv(1,nx,i,xg,cx);
        for(int ki=0; ki<nx; ki++)
        {
	  // if(i==1 && j==1) System.out.print("cx["+ki+"]="+cx[ki]+" \n");
          if(ki==0 || ki==nx-1)
          {
	    ut[row][cs] = ut[row][cs] - ccx*cx[ki]*uana(xg[ki],yg[j],zg[k]);
          }	 
          else
	  {
	    ut[row][S(ki,j,k,ny,nz)] = ut[row][S(ki,j,k,ny,nz)] + ccx*cx[ki];
	  }
        } // ki

        // ccy*uy  = d uana(x,y,z)/dy 
        new nuderiv(1,ny,j,yg,cy);
        for(int kj=0; kj<ny; kj++)
        {
	  // if(i==1 && j==1) System.out.print("cy["+kj+"]="+cy[kj]+" \n");
          if(kj==0 || kj==ny-1)
          {
	    ut[row][cs] = ut[row][cs] - ccy*cy[kj]*uana(xg[i],yg[kj],zg[k]);
          }
          else
          {
	    ut[row][S(i,kj,k,ny,nz)] = ut[row][S(i,kj,k,ny,nz)] + ccy*cy[kj];
          }
        } // kj

        // ccz*uz  = d uana(x,y,z)/dz 
        new nuderiv(1,nz,k,zg,cz);
        for(int kk=0; kk<nz; kk++)
        {
	  // if(i==1 && j==1) System.out.print("cz["+kk+"]="+cz[kk]+" \n");
          if(kk==0 || kk==nz-1)
          {
	    ut[row][cs] = ut[row][cs] - ccz*cz[kk]*uana(xg[i],yg[j],zg[kk]);
          }
          else
          {
	    ut[row][S(i,j,kk,ny,nz)] = ut[row][S(i,j,kk,ny,nz)] + ccz*cz[kk];
          }
        } // kk


        // ccc*u 
        ut[row][S(i,j,k,ny,nz)] = ut[row][S(i,j,k,ny,nz)] + ccc;


        // RHS constant 
        ut[row][cs] = ut[row][cs] + f(xg[i],yg[j],zg[k]);
        // end terms 
      } // k
    } // j 
  } // i 
} // end init_matrix 

void check_soln(int nx, int ny, int nz, double us[], double xg[],
                double yg[], double zg[],
                double ccxx, double ccyy, double cczz,
                double ccxy, double ccxz, double ccyz,
                double ccx, double ccy, double ccz, double ccc)
                  // check when solution unknown 
{
  double cxx[] = new double[30];  // nuderiv coefficients 
  double cyy[] = new double[30];
  double czz[] = new double[30];
  double cx[] = new double[30]; 
  double cy[] = new double[30];
  double cz[] = new double[30];
  int row;
  int cs = (nx-2)*(ny-2)*(nz-2);
  double smaxerr, err, uxx=0.0, uyy=0.0, uzz=0.0, uxy=0.0, uxz=0.0, uyz=0.0;
  double ux=0.0, uy=0.0, uz=0.0;

  // check ccxx*uxx(x,y,z)+ccyy*uyy(x,y,z)+cczz*uzz(x,y,z)+
  //       ccxy*uxy(x,y,z)+ccxy*uxy(x,y,z)+ccyz*uyz(x,y,z)+
  // ccx*ux(x,y,z)+ccy*uy(x,y,z)+ccz*uz(x,y,z)+ccc*uana(x,y,z)=f(x,y,z)
  // against computed solution us[S(i,j,k,ny,nz)] at each i,j,k 
  System.out.print("check computed solution against PDE \n");
  smaxerr = 0.0;

  for(int i=1; i<nx-1; i++)  // through dof equations 
  {
    new nuderiv(1,nx,i,xg,cx);
    new nuderiv(2,nx,i,xg,cxx);
    for(int j=1; j<ny-1; j++)
    {
      new nuderiv(1,ny,j,yg,cy);
      new nuderiv(2,ny,j,yg,cyy);
      for(int k=1; k<nz-1; k++)
      {
        new nuderiv(1,nz,k,zg,cz);
        new nuderiv(2,nz,k,zg,czz);

        ux  = 0.0;
        uxx = 0.0;
        for(int ki=0; ki<nx; ki++) // compute numerical derivative at (i,j) 
        {
          if(ki==0 || ki==nx-1)
          {
            ux  = ux  + cx [ki]*uana(xg[ki],yg[j],zg[k]); // boundary 
            uxx = uxx + cxx[ki]*uana(xg[ki],yg[j],zg[k]); // boundary 
          }
          else
          {
            ux  = ux  + cx [ki]*us[S(ki,j,k,ny,nz)];
            uxx = uxx + cxx[ki]*us[S(ki,j,k,ny,nz)];
          }
        } // ki

        uy  = 0.0;
        uyy = 0.0;
        for(int kj=0; kj<ny; kj++)
        {
          if(kj==0 || kj==ny-1)
          {
            uy  = uy  + cy [kj]*uana(xg[i],yg[kj],zg[k]); // boundary 
            uyy = uyy + cyy[kj]*uana(xg[i],yg[kj],zg[k]); // boundary 
          }
          else
          {
            uy  = uy  + cy [kj]*us[S(i,kj,k,ny,nz)];
            uyy = uyy + cyy[kj]*us[S(i,kj,k,ny,nz)];
          }
        } // kj

        uz  = 0.0;
        uzz = 0.0;
        for(int kk=0; kk<nz; kk++)
        {
          if(kk==0 || kk==nz-1)
          {
            uz  = uz  + cz [kk]*uana(xg[i],yg[j],zg[kk]); // boundary 
            uzz = uzz + czz[kk]*uana(xg[i],yg[j],zg[kk]); // boundary 
          }
          else
          {
            uz  = uz  + cz [kk]*us[S(i,j,kk,ny,nz)];
            uzz = uzz + czz[kk]*us[S(i,j,kk,ny,nz)];
          }
        } // kk


        uxy = 0.0;
        for(int ki=0; ki<nx; ki++) // compute numerical derivative at (i,j,k)
        {
          for(int kj=0; kj<ny; kj++)
          {
            if(ki==0 || ki==nx-1 || kj==0 || kj==ny-1)
            {
              uxy  = uxy  + cx[ki]*cy[kj]*uana(xg[ki],yg[kj],zg[k]); // bound
            }
            else
            {
	      uxy = uxy + cx[ki]*cy[kj]*us[S(ki,kj,k,ny,nz)];
            }
          } // kj
        } // ki

        uxz = 0.0;
        for(int ki=0; ki<nx; ki++) // compute numerical derivative at (i,j,k)
        {
          for(int kk=0; kk<nz; kk++)
          {
            if(ki==0 || ki==nx-1 || kk==0 || kk==nz-1)
            {
              uxz  = uxz  + cx[ki]*cz[kk]*uana(xg[ki],yg[j],zg[kk]); // bound
            }
            else
            {
              uxz = uxz + cx[ki]*cz[kk]*us[S(ki,j,kk,ny,nz)];
            }
          } // kk
        } // ki

        uyz = 0.0;
        for(int kj=0; kj<ny; kj++) // compute numerical derivative at (i,j,k) 
        {
          for(int kk=0; kk<nz; kk++)
          {
            if(kj==0 || kj==ny-1 || kk==0 || kk==nz-1)
            {
              uyz  = uyz  + cy[kj]*cz[kk]*uana(xg[i],yg[kj],zg[kk]); // bound
            }
            else
            {
              uyz = uyz + cy[kj]*cz[kk]*us[S(i,kj,kk,ny,nz)];
            }
          } // kk
        } // kj

        err = ccxx*uxx + ccyy*uyy + cczz*uzz + ccxy*uxy + ccxz*uxz + ccyz*uyz +
              ccx*ux + ccy*uy + ccz*uz + ccc*us[S(i,j,k,ny,nz)]
              - f(xg[i],yg[j],zg[k]); // minus right hand side 
        if(Math.abs(err) > smaxerr) smaxerr = Math.abs(err);
      } // k z
    } // j y 
  } // i x 

  System.out.print(" \n");
  System.out.print("check computed soln against PDE max error="+smaxerr+" \n");
  System.out.print(" \n");
  System.out.print("x="+xg[nx-2]+", y="+yg[ny-2]+", z="+zg[ny-2]+" \n");
  System.out.print("uxx="+uxx+", uyy="+uyy+", uzz="+uzz+" \n");
  System.out.print("ux="+ux+", uy="+uy+", uz="+uz+" \n");
  System.out.print("us="+us[S(nx-2,ny-2,nz-2,ny,nz)]+" \n");
  System.out.print("f(x,y,z)="+f(xg[nx-2],yg[ny-2],zg[nz-2])+" \n");
  System.out.print(" \n");
} // end Check_Soln 


double uana(double x,double y,double z) // ub(x,y,z) for boundary 
{
  return  x*x+y*y+z*z+x*y+x*z+y*z+1;
} // end uana

double f(double x,double y,double z) // RHS f(x,y,z)
{
 return  10*(x*x+y*y+z*z+x*y+x*z+y*z)+31*x+32*y+33*z+37;
} // end f

  public static void main (String[] args)
  {
    new pde23js_eq();
  }

} // end class pde23js_eq