// pde23js_la  Galerkin FEM  
// solve ccxx*uxx+cyy*uyy+cczz*uzz+ccxy*uxy+ccxz*uxz+
//       ccyz*uyz+ccx*ux+ccy*uy+ccz*uzcc*u=F(x,y,z)
//       given f(x,y,z) =
//       Analytic solution u(x,y,z) = boundary and check
//       inside nx by ny by nz grid, Gauss Legendre integration

class pde23js_la
{
  laphi L = new laphi();
pde23js_la()
{
  double start, now;
  boolean debug = false;
  double ccxx = 1.0; // PDE definition
  double ccyy = 2.0;
  double cczz = 3.0;
  double ccxy = 4.0;
  double ccxz = 5.0;
  double ccyz = 6.0;
  double ccx  = 7.0;
  double ccy  = 8.0;
  double ccz  = 9.0;
  double cc   = 10.0;
  final int nx = 5;
  final int ny = 6;
  final int nz = 7;
  int nxyz = nx*ny*nz;
  double x, y, z, hx, hy, hz, tmp;
  int npx = 7; // Gauss Legendre integration order
  int npy = 8;
  int npz = 9;
  double xmin = 0.0; // problem parameters 
  double xmax = 1.0;
  double ymin = 0.1;
  double ymax = 1.1;
  double zmin = 0.2;
  double zmax = 1.2;
  double val, err, avgerr, maxerr;
  double xg[] = new double[nx]; 
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double xx[] = new double[npx+1]; // for Gauss-Legendre 
  double wx[] = new double[npx+1]; 
  double yy[] = new double[npy+1];
  double wy[] = new double[npy+1];
  double zz[] = new double[npz+1];
  double wz[] = new double[npz+1];
  double kg[][] = new double[nxyz][nxyz];  
  double fg[] = new double[nxyz];
  double ug[] = new double[nxyz];
  double Ua[] = new double[nxyz];

  System.out.print("pde23js_la.java running \n");
  System.out.print("solve ccxx*uxx+cyy*uyy+cczz*uzz+ccxy*uxy+ccxz*uxz+\n");
  System.out.print("ccyz*uyz+ccx*ux+ccy*uy+ccz*uz+cc*u=F(x,y,z)\n");
  System.out.print("given  F(x,y,z) \n");
  System.out.print("  xmin<=x<=xmax ymin<=y<=ymax Boundaries \n");
  System.out.print("  boundary and possibly check uana(x,y,z) \n");
  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+", cc="+cc+" \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");
  System.out.print("npx="+npx+", npy="+npy+", npz="+npz+" \n");
  start = System.currentTimeMillis()/1000.0;

