// pde24js_la  Galerkin FEM  
// solve ccxx*uxx+ccyy*uyy+cczz*uzz+cctt*utt+ccxy*uxy+ccxz*uxz+
//       ccxt*uxt+ccyz*uyz+ccyt*uyt+cczt*uzt+
//       ccx*ux+ccy*uy+ccz*uz+cct*ut+cc*u=F(x,y,z,t)
//       given f(x,y,z,t) =
//       Analytic solution u(x,y,z,t) = boundary and check
//       inside nx by ny by nz by nt grid, Gauss Legendre integration

class pde24js_la
{
  laphi L = new laphi();
pde24js_la()
{
  double start, prev, now;
  boolean debug = false;
  double ccxx = 1.0; // PDE definition
  double ccyy = 2.0;
  double cczz = 3.0;
  double cctt = 4.0;
  double ccxy = 5.0;
  double ccxz = 6.0;
  double ccxt = 7.0;
  double ccyz = 8.0;
  double ccyt = 9.0;
  double cczt = 10.0;
  double ccx  = 11.0;
  double ccy  = 12.0;
  double ccz  = 13.0;
  double cct  = 14.0;
  double cc   = 15.0;
  final int nx = 5;
  final int ny = 5;
  final int nz = 5;
  final int nt = 5;
  int nxyzt = nx*ny*nz*nt;
  double x, y, z, t, hx, hy, hz, ht, tmp;
  int npx = 7; // Gauss Legendre integration order
  int npy = 7;
  int npz = 7;
  int npt = 7;
  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 tmin = 0.3;
  double tmax = 1.3;
  double val, err, avgerr, maxerr;
  double xg[] = new double[nx]; 
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double tg[] = new double[nt];
  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 tt[] = new double[npt+1];
  double wt[] = new double[npt+1];
  double kg[][] = new double[nxyzt][nxyzt];  
  double fg[] = new double[nxyzt];
  double ug[] = new double[nxyzt];
  double Ua[] = new double[nxyzt];

  System.out.print("pde24js_la.js running \n");
  System.out.print("solve ccxx*uxx+ccyy*uyy+cczz*uzz+cctt*utt+ccxy*uxy+\n");
  System.out.print("ccxz*uxz+ccxt*uxt+ccyz*uyz+ccyt*uyt+cczt*uzt+ \n");
  System.out.print("ccx*ux+ccy*uy+ccz*uz+cct*ut+cc*u=F(x,y,z,t)\n");
  System.out.print("given  F(x,y,z,t) \n");
  System.out.print("  xmin<=x<=xmax ymin<=y<=ymax Boundaries \n");
  System.out.print("  zmin<=z<=zmax tmin<=t<=tmax Boundaries \n");
  System.out.print("  boundary and possibly check uana(x,y,z,t) \n");
  System.out.print("ccxx="+ccxx+", ccyy="+ccyy+", cczz="+cczz+" \n");
  System.out.print("cctt="+cctt+", ccxy="+ccxy+", ccxz="+ccxz+" \n");
  System.out.print("ccxt="+ccxt+", ccyz="+ccyz+", ccyt="+ccyt+" \n");
  System.out.print("cczt="+cczt+", ccx=" +ccx+ ", ccy=" +ccy+ " \n");
  System.out.print("ccz="+ccz+  ", cct=" +cct+ ", 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("tmin="+tmin+", tmax="+tmax+", nt="+nt+" \n");
  now = System.currentTimeMillis()/1000.0;
  prev = now;
  start = now;

  hx = (xmax-xmin)/(nx-1);
  for(int i=0; i<nx; i++)
  {
    xg[i] = xmin+i*hx;
  }  
  hy = (ymax-ymin)/(ny-1);
  for(int i2=0; i2<ny; i2++)
  {
    yg[i2] = ymin+i2*hy;
  }  
  hz = (zmax-zmin)/(nz-1);
  for(int i3=0; i3<ny; i3++)
  {
    zg[i3] = zmin+i3*hz;
  }  
  ht = (tmax-tmin)/(nt-1);
  for(int i4=0; i4<nt; i4++)
  {
    tg[i4] = tmin+i4*ht;
  }  
  System.out.print("xg, yg, zg, tg initialized \n");
  // put solution in Ua vector 
  for(int i=0; i<nx; i++)
  {
    x = xg[i];
    for(int i2=0; i2<ny; i2++)
    {
      y = yg[i2];
      for(int i3=0; i3<nz; i3++)
      {
        z = zg[i3];
        for(int i4=0; i4<nt; i4++)
        {
          t = tg[i4];
          Ua[S(i,i2,i3,i4,ny,nz,nt)] = uana(x,y,z,t);
          if(debug)
            System.out.print("solution at i="+i+",x="+x+", i2="+i2+
                     ",y="+y+",i3="+i3+", z="+z+",i4="+i4+", t="+t+
                     " is "+Ua[S(i,i2,i3,i4,ny,nz,nt)]+"\n");
	} // i4
      } // i3
    } // i2
  } // i
  System.out.print("Ua initialized \n");

  // initialize kg and fg 
  for(int i=0; i<nx; i++)
  {
    for(int i2=0; i2<ny; i2++)
    {
      for(int i3=0; i3<nz; i3++)
      {
        for(int i4=0; i4<nt; i4++)
        {
          for(int j=0; j<nx; j++)
          {
            for(int j2=0; j2<ny; j2++)
            {
              for(int j3=0; j3<nz; j3++)
              {
                for(int j4=0; j4<nt; j4++)
                {
	          kg[S(i,i2,i3,i4,ny,nz,nt)][S(j,j2,j3,j4,ny,nz,nt)] = 0.0;
		} // j4
	      } // j3
	    } // j2
	  } // j
	} // i4
      } // i3
    } // i2
  } // i
  for(int i=0; i<nx; i++)
  {
    for(int i2=0; i2<ny; i2++)
    {
      for(int i3=0; i3<nz; i3++)
      {
        for(int i4=0; i4<nt; i4++)
        {
	  fg[S(i,i2,i3,i4,ny,nz,nt)] = 0.0;
	} // i4
      } // i3
    } // i2
  } // i
  System.out.print("kg and fg initialized \n");

  // set boundary conditions, all faces, any value at max 
  for(int i=0; i<nx; i++) // all x,y,z at tmin, tmax
  {
    x = xg[i];
    for(int i2=0; i2<ny; i2++)
    {
      y = yg[i2];
      for(int i3=0; i3<nz; i3++)
      {
        z = zg[i3];
        for(int i4=0; i4<nt; i4=i4+nt-1) // i4=0, i4=nt-1
        {
          t = tg[i4];
          kg[S(i,i2,i3,i4,ny,nz,nt)][S(i,i2,i3,i4,ny,nz,nt)] = 1.0;
          fg[S(i,i2,i3,i4,ny,nz,nt)] = uana(x,y,z,t);
	} // i4
      } // i3
    } // i2
  } // i
  for(int i=0; i<nx; i++) // all x,y,t at zmin, zmax
  {
    x = xg[i];
    for(int i2=0; i2<ny; i2++)
    {
      y = yg[i2];
      for(int i4=0; i4<nt; i4++)
      {
        t = tg[i4];
        for(int i3=0; i3<nz; i3=i3+nz-1) // i3=0, i3=nz-1
        {
          z = zg[i3];
          kg[S(i,i2,i3,i4,ny,nz,nt)][S(i,i2,i3,i4,ny,nz,nt)] = 1.0;
          fg[S(i,i2,i3,i4,ny,nz,nt)] = uana(x,y,z,t);
	} // i4
      } // i3
    } // i2
  } // i
  for(int i=0; i<nx; i++) // all x,z,t at ymin, ymax
  {
    x = xg[i];
    for(int i3=0; i3<nz; i3++)
    {
      z = zg[i3];
      for(int i4=0; i4<nt; i4++)
      {
        t = tg[i4];
        for(int i2=0; i2<ny; i2=i2+ny-1) // i2=0, i2=ny-1
        {
          y = yg[i2];
          kg[S(i,i2,i3,i4,ny,nz,nt)][S(i,i2,i3,i4,ny,nz,nt)] = 1.0;
          fg[S(i,i2,i3,i4,ny,nz,nt)] = uana(x,y,z,t);
	} // i2
      } // i4
    } // i3
  } // i
  for(int i2=0; i2<ny; i2++) // all y,z,t at xmin, xmax
  {
    y = yg[i2];
    for(int i3=0; i3<nz; i3++)
    {
      z = zg[i3];
      for(int i4=0; i4<nt; i4++)
      {
        t = tg[i4];
        for(int i=0; i<nx; i=i+nx-1) // i=0, i=nx-1
        {
          x = xg[i];
          kg[S(i,i2,i3,i4,ny,nz,nt)][S(i,i2,i3,i4,ny,nz,nt)] = 1.0;
          fg[S(i,i2,i3,i4,ny,nz,nt)] = uana(x,y,z,t);
        } // i
      } // i4
    } // i3
  } // i2
  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;
  }
  for(int k=0; k<=npt; k++)
  {
    tt[k] = 0.0;
    wt[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");
  new gaulegf(tmin, tmax, tt, wt, npt); // tt and wt set up for integrations 
  System.out.print("tt[1]="+tt[1]+", tt[2]="+tt[2]+
                 ", wt[1]="+wt[1]+", wt[2]="+wt[2]+"\n");
  System.out.print("calling F( xmax="+xmax+", ymax="+ymax+",zmax="+zmax+
                   ",tmax="+tmax+")\n");
  val = F(xmax, ymax, zmax, tmax);
  System.out.print("F(xmax,ymax,zmax,tmax)="+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],tt[2],1,1,1,1,1,1,1,1,
             ccxx,ccyy,cczz,cctt,ccxy,ccxz,ccxt,ccyz,ccyt,cczt,
             ccx,ccy,ccz,cct,cc,
             nx,ny,nz,nt, xg,yg,zg,tg);
  System.out.print("galk(xx[2],yy[2])="+val+"\n");
  val = galf(xx[2],yy[2],zz[2],tt[2],1,1,1,1,nx,ny,nz,nt,xg,yg,zg,tg);
  System.out.print("galf(xx[2],yy[2],zz[2],tt[2])="+val+"\n");

  now = System.currentTimeMillis()/1000.0;
  System.out.println("setup time "+(now-prev)+" seconds");
  prev = now;

  // compute entries in stiffness matrix 
  System.out.print("compute stiffness matrix  \n");
  for(int i=1; i<nx-1; i++)
  {
    for(int i2=1; i2<ny-1; i2++)
    {
      for(int i3=1; i3<nz-1; i3++)
      {
        for(int i4=1; i4<nt-1; i4++)
        {
          for(int j=0; j<nx; j++)
          {
            for(int j2=0; j2<ny; j2++)
            {
              for(int j3=0; j3<nz; j3++)
	      {
                for(int j4=0; j4<nt; j4++)
                {
                  // compute integral for k(i,i2,i3,i4)  
                  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,i2,i3,i4,j,j2,j3,j4,
                                     ccxx,ccyy,cczz,cctt,ccxy,ccxz,
                                     ccxt,ccyz,ccyt,cczt,
                                     ccx,ccy,ccz,cct,cc,
                                     nx,ny,nz,nt,xg,yg,zg,tg);
			} // end pt

  if(i==1 && i2==1 && i3==1 && i4==1 && j==0 && j2==0 && j3==0 && j4==0 &&
     px==1 && py==1 && pz==1) 
  {
    now = System.currentTimeMillis()/1000.0;
    System.out.println("pt time "+(now-prev)+" seconds");
    prev = now;
  }
		      } // end pz
                    } // end py
                  } // end px
                  if(debug)
                  System.out.print("Legendre integration="+val+", at i="+i+
                  ", j="+j+", i2="+i2+", j2="+j2+", i3="+i3+", j3="+j3+"\n");
                  kg[S(i,i2,i3,i4,ny,nz,nt)][S(j,j2,j3,j4,ny,nz,nt)] = val;

  if(i==1 && i2==1 && i3==1 && i4==1 && j==0 && j2==0 && j3==0 && j4==0) 
  {
    now = System.currentTimeMillis()/1000.0;
    System.out.println("first integration time "+(now-prev)+" seconds");
    prev = now;
  }

		} // end j4
	      } // end j3
            } // end j2 
          } // end j

  if(i==1 && i2==1 && i3==1 && i4==1) 
  {
    now = System.currentTimeMillis()/1000.0;
    System.out.println("first point time "+(now-prev)+" seconds");
    prev = now;
  }

	} // end i4
      } // end i3 
    } // end i2 
  } // end i 

  // compute integral for f(i)  
  for(int i=1; i<nx-1; i++)
  {
    for(int i2=1; i2<ny-1; i2++)
    {
      for(int i3=1; i3<nz-1; i3++)
      {
        for(int i4=1; i4<nt-1; i4++)
        {
          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,i2,i3,i4,nx,ny,nz,nt,xg,yg,zg,tg);
		} // end pt
              } // end pz
            } // end py
          } // end px
          if(debug)
            System.out.print("Legendre integration="+val+", f at i="+i+
                  ", i2="+i2+", i3="+i3+", i4="+i4+"\n");
          fg[S(i,i2,i3,i4,ny,nz,nt)] = val;
	} // end i4
      } // end i3
    } // end i2
  } // end i
  System.out.print("k computed stiffness matrix \n");
  System.out.print("f computed forcing function \n");

  now = System.currentTimeMillis()/1000.0;
  System.out.println("system setup time "+(now-prev)+" seconds");
  prev = now;

  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 i2=0; i2<ny; i2++)
    {
      for(int i3=0; i3<nz; i3++)
      {
        for(int i4=0; i4<nt; i4++)
        {
          err = Math.abs(ug[S(i,i2,i3,i4,ny,nz,nt)]-
                         Ua[S(i,i2,i3,i4,ny,nz,nt)]);
          if(err>maxerr) maxerr = err;
          avgerr = avgerr + err;
          System.out.print("ug["+i+","+i2+","+i3+"]="+
                           ug[S(i,i2,i3,i4,ny,nz,nt)]+
                           ", Ua="+Ua[S(i,i2,i3,i4,ny,nz,nt)]+", err="+
                           (ug[S(i,i2,i3,i4,ny,nz,nt)]-
                           Ua[S(i,i2,i3,i4,ny,nz,nt)])+"\n");
	} // i4
      } // i3
    } // i2
  } // i
  System.out.print("                  maxerr="+maxerr+", avgerr="+
                 (avgerr/(nxyzt))+"\n");

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

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

