// fem_check44_pla.java  Galerkin FEM
// solve uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) +
//     4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) +
//     7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) +
//     10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) =
//     F(x,y,z,t)
//     F(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +
//                  180 z^2*t + 200 y^2z*t + 44 t^3 +
//                  110 y^2z^2t + 66 xy + 12t^4 +
//                  24 z^4 + 36 y^4 + 48 x^4 +
//                  60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t
//
// boundary conditions computed using u(x,y,z,t)
// analytic solution is u(x,y,z,t) = t^4 + 2 z^4 + 3 y^4 + 4 x^4 +
//                                   5 y^2 z^2 t^2 + 6 x y t +
//                                   7 x z + 8 t + 9
// gauss-legendre integration used
// solve for Uijkl numerical approximation U(x_i,y_j,z_k,t_l)
// kg * ug = fg, solve simultaneous equations for U

// needs psimeq, laphi, gaulegf

import java.util.concurrent.*; // uses barrier.await() and others

class fem_check44_pla
{
  final int NP = 4; // number of processors
  CyclicBarrier barrieri1;
  CyclicBarrier barrieri2;

  final int nx = 5;
  final int ny = 5;
  final int nz = 5;
  final int nt = 5;
  final int nxyzt = nx*ny*nz*nt;
  double kg[] = new double[nxyzt*nxyzt]; // stifness matrix  
  double fg[] = new double[nxyzt];       // RHS
  double ug[] = new double[nxyzt];       // computed solution
  double Ua[] = new double[nxyzt];       // check solution
  double xg[] = new double[nx];
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double tg[] = new double[nt];
  laphi L = new laphi();
  final int npx = 5; // Gauss Legendre integration order
  final int npy = 5;
  final int npz = 5;
  final int npt = 5;
  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 tt[] = new double[npt+1];
  double wt[] = new double[npt+1]; // for Gauss-Legendre 

  class Integrate extends Thread
  {
    int myid;
    double val;

    Integrate(int id)
    {
      myid = id;
    }