  hx = (xmax-xmin)/(nx-1);
  for(int i=0; i<nx; i++)
  {
    xg[i] = xmin+i*hx;
  }  
  hy = (ymax-ymin)/(ny-1);
  for(int ii=0; ii<ny; ii++)
  {
    yg[ii] = ymin+ii*hy;
  }  
  hz = (zmax-zmin)/(nz-1);
  for(int iii=0; iii<ny; iii++)
  {
    zg[iii] = zmin+iii*hz;
  }  
  System.out.print("xg, yg, zg initialized \n");
  // 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];
        Ua[S(i,ii,iii,ny,nz)] = uana(x,y,z);
        if(debug)
          System.out.print("solution at i="+i+",x="+x+", ii="+ii+",y="+y+
		       ",iii="+iii+", z="+z+" is "+Ua[S(i,ii,iii,ny,nz)]+"\n");
      } // iii
    } // ii
  } // i
  System.out.print("Ua initialized \n");

  // initialize kg and fg 
  for(int i=0; i<nx; i++)
  {
    for(int ii=0; ii<ny; ii++)
    {
      for(int iii=0; iii<nz; iii++)
      {
        for(int j=0; j<nx; j++)
        {
          for(int jj=0; jj<ny; jj++)
          {
            for(int jjj=0; jjj<nz; jjj++)
            {
	      kg[S(i,ii,iii,ny,nz)][S(j,jj,jjj,ny,nz)] = 0.0;
	    } // jjj
	  } // jj
	} // j
      } // iii
    } // ii
  } // i
  for(int i=0; i<nx; i++)
  {
    for(int ii=0; ii<ny; ii++)
    {
      for(int iii=0; iii<nz; iii++)
      {
	  fg[S(i,ii,iii,ny,nz)] = 0.0;
      } // iii
    } // ii
  } // i
  System.out.print("kg and fg initialized \n");

  // set boundary conditions, six faces 
  for(int i=0; i<nx; i++) // all x,y at zmin, zmax
  {
    x = xg[i];
    for(int ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      for(int iii=0; iii<nz; iii=iii+nz-1) // iii=0, iii=nz-1
      {
        z = zg[iii];
        kg[S(i,ii,iii,ny,nz)][S(i,ii,iii,ny,nz)] = 1.0;
        fg[S(i,ii,iii,ny,nz)] = uana(x,y,z);
        System.out.print("boundary i="+i+",x="+x+", ii="+ii+",y="+y+
 		 ",iii="+iii+",z="+z+" is "+fg[S(i,ii,iii,ny,nz)]+"\n");
      } // iii
    } // ii
  } // i
  for(int i=0; i<nx; i++) // all x,z at ymin, ymax
  {
    x = xg[i];
    for(int iii=0; iii<nz; iii++)
    {
      z = zg[iii];
      for(int ii=0; ii<ny; ii=ii+ny-1) // ii=0, ii=ny-1
      {
        y = yg[ii];
        kg[S(i,ii,iii,ny,nz)][S(i,ii,iii,ny,nz)] = 1.0;
        fg[S(i,ii,iii,ny,nz)] = uana(x,y,z);
        System.out.print("boundary i="+i+",x="+x+", ii="+ii+",y="+y+
 		 ",iii="+iii+",z="+z+" is "+fg[S(i,ii,iii,ny,nz)]+"\n");
      } // ii
    } // iii
  } // i
  for(int ii=0; ii<ny; ii++) // all y,z at xmin, xmax
  {
    y = yg[ii];
    for(int iii=0; iii<nz; iii++)
    {
      z = zg[iii];
      for(int i=0; i<nx; i=i+nx-1) // i=0, i=nx-1
      {
        x = xg[i];
        kg[S(i,ii,iii,ny,nz)][S(i,ii,iii,ny,nz)] = 1.0;
        fg[S(i,ii,iii,ny,nz)] = uana(x,y,z);
        System.out.print("boundary i="+i+",x="+x+", ii="+ii+",y="+y+
 		 ",iii="+iii+",z="+z+" is "+fg[S(i,ii,iii,ny,nz)]+"\n");
      } // i
    } // iii
  } // ii
  System.out.print("kg and fg boundary set \n");

  // initialize xx, wx, yy, wy
  for(int k=0; k<=npx; k++)
  {
    xx[k] = 0.0;
    wx[k] = 0.0;
  }
  for(int k=0; k<=npy; k++)
  {
    yy[k] = 0.0;
    wy[k] = 0.0;
  }
  for(int k=0; k<=npz; k++)
  {
    zz[k] = 0.0;
    wz[k] = 0.0;
  }
  System.out.print("calling gauleg xmin="+xmin+", xmax="+xmax+
                     ", npx="+npx+"\n");
  new gaulegf(xmin, xmax, xx, wx, npx); // xx and wx set up for integrations 
  System.out.print("xx[1]="+xx[1]+", xx[2]="+xx[2]+
                     ", wx[1]="+wx[1]+", wx[2]="+wx[2]+"\n");
  System.out.print("calling gauleg ymin="+ymin+", ymax="+ymax+
                 ", npy="+npy+"\n");
  new gaulegf(ymin, ymax, yy, wy, npy); // yy and wy set up for integrations 
  System.out.print("yy[1]="+yy[1]+", yy[2]="+yy[2]+
                 ", wy[1]="+wy[1]+", wy[2]="+wy[2]+"\n");
  System.out.print("calling gauleg zmin="+zmin+", zmax="+zmax+
                 ", npz="+npz+"\n");
  new gaulegf(zmin, zmax, zz, wz, npz); // zz and wz set up for integrations 
  System.out.print("zz[1]="+zz[1]+", zz[2]="+zz[2]+
                 ", wz[1]="+wz[1]+", wz[2]="+wz[2]+"\n");
  System.out.print("calling F( xmax="+xmax+", ymax="+ymax+",zmax="+zmax+"\n");
  val = F(xmax, ymax, zmax);
  System.out.print("F(xmax,ymax,zmax)="+val+"\n");
  System.out.print("calling phi(xmin="+xmin+",j=0,nx-1="+(nx-1)+",xg\n");
  val = L.phi(xmin, 0, nx-1, xg);
  System.out.print("phi(xmin,j=0,nx-i,xg)="+val+"\n");
  System.out.print("calling phip(xmin="+xmin+",j=0,nx-1="+(nx-1)+",xg\n");
  val = L.phip(xmin, 0, nx-1, xg);
  System.out.print("phip(xmin,j=0,nx-i,xg)="+val+"\n");
  System.out.print("calling phipp(xmin="+xmin+",j=0,nx-1="+(nx-1)+",xg\n");
  val = L.phipp(xmin, 0, nx-1, xg);
  System.out.print("phipp(xmin,j=0,nx-i,xg)="+val+"\n");

  val = galk(xx[2],yy[2],zz[2],1,1,1,1,1,1,
             ccxx,ccyy,cczz,ccxy,ccxz,ccyz,ccx,ccy,ccz,cc,
             nx,ny,nz,xg,yg,zg);
  System.out.print("galk(xx[2],yy[2])="+val+"\n");
  val = galf(xx[2],yy[2],zz[2],1,1,1,nx,ny,nz,xg,yg,zg);
  System.out.print("galf(xx[2],yy[2])="+val+"\n");
  // compute entries in stiffness matrix 
  System.out.print("compute stiffness matrix  \n");
  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++)
            {
              // compute integral for k(i,ii,iii)  
              val = 0.0;
              for(int px=1; px<=npx; px++)
              {
                for(int py=1; py<=npy; py++)
                {
                  for(int pz=1; pz<=npz; pz++)
                  {
	            val = val + wx[px]*wy[py]*wz[pz]*
		        galk(xx[px],yy[py],zz[pz],i,ii,iii,j,jj,jjj,
                             ccxx,ccyy,cczz,ccxy,ccxz,ccyz,ccx,ccy,ccz,cc,
                             nx,ny,nz,xg,yg,zg);
		  } // end pz
                } // end py
              } // end px
              if(debug)
                System.out.print("Legendre integration="+val+", at i="+i+
                ", j="+j+", ii="+ii+", jj="+jj+", iii="+iii+", jjj="+jjj+"\n");
              kg[S(i,ii,iii,ny,nz)][S(j,jj,jjj,ny,nz)] = val;
	    } // end jjj
          } // end iii 
        } // end jj 
      } // end ii 
    } // end j 
  } // end 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++)
      {
        val = 0.0;
        for(int px=1; px<=npx; px++)
        {
          for(int py=1; py<=npy; py++)
          {
            for(int pz=1; pz<=npz; pz++)
            {
	      val = val + wx[px]*wy[py]*wz[pz]*
		    galf(xx[px],yy[py],zz[pz],i,ii,iii,nx,ny,nz,xg,yg,zg);
            } // end pz
          } // end py
        } // end px
        if(debug)
          System.out.print("Legendre integration="+val+", f at i="+i+
                         ", ii="+ii+"\n");
        fg[S(i,ii,iii,ny,nz)] = val;
      } // end iii
    } // end ii
  } // end i
  System.out.print("k computed stiffness matrix \n");
  System.out.print("f computed forcing function \n");
  new simeq(kg, fg, ug);

  avgerr = 0.0;
  maxerr = 0.0;
  System.out.print("ug computed Galerkin, Ua analytic, error \n");
  for(int i=0; i<nx; i++)
  {
    for(int ii=0; ii<ny; ii++)
    {
      for(int iii=0; iii<nz; iii++)
      {
        err = Math.abs(ug[S(i,ii,iii,ny,nz)]-Ua[S(i,ii,iii,ny,nz)]);
        if(err>maxerr) maxerr = err;
        avgerr = avgerr + err;
        System.out.print("ug["+i+","+ii+","+iii+"]="+ug[S(i,ii,iii,ny,nz)]+
                         ", Ua="+Ua[S(i,ii,iii,ny,nz)]+", err="+
                         (ug[S(i,ii,iii,ny,nz)]-Ua[S(i,ii,iii,ny,nz)])+"\n");
      } // iii
    } // ii
  } // i
  System.out.print("                  maxerr="+maxerr+", avgerr="+
                 (avgerr/(nxyz))+"\n");
  now = System.currentTimeMillis()/1000.0;
  System.out.println("total time "+(now-start)+" seconds");
  System.out.print("pde23js_la.java finished\n");
} // end pde23js_la

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