// from uana23_fem.js  problem definition
double F(double x, double y, double z, double t) // RHS  
{
    return 15.0*(x*x+y*y+z*z+t*t+x*y+x*z+x*t+y*z+y*t+z*t)+
	 61.0*x + 62.0*y + 63.0*z + 64.0*t + 80.0;
} // end f

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

// needs laphi.js available
// for finite element method (not needed if coefficients available)
double galk(double x, double y, double z, double t,
            int i, int i2, int i3, int i4,
            int j, int j2, int j3, int j4,
            double ccxx, double ccyy, double cczz, double cctt,
            double ccxy, double ccxz, double ccxt, double ccyz,  
            double ccyt, double cczt, double ccx,  double ccy,
            double ccz,  double cct,  double cc,
            int nx, int ny, int nz, int nt,
            double xg[], double yg[], double zg[], double tg[])
{           // Galerkin k stiffness function for this problem 
  return (ccxx * L.phipp(x,j ,nx-1,xg)*  L.phi(y,j2,ny-1,yg)*
                   L.phi(z,j3,nz-1,zg)*  L.phi(t,j4,nt-1,tg)+
          ccyy *   L.phi(x,j ,nx-1,xg)*L.phipp(y,j2,ny-1,yg)*
                   L.phi(z,j3,nz-1,zg)*  L.phi(t,j4,nt-1,tg)+
          cczz *   L.phi(x,j ,nx-1,xg)*  L.phi(y,j2,ny-1,yg)*
                 L.phipp(z,j3,nz-1,zg)*  L.phi(t,j4,nt-1,tg)+
          cctt *   L.phi(x,j ,nx-1,xg)*  L.phi(y,j2,ny-1,yg)*
                   L.phi(z,j3,nz-1,zg)*L.phipp(t,j4,nt-1,tg)+
          ccxy *  L.phip(x,j ,nx-1,xg)* L.phip(y,j2,ny-1,yg)*
                   L.phi(z,j3,nz-1,zg)*  L.phi(t,j4,nt-1,tg)+
          ccxz *  L.phip(x,j ,nx-1,xg)*  L.phi(y,j2,ny-1,yg)*
                  L.phip(z,j3,nz-1,zg)*  L.phi(t,j4,nt-1,tg)+
          ccxt *  L.phip(x,j ,nx-1,xg)*  L.phi(y,j2,ny-1,yg)*
                   L.phi(z,j3,nz-1,zg)* L.phip(t,j4,nt-1,tg)+
          ccyz *   L.phi(x,j ,nx-1,xg)* L.phip(y,j2,ny-1,yg)*
                  L.phip(z,j3,nz-1,zg)*  L.phi(t,j4,nt-1,tg)+
          ccyt *   L.phi(x,j ,nx-1,xg)* L.phip(y,j2,ny-1,yg)*
                   L.phi(z,j3,nz-1,zg)* L.phip(t,j4,nt-1,tg)+
          cczt *   L.phi(x,j ,nx-1,xg)*  L.phi(y,j2,ny-1,yg)*
                  L.phip(z,j3,nz-1,zg)* L.phip(t,j4,nt-1,tg)+
          ccx  *  L.phip(x,j ,nx-1,xg)*  L.phi(y,j2,ny-1,yg)*
                   L.phi(z,j3,nz-1,zg)*  L.phi(t,j4,nt-1,tg)+
          ccy  *   L.phi(x,j ,nx-1,xg)* L.phip(y,j2,ny-1,yg)*
                   L.phi(z,j3,nz-1,zg)*  L.phi(t,j4,nt-1,tg)+
          ccz  *   L.phi(x,j ,nx-1,xg)*  L.phi(y,j2,ny-1,yg)*
                  L.phip(z,j3,nz-1,zg)*  L.phi(t,j4,nt-1,tg)+
          cct  *   L.phi(x,j ,nx-1,xg)*  L.phi(y,j2,ny-1,yg)*
                   L.phi(z,j3,nz-1,zg)* L.phip(t,j4,nt-1,tg)+
          cc   *   L.phi(x,j ,nx-1,xg)*  L.phi(y,j2,ny-1,yg)*
                   L.phi(z,j3,nz-1,zg)*  L.phi(t,j4,nt-1,tg))*
                   L.phi(x,i ,nx-1,xg)*  L.phi(y,i2,ny-1,yg)*
                   L.phi(z,i3,nz-1,zg)*  L.phi(t,i4,nt-1,tg);
} // end galk used by integration 

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

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

} // end class pde24js_la