// fem_bihar3d_la.java  Galerkin FEM biharmonic 4th order PDE in three dimensions
//
// solve  uxxxx(x,y,z) + uyyyy(x,y,z) + uzzzz(x,y,z) + 2*uxxyy(x,y,z) + 
//        2*uxxzz(x,y,z) + 2*uyyzz(x,y,z) = f(x,y,z)
//
//        f(x,y,z) := 0.24  ,
//
//        in the cube   xmin<x<xmax and ymin<y<ymax and zmin<z<xmax  ,
//
//  Dirichlet boundary conditions are computed using:
//  ub(x,y,z) := 0.01*x*x*x*x + 0.02*y*y*y*y + 0.03*z*z*z*z - 0.04*x*x*y*y -
//               0.05*y*y*z*z - 0.06*x*x*z*z + 0.1*x*x*x*y + 
//               0.3*y*y*y*z + 0.2*z*z*z*x - x*y - y*z - x*z + 1
 
import java.io.*;
//                uses laphi.java and simeq.java

class fem_bihar3d_la
{
  // i, j, k, nx, ii, jj, kk, ny, iii, jjj, kkk, xg, yg, zg
  //                needed by galk or galf, available at call 
    int i, j, k, ii, jj, kk, iii, jjj, kkk, nx, ny, nz;
  double xg[]; // X grig can esily be made non uniform
  double yg[]; // Y grid can esily be made non uniform
  double zg[]; // Y grid can esily be made non uniform
  laphi L = new laphi();
    double hx, hy, hz;
  double ug[]; // solution, computer indexing
  double maxerr;

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

  void run(double xmax, double ymax, double zmax, 
           int anx, int any, int anz, int npx, int npy, int npz)
  {
    nx = anx;
    ny = any;
    nz = anz;
    int nxyz = nx*ny*nz;
    xg = new double[nx];
    yg = new double[ny];
    zg = new double[nz];
    double x, y, z;
    // npx, npy, npz  Gauss Legendre integration order
    double xx[] = new double[npx+1];
    double wx[] = new double[npx+1]; // for Gauss-Legendre 
    double yy[] = new double[npy+1];
    double wy[] = new double[npy+1]; // for Gauss-Legendre 
    double zz[] = new double[npz+1];
    double wz[] = new double[npz+1]; // for Gauss-Legendre 
    double val, err, avgerr;
    int px, py, pz;

    double xmin = 0.0; // problem parameters 
    double ymin = 0.0;
    double zmin = 0.0;

    double kg[][] = new double[nxyz][nxyz]; // stifness matrix  
    double fg[] = new double[nxyz];  // RHS
    ug = new double[nxyz]; // computed solution
    double Ua[] = new double[nxyz]; // analytic solution, for testing

    System.out.println("fem_bihar3d_la.java running.");
    System.out.println("solve  uxxxx(x,y,z) + uyyyy(x,y,z) + uzzzz(x,y,z) + 2*uxxyy(x,y,z) +");
    System.out.println("       2*uxxzz(x,y,z) + 2*uyyzz(x,y,z) = f(x,y,z)");
    System.out.println(" ");
    System.out.println("       f(x,y) := 0.24");
    System.out.println(" ");
    System.out.println("       in the cube   xmin<x<xmax and ymin<y<ymax and zmin<z<xmax  ,");
    System.out.println(" ");
    System.out.println("Dirichlet boundary conditions are computed using:");
    System.out.println("ub(x,y,z) := 0.01*x*x*x*x + 0.02*y*y*y*y + 0.03*z*z*z*z - 0.04*x*x*y*y -");
    System.out.println("             0.05*y*y*z*z - 0.06*x*x*z*z + 0.1*x*x*x*y + ");
    System.out.println("             0.3*y*y*y*z + 0.2*z*z*z*x - x*y - y*z - x*z + 1");
    System.out.println(" ");
    System.out.println("xmin="+xmin+", xmax="+xmax+
                       ", ymin="+ymin+", ymax="+ymax+
                       ", zmin="+zmin+", zmax="+zmax);
    System.out.println("nx="+nx+", ny="+ny+", nz="+nz);

    System.out.println("x grid and analytic solution at ymin, zmin  ");
    hx = (xmax-xmin)/(double)(nx-1);
    for(i=0; i<nx; i++)
    {
      xg[i] = xmin+i*hx;
      System.out.println("i="+i+", Ua("+xg[i]+")="+ub(xg[i],ymin,zmin));
    }  
    System.out.println(" ");
    System.out.println("y grid and analytic solution at xmin, zmin ");
    hy = (ymax-ymin)/(double)(ny-1);
    for(ii=0; ii<ny; ii++)
    {
      yg[ii] = ymin+ii*hy;
      System.out.println("ii="+ii+", Ua("+yg[ii]+")="+ub(xmin, yg[ii], zmin));
    }  
    System.out.println(" ");
    System.out.println("z grid and analytic solution at xmin,ymin ");
    hz = (zmax-zmin)/(double)(nz-1);
    for(iii=0; iii<nz; iii++)
    {
      zg[iii] = zmin+iii*hz;
      System.out.println("iii="+iii+", Ua("+zg[iii]+")="+ub(xmin, ymin, zg[iii]));
    }  
    System.out.println(" ");
    // put solution in Ua vector 
    for(i=0; i<nx; i++)
    {
      x = xmin+(double)i*hx;
      for(ii=0; ii<ny; ii++)
      {
	y = ymin+(double)ii*hy;
        for(iii=0; iii<nz; iii++)
        {
	  z = zmin+(double)iii*hz;
          Ua[s(i,ii,iii)] = ub(x,y,z);
          // System.out.println("solution at i="+i+",x="+x+", ii="+ii+",y="+y+
          //                    ",z="+z+" is"+Ua[s(i,ii,iii)]);
	}
      }
    }
    System.out.println(" ");

    // initialize k and f 
    for(i=0; i<nx; i++)
    {
      for(j=0; j<nx; j++)
      {
        for(ii=0; ii<ny; ii++)
        {
          for(jj=0; jj<ny; jj++)
          {
            for(iii=0; iii<nz; iii++)
            {
              for(jjj=0; jjj<nz; jjj++)
              {
		kg[s(i,ii,iii)][(s(j,jj,jjj))] = 0.0;
	      }
	    }
	  }
        }
      }
    }
    for(i=0; i<nx; i++)
    {
      for(ii=0; ii<ny; ii++)
      {
        for(iii=0; iii<nz; iii++)
        {
	  fg[s(i,ii,iii)] = 0.0;
	}
      }
    }
    // set boundary conditions, six sides 
    for(i=0; i<nx; i++) // x, y, zmin
    {
      x = xg[i];
      for(ii=0; ii<ny; ii++)
      {
        y = yg[ii];
        iii = 0;
        z = zg[iii];
        kg[s(i,ii,iii)][s(i,ii,iii)] = 1.0;
        fg[s(i,ii,iii)] = ub(x,y,z);
      }
    }
    for(i=0; i<nx; i++) // x, y, zmax
    {
      x = xg[i];
      for(ii=0; ii<ny; ii++)
      {
        y = yg[ii];
        iii = nz-1;
        z = zg[iii];
        kg[s(i,ii,iii)][s(i,ii,iii)] = 1.0;
        fg[s(i,ii,iii)] = ub(x,y,z);
      }
    }
    i = 0;
    x = xg[i];
    for(ii=0; ii<ny; ii++) // xmin, y, z
    {
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        kg[s(i,ii,iii)][s(i,ii,iii)] = 1.0;
        fg[s(i,ii,iii)] = ub(x,y,z);
      }
    }
    i = nx-1;
    x = xg[i];
    for(ii=0; ii<ny; ii++) // xmax, y, z
    {
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        kg[s(i,ii,iii)][s(i,ii,iii)] = 1.0;
        fg[s(i,ii,iii)] = ub(x,y,z);
      }
    }
    for(i=0; i<nx; i++) // x, ymin, z
    {
      x = xg[i];
      ii = 0;
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        kg[s(i,ii,iii)][s(i,ii,iii)] = 1.0;
        fg[s(i,ii,iii)] = ub(x,y,z);
      }
    }
    for(i=0; i<nx; i++) // x, ymax, z
    {
      x = xg[i];
      ii = ny-1;
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        kg[s(i,ii,iii)][s(i,ii,iii)] = 1.0;
        fg[s(i,ii,iii)] = ub(x,y,z);
      }
    }

    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("xx[1]="+xx[1]+", xx[2]="+xx[2]+
                       ", wx[1]="+wx[1]+", wx[2]="+wx[2]);
    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("yy[1]="+yy[1]+", yy[2]="+yy[2]+
                       ", wy[1]="+wy[1]+", wy[2]="+wy[2]);
    new gaulegf(zmin, zmax, zz, wz, npz); // zz and wz set up for integrations 
    System.out.println("zz[1]="+zz[1]+", zz[2]="+zz[2]+
                       ", wz[1]="+wz[1]+", wz[2]="+wz[2]);
    i=1; ii=1; iii=1; j=1; jj=1; jjj=1; 
    val = galk(xx[1],yy[1],zz[1]);
    System.out.println("galk(xx[1],yy[1],zz[1])="+val);
    val = galf(xx[2],yy[2],zz[2]);
    System.out.println("galf(xx[2],yy[2],zz[2])="+val);


    // compute entries in stiffness matrix 
    System.out.println("compute stiffness matrix  ");
    for(i=1; i<nx-1; i++)
    {
      for(j=0; j<nx; j++)
      {
        for(ii=1; ii<ny-1; ii++)
        {
          for(jj=0; jj<ny; jj++)
          {
            for(iii=1; iii<nz-1; iii++)
            {
              for(jjj=0; jjj<nz; jjj++)
              {
                // compute integral for k(i,j,k)  
                val = 0.0;
                for(px=1; px<=npx; px++)
                {
                  for(py=1; py<=npy; py++)
                  {
                    for(pz=1; pz<=npz; pz++)
                    {
			val = val + wx[px]*wy[py]*wz[pz]*galk(xx[px],yy[py],zz[pz]);
		    }
	          }
                }
                // System.out.println("Legendre integration="+val+", at i="+i+
                //                    ", j="+j+", ii="+ii+", jj="+jj);
                kg[s(i,ii,iii)][s(j,jj,jjj)] = val;
	      } // end jjj loop
	    } // end iii loop
          } // end jj loop 
        } // end ii loop 
      } // end j loop 
    } // end i loop 

    // compute integral for f(i,j,k)  
    for(i=1; i<nx-1; i++)
    {
      for(ii=1; ii<ny-1; ii++)
      {
        for(iii=1; iii<nz-1; iii++)
        {
          val = 0.0;
          for(px=1; px<=npx; px++)
          {
            for(py=1; py<=npy; py++)
            {
              for(pz=1; pz<=npz; pz++)
              {
		  val = val + wx[px]*wy[py]*wz[pz]*galf(xx[px],yy[py],zz[pz]);
	      }
            }
          }
          // System.out.println("Legendre integration="+val+", f at i="+i+
          //                    ", ii="+ii);
          fg[s(i,ii,iii)] = val;
	}
      }
    }
    System.out.println("k computed stiffness matrix, see above  ");
    System.out.println("f computed forcing function, see above  ");
    new simeq(kg, fg, ug);

    write_soln(); // to file x y soln (move z later)

    // check_soln(xmax, ymax, zmax, nx, ny, nz); // independent of nowing analytic solution

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("ug computed Galerkin, Ua analytic, error  ");
    for(i=0; i<nx; i++)
    {
      for(ii=0; ii<ny; ii++)
      {
        for(iii=0; iii<nz; iii++)
        {
	  err = Math.abs(ug[s(i,ii,iii)]-Ua[s(i,ii,iii)]);
          if(err>maxerr) maxerr = err;
          avgerr = avgerr + err;
          System.out.println("ug["+i+","+ii+","+iii+"]="+ug[s(i,ii,iii)]+", Ua="+Ua[s(i,ii,iii)]+
                             ", err="+(ug[s(i,ii,iii)]-Ua[s(i,ii,iii)]));
	}
      }
    }
    System.out.println("xmax="+xmax+", ymax="+ymax+", zmax="+zmax);
    System.out.println("nx="+nx+", ny="+ny+", nz="+nz+", npx="+npx+", npy="+npy+", npz="+npz);
    System.out.println("                  maxerr="+maxerr+", avgerr="+
                       avgerr/(double)(nx*ny*nz));
    System.out.println("\n");
  } // end fem_bihar3d_la constructor

  // PDE problem definition functions  f and ub

  double F(double x, double y, double z)
  {
    return 0.24; // kept general for future use
  }

  double ub(double x, double y, double z)
  {
    return 0.01*x*x*x*x + 0.02*y*y*y*y + 0.03*z*z*z*z - 0.04*x*x*y*y -
	   0.05*y*y*z*z - 0.06*x*x*z*z + 0.1*x*x*x*y + 0.2*z*z*z*x -
           0.03*y*y*y*z - x*y - y*z - x*z + 1.0;
  }

  double u(double x, double y, double z)
  {
    return ub(x,y,z);
  }

  double uxxxx(double x, double y, double z)
  {
    return 0.24;
  }

  double uyyyy(double x, double y, double z)
  {
    return 0.48;
  }

  double uzzzz(double x, double y, double z)
  {
    return 0.72;
  }

  double uxxyy(double x, double y, double z)
  {
    return -0.16;
  }

  double uxxzz(double x, double y, double z)
  {
    return -0.24;
  }

  double uyyzz(double x, double y, double z)
  {
    return -0.20;
  }

  double galk(double x, double y, double z)
  {                        // Galerkin k stiffness function for this problem 
    // solve  uxxxx(x,y,z) + uyyyy(x,y,z) + uzzzz(x,y,z) 
    //        2*uxxyy(x,y,z) + 2*uxxzz(x,y,z) + 2*uyyzz(x,y,z) = f(x,y,z)
    return ( L.phipppp(x,j,nx-1,xg)*     L.phi(y,jj,ny-1,yg)*     L.phi(z,jjj,nz-1,zg) +
                 L.phi(x,j,nx-1,xg)* L.phipppp(y,jj,ny-1,yg)*     L.phi(z,jjj,nz-1,zg) +
                 L.phi(x,j,nx-1,xg)*     L.phi(y,jj,ny-1,yg)* L.phipppp(z,jjj,nz-1,zg) +
        2.0*   L.phipp(x,j,nx-1,xg)*   L.phipp(y,jj,ny-1,yg)*     L.phi(z,jjj,nz-1,zg) +
        2.0*   L.phipp(x,j,nx-1,xg)*     L.phi(y,jj,ny-1,yg)*   L.phipp(z,jjj,nz-1,zg) +
        2.0*     L.phi(x,j,nx-1,xg)*   L.phipp(y,jj,ny-1,yg)*   L.phipp(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) // 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 

  void write_soln()
  {
    try
    {
      FileWriter fout = new FileWriter("fem_bihar3d_la.dat");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing fem_bihar3d_la.dat");

      for(int i=0; i<nx; i++)
      {
        for(int ii=0; ii<ny; ii++)
        {
	  iii = nz-1;
          fileout.write(" "+xg[i]+" "+yg[ii]+" "+ug[s(i,ii,iii)]+" \n");
        } // end y
        fileout.write(" \n");
      } // end x
      fileout.close();
      System.out.println("finished writing fem_bihar3d_la.dat");
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("fem_bihar3d_la.dat  not opened");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("fem_bihar3d_la.java output-file error");
    } // end try
  } // end write_soln


  void check_soln(double xmax, double ymax, int nx, int ny)
                  // for all interior points
  {

  } // 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

  int s(int i, int j, int k) // subscript for kg and Ua
  {
    return i*ny*nz+ii*nz+iii;
  }

  public static void main (String[] args)
  {
    fem_bihar3d_la F2 = new fem_bihar3d_la();
  
    F2.run(2.0, 2.0, 2.0, 6, 6, 6, 6, 6 ,6); // xmax, ymax, zmax, nx, ny, nz, npx, npy, npz

    try
    {
      Runtime runtime = Runtime.getRuntime();
      Process proc = runtime.exec("gnuplot -persist fem_bihar3d_la.plot");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("fem_bihar3d_la.java runtime error");
    } // end try
  }
} // end class fem_bihar3d_la