// from uana23_fem.js  problem definition
double F(double x, double y, double z) // RHS  
{
  return 10.0*x*x + 10.0*x*y + 10.0*x*z + 10.0*y*y + 10.0*y*z + 10.0*z*z+
	 31.0*x + 32.0*y + 33.0*z + 37.0;
} // end f

//       Analytic solution
double uana(double x, double y, double z)
{
  return x*x+y*y+z*z+x*y+x*z+y*z+1.0;
} // end u

// needs laphi.java available
// for finite element method (not needed if coefficients available)
double galk(double x, double y, double z, int i, int ii, int iii,
            int j, int jj, int jjj,
            double ccxx, double ccyy, double cczz, double ccxy,
            double ccxz, double ccyz, double ccx,  double ccy,
            double ccz,  double cc,
            int nx, int ny, int nz, double xg[], double yg[], double zg[])
{           // Galerkin k stiffness function for this problem 
  return (ccxx * L.phipp(x,j,nx-1,xg)*  L.phi(y,jj,ny-1,yg)*  L.phi(z,jjj,nz-1,zg) +
          ccyy *   L.phi(x,j,nx-1,xg)*L.phipp(y,jj,ny-1,yg)*  L.phi(z,jjj,nz-1,zg) +
          cczz *   L.phi(x,j,nx-1,xg)*  L.phi(y,jj,ny-1,yg)*L.phipp(z,jjj,nz-1,zg) +
          ccxy *  L.phip(x,j,nx-1,xg)* L.phip(y,jj,ny-1,yg)*  L.phi(z,jjj,nz-1,zg) +
          ccxz *  L.phip(x,j,nx-1,xg)*  L.phi(y,jj,ny-1,yg)* L.phip(z,jjj,nz-1,zg) +
          ccyz *   L.phi(x,j,nx-1,xg)* L.phip(y,jj,ny-1,yg)* L.phip(z,jjj,nz-1,zg) +
          ccx  *  L.phip(x,j,nx-1,xg)*  L.phi(y,jj,ny-1,yg)*  L.phi(z,jjj,nz-1,zg) +
          ccy  *   L.phi(x,j,nx-1,xg)* L.phip(y,jj,ny-1,yg)*  L.phi(z,jjj,nz-1,zg) +
          ccz  *   L.phi(x,j,nx-1,xg)*  L.phi(y,jj,ny-1,yg)* L.phip(z,jjj,nz-1,zg) +
          cc   *   L.phi(x,j,nx-1,xg)*  L.phi(y,jj,ny-1,yg)*  L.phi(z,jjj,nz-1,zg) )*
                   L.phi(x,i,nx-1,xg)*  L.phi(y,ii,ny-1,yg)*  L.phi(z,iii,nz-1,zg);
} // end galk used by integration 

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

public static void main (String[] args) // "main" required
{
  new pde23js_la(); // construct and execute
} // end main

} // end class pde23js_la