    public void run()
    {
      try
      {
	System.out.println("Integrate "+myid+" at i1 await");
	barrieri1.await();
	System.out.println("Integrate "+myid+" after i1 await");
	int i = myid; // could be many
        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++)
                {
                  for(int iiii=1; iiii<nt-1; iiii++)
                  {
                    for(int jjjj=0; jjjj<nt; jjjj++)
                    {
                      // compute integral for k(i,j,k,l)  
                      val = 0.0;
                      for(int px=1; px<=npx; px++)
                      {
                        for(int py=1; py<=npy; py++)
                        {
                          for(int pz=1; pz<=npz; pz++)
                          {
                            for(int pt=1; pt<=npt; pt++)
                            {
		              val = val + wx[px]*wy[py]*wz[pz]*wt[pt]*
				          galk(xx[px],yy[py],zz[pz],tt[pt],
                                               i,j,ii,jj,iii,jjj,iiii,jjjj);
			    } // pt
		          } // pz
		        } // py
		      } // px
                      kg[ss(i,ii,iii,iiii,j,jj,jjj,jjjj)] = val;
		    } // end jjjj loop
		  } // end iiii loop
	        } // end jjj loop
              } // end iii loop
            } // end jj loop 
          } // end ii loop 
          System.out.println("Legendre integration="+val+", at i="+i+
       	                     ", j="+j+", ii="+(ny-1)+", jj="+ny);
          System.out.println("        , iii="+(nz-1)+", jjj="+nz+
                             ", iiii="+(nt-1)+", jjjj="+nt);
        } // end j loop
 
	// i from above
        // compute integral for f(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;
              for(int px=1; px<=npx; px++)
              {
                for(int py=1; py<=npy; py++)
                {
	          for(int pz=1; pz<=npz; pz++)
	          {
	            for(int pt=1; pt<=npt; pt++)
	            {
	              val = val + wx[px]*wy[py]*wz[pz]*wt[pt]*
			  galf(xx[px],yy[py],zz[pz],tt[pt],
                               i, ii, iii, iiii);
		    } // pt
	          } // pz
	        } // py
              } //px
              fg[s(i,ii,iii,iiii)] = val;
	    } // end iiii loop
          } // end iii loop
        } // end ii loop
        System.out.println("Legendre integration="+val+", f at i="+i+
                           ", ii="+(ny-1)+", iii="+(nz-1)+", iiii="+(nt-1));
        // end i, might be loop

	barrieri2.await();
	System.out.println("Integrate "+myid+" after i2 finished");
      }
      catch(InterruptedException e)
      {
	System.out.println("InterruptedException in master");
      }
      catch(BrokenBarrierException e)
      {
	System.out.println("BrokenBarrierException in master");
      }
    } // end Integrate run
  } // end Integrate

  fem_check44_pla()
  {
    double x, y, z, t, hx, hy, hz, ht;
    double val, err, avgerr, maxerr;
    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 t_start, t_init, t_kf, t_simeq, t_end;

    System.out.println("fem_check44_pla.c running");
    System.out.println("Given uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) +");
    System.out.println("    4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) +");
    System.out.println("    7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) +");
    System.out.println("    10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) =");
    System.out.println("  F(x,y,z,t)");
    System.out.println("  F(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +  ");
    System.out.println("               180 z^2*t + 200 y^2z*t + 44 t^3 +   ");
    System.out.println("               110 y^2z^2t + 66 xy + 12t^4 +       ");
    System.out.println("               24 z^4 + 36 y^4 + 48 x^4 +          ");
    System.out.println("               60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t");
    System.out.println("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries");
    System.out.println("Analytic solution u(x,y,z,t) =  t^4 + 2 z^4 + 3 y^4 +");
    System.out.println("     4 x^4 + 5 y^2 z^2 t^2 + 6 x y t + 7 x z + 8 t + 9");
    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);

    t_start = System.currentTimeMillis();
    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+i*hx;
      Ua[i] = uana(xg[i], ymin, zmin, tmin); // temp 
      System.out.println("i="+i+", Ua("+xg[i]+")="+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+ii*hy;
      Ua[ii] = uana(xmin, yg[ii], zmin, tmin); // temp 
      System.out.println("ii="+ii+", Ua("+yg[ii]+")="+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+iii*hz;
      Ua[iii] = uana(xmin, ymin, zg[iii], tmin); // temp 
      System.out.println("iii="+iii+", Ua("+zg[iii]+")="+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+iiii*ht;
      Ua[iiii] = uana(xmin, ymin, zmin, tg[iiii]); // temp 
      System.out.println("iiii="+iiii+", Ua("+tg[iiii]+")="+Ua[iiii]);
    } // iiii
    System.out.println(" ");

    // put solution in Ua vector 
    for(int i=0; i<nx; i++)
    {
      x = xmin+(double)i*hx;
      for(int ii=0; ii<ny; ii++)
      {
	y = ymin+(double)ii*hy;
        for(int iii=0; iii<nz; iii++)
        {
	  z = zmin+(double)iii*hz;
          for(int iiii=0; iiii<nt; iiii++)
          {
	    t = tmin+(double)iiii*ht;
            Ua[s(i,ii,iii,iiii)] = uana(x,y,z,t);
	  }
	}
      }
    }
    System.out.println(" ");

    // initialize k and f 
    for(int i=0; i<nx; i++)
    {
      for(int j=0; j<nx; j++)
      {
        for(int ii=0; ii<ny; ii++)
        {
          for(int jj=0; jj<ny; jj++)
          {
            for(int iii=0; iii<nz; iii++)
            {
              for(int jjj=0; jjj<nz; jjj++)
              {
                for(int iiii=0; iiii<nt; iiii++)
                {
                  for(int jjjj=0; jjjj<nt; jjjj++)
                  {
		    kg[ss(i,ii,iii,iiii,j,jj,jjj,jjjj)] = 0.0;
		  }
		}
	      }
	    }
	  }
        }
      }
    }

    for(int i=0; i<nx; i++)
    {
      for(int ii=0; ii<ny; ii++)
      {
        for(int iii=0; iii<nz; iii++)
        {
          for(int iiii=0; iiii<nt; iiii++)
          {
	    fg[s(i,ii,iii,iiii)] = 0.0;
	  }
	}
      }
    }

    // set boundary conditions, overkill faces 
    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];
            kg[ss(i,ii,iii,iiii,i,ii,iii,iiii)] = 1.0;
	    fg[s(i,ii,iii,iiii)] = uana(x,y,z,t);
	  } // iiii
        } // iii
      } // ii
    } // i

    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("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("calling gauleg zmin="+zmin+", zmax="+zmax+
                       ", npz="+npz);
    new gaulegf(zmin, zmax, zz, wz, npz); // zz and wz set up for integrations 
    System.out.println("calling gauleg tmin="+tmin+", tmax="+tmax+
                       ", npt="+npt);
    new gaulegf(tmin, tmax, tt, wt, npt); // tt and wt set up for integrations 

    // debug
    System.out.println("xx[1]="+xx[1]+", xx[2]="+xx[2]+
                       ", wx[1]="+wx[1]+", wx[2]="+wx[2]);
    System.out.println("yy[1]="+yy[1]+", yy[2]="+yy[2]+
                       ", wy[1]="+wy[1]+", wy[2]="+wy[2]);
    System.out.println("zz[1]="+zz[1]+", zz[2]="+zz[2]+
                       ", wz[1]="+wz[1]+", wz[2]="+wz[2]);
    System.out.println("tt[1]="+tt[1]+", tt[2]="+tt[2]+
                       ", wt[1]="+wt[1]+", wt[2]="+wt[2]);
    val = galk(xx[2],yy[2],zz[2],tt[2],1,1,1,1,1,1,1,1);
    System.out.println("galk(xx[2],yy[2],zz[2],tt[2])="+val);
    val = galf(xx[2],yy[2],zz[2],tt[2],1,1,1,1);
    System.out.println("galf(xx[2],yy[2],zz[2],tt[2])="+val);

    t_init = System.currentTimeMillis();
    System.out.println("initialization took "+((t_init-t_start)/1000.0)+" seconds");
    // compute entries in stiffness matrix 
    System.out.println("compute stiffness matrix  ");
    Integrate sthread[] = new Integrate[nx]; // allocate, overspace OK
    for(int i=1; i<nx-1; i++) sthread[i] = new Integrate(i); // construct
    try
    {
      barrieri1 = new CyclicBarrier(nx-2+1); // 
      barrieri2 = new CyclicBarrier(nx-2+1); // 
      for(int i=1; i<nx-1; i++) sthread[i].start(); // start 'run'
      System.out.println("master started threads ");
      barrieri1.await(); // let threads integrate
      System.out.println("master waiting for threads to integrate ");
      barrieri2.await(); // all threads finished integration
      System.out.println("master afyer all threads finished ");
    }
    catch(InterruptedException e)
    {
      System.out.println("InterruptedException in master");
    }
    catch(BrokenBarrierException e)
    {
      System.out.println("BrokenBarrierException in master");
    }

    System.out.println("k computed stiffness matrix");
    System.out.println("f computed forcing function");
    t_kf = System.currentTimeMillis();
    System.out.println("kf build took "+((t_kf-t_init)/1000.0)+" seconds");

    System.out.println("solve "+nxyzt+" equations in "+nxyzt+" unknowns");
    new psimeq(NP, nxyzt, kg, fg, ug);
    System.out.println("equations solved");
    t_simeq = System.currentTimeMillis();
    System.out.println("simeq took "+((t_simeq-t_kf)/1000.0)+" seconds");

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("ug computed Galerkin, Ua analytic, error  ");
    for(int i=0; i<nx; i++)
    {
      for(int ii=0; ii<ny; ii++)
      {
        for(int iii=0; iii<nz; iii++)
        {
          for(int iiii=0; iiii<nt; iiii++)
          {
	      err = Math.abs(ug[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+"]="+
                                ug[s(i,ii,iii,iiii)]+
			 ", Ua="+Ua[s(i,ii,iii,iiii)]+
			 ", err="+(ug[s(i,ii,iii,iiii)]-Ua[s(i,ii,iii,iiii)]));
	  }
	}
      }
    }
    System.out.println("                         maxerr="+maxerr+", avgerr="+
                       avgerr/(double)(nxyzt));
    t_end = System.currentTimeMillis();
    System.out.println("total took "+((t_end-t_start)/1000.0)+" seconds");

    System.out.println(" ");
  } // end fem_check44_pla constructor

  // utility functions to compute indices
  int s(int i, int ii, int iii, int iiii)
  {
    return i*ny*nz*nt+ii*nz*nt+iii*nt+iiii;
  } // end s

  int ss(int i, int ii, int iii, int iiii, int j, int jj, int jjj, int jjjj)
  {
    return (i*ny*nz*nt+ii*nz*nt+iii*nt+iiii)*nx*ny*nz*nt+
	   (j*ny*nz*nt+jj*nz*nt+jjj*nt+jjjj);
  } // end ss

  // PDE problem definition functions
  double F(double x, double y, double z, double t)
                                           // forcing function, F(x,y,z,t)
  {
    return 706.0 + 120.0*t*t + 140.0*y*y*z + 768.0*x + 180.0*z*z*t +
           200.0*y*y*z*t + 44.0*t*t*t + 110.0*y*y*z*z*t + 66.0*x*y +
           12.0*t*t*t*t + 24.0*z*z*z*z + 36.0*y*y*y*y + 48.0*x*x*x*x +
           60.0*y*y*z*z*t*t + 72.0*x*y*t + 84.0*x*z + 96.0*t;
  }

  double uana(double x, double y, double z, double t)
                                      // analytic solution and boundaries
  {
    return t*t*t*t + 2.0*z*z*z*z + 3.0*y*y*y*y + 4.0*x*x*x*x +
           5.0*y*y*z*z*t*t + 6.0*x*y*t + 7.0*x*z + 8.0*t + 9.0;
  }

  double galk(double x, double y, double z, double t,
              int i, int j, int ii, int jj, int iii, int jjj, int iiii, int jjjj)
                       // Galerkin k stiffness function for this problem
  {
    return
          (L.phipppp(x,j,nx-1,xg)*    L.phi(y,jj,ny-1,yg)*    L.phi(z,jjj,nz-1,zg)*
                                                              L.phi(t,jjjj,nt-1,tg)+
       2.0*    L.phi(x,j,nx-1,xg)*L.phipppp(y,jj,ny-1,yg)*    L.phi(z,jjj,nz-1,zg)*
                                                              L.phi(t,jjjj,nt-1,tg)+
       3.0*    L.phi(x,j,nx-1,xg)*    L.phi(y,jj,ny-1,yg)*L.phipppp(z,jjj,nz-1,zg)*
	                                                      L.phi(t,jjjj,nt-1,tg)+
       4.0*    L.phi(x,j,nx-1,xg)*    L.phi(y,jj,ny-1,yg)*    L.phi(z,jjj,nz-1,zg)*
	                                                  L.phipppp(t,jjjj,nt-1,tg)+
       5.0*   L.phip(x,j,nx-1,xg)*   L.phip(y,jj,ny-1,yg)*    L.phi(z,jjj,nz-1,zg)*
	                                                     L.phip(t,jjjj,nt-1,tg)+
       6.0*    L.phi(x,j,nx-1,xg)*  L.phipp(y,jj,ny-1,yg)*  L.phipp(z,jjj,nz-1,zg)*
	                                                      L.phi(t,jjjj,nt-1,tg)+
       7.0*    L.phi(x,j,nx-1,xg)*    L.phi(y,jj,ny-1,yg)*   L.phip(z,jjj,nz-1,zg)*
	                                                    L.phipp(t,jjjj,nt-1,tg)+
       8.0* L.phippp(x,j,nx-1,xg)*    L.phi(y,jj,ny-1,yg)*    L.phi(z,jjj,nz-1,zg)*
	                                                      L.phi(t,jjjj,nt-1,tg)+
       9.0*    L.phi(x,j,nx-1,xg)*  L.phipp(y,jj,ny-1,yg)*    L.phi(z,jjj,nz-1,zg)*
	                                                     L.phip(t,jjjj,nt-1,tg)+
      10.0*    L.phi(x,j,nx-1,xg)*    L.phi(y,jj,ny-1,yg)*   L.phip(z,jjj,nz-1,zg)*
	                                                     L.phip(t,jjjj,nt-1,tg)+
      11.0*    L.phi(x,j,nx-1,xg)*    L.phi(y,jj,ny-1,yg)*    L.phi(z,jjj,nz-1,zg)*
	                                                     L.phip(t,jjjj,nt-1,tg)+
      12.0*    L.phi(x,j,nx-1,xg)*    L.phi(y,jj,ny-1,yg)*    L.phi(z,jjj,nz-1,zg)*
	                                                      L.phi(t,jjjj,nt-1,tg))*
              (L.phi(x,i,nx-1,xg)*    L.phi(y,ii,ny-1,yg)*    L.phi(z,iii,nz-1,zg)*
                                                              L.phi(t,iiii,nt-1,tg));
  } // end galk used by integration

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

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