// boltzmann_7la.java  Galerkin FEM
//
// solve one version of Boltzmann Equation in 7 independent variables
// I call it 7 dimensions, x,y,z in 3 space, t time domain,
// p,q,r momentum in x,y,z directions usually called px, py, pz
//
// I call the desired solution U(x,y,z,t,p,q,r) here, typically
// the probability of molecules f(x,y,z,t,px,py,pz) being
// in phase space volume dx,dy,dz  at x,y,z 
// in momentum space dpx,dpy,dpz  at px,py,pz at time t
//  
// partial derivative with respect to x is Ux, etc
// Ux(x,y,z,t,p,q,r) + Uy(x,y,z,t,p,q,r) + Uz(x,y,z,t,p,q,r) +
// Ut(x,y,z,t,p,q,r) + Sp*Up(x,y,z,t,p,q,r) + Sq*Uq(x,y,z,t,p,q,r) +
// Sr*Ur(x,y,z,t,p,q,r)  = Sf(x,y,z,t,p,q,r)
// Where Sf is the forcing function, Sp, Sq, Sr probabilities
// taken as constant here, yet may be a function of x,y,z,t
//
// boundary conditions, computed using Ub(x,y,z,t,p,q,r)
// were chosen, arbitrarily,
// for purpose of comparing different solution methods.
//
// gauss-legendre integration used
// solve for Uijkluvw numerical approximation of
// U(x_i,y_ii,z_iii,t_iiii,p_iiiii,q_iiiiii,r_iiiiiii)
// kg * ug = fg, stiffness matrix times unknown U equal forcing function
// solve simultaneous equations for U

import java.io.*;

class boltzmann_7la
{
  // i, j, nx, xg, ii, jj, ny, yg ... needed by galk or galf, available at call 
  final int nx = 3; // do not need to be the same
  final int ny = 3;
  final int nz = 3;
  final int nt = 3;
  final int np = 3;
  final int nq = 3;
  final int nr = 3;
  final int nxyztpqr = nx*ny*nz*nt*np*nq*nr;
  double xg[] = new double[nx];
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double tg[] = new double[nt];
  double pg[] = new double[np];
  double qg[] = new double[nq];
  double rg[] = new double[nr];
  int i, j, ii, jj, iii, jjj, iiii, jjjj;
  int iiiii, jjjjj, iiiiii, jjjjjj, iiiiiii, jjjjjjj;
  double x, y, z, t, p, q, r, hx, hy, hz, ht, hp, hq, hr;
  laphi L = new laphi();

  boltzmann_7la()
  {
    final int npx = 4; // Gauss Legendre integration order
    final int npy = 4; // do not need to be the same
    final int npz = 4;
    final int npt = 4;
    final int npp = 4;
    final int npq = 4;
    final int npr = 4;
    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 
    double pp[] = new double[npp+1];
    double wp[] = new double[npp+1]; // for Gauss-Legendre 
    double qq[] = new double[npq+1];
    double wq[] = new double[npq+1]; // for Gauss-Legendre 
    double rr[] = new double[npr+1];
    double wr[] = new double[npr+1]; // for Gauss-Legendre 
    double val, err, avgerr, maxerr;
    int px, py, pz, pt, pa, pq, pr; // pp used, thus pa

    double xmin = 0.1; // problem parameters 
    double xmax = 0.7;
    double ymin = 0.1;
    double ymax = 0.7;
    double zmin = 0.1;
    double zmax = 0.7;
    double tmin = 0.1;
    double tmax = 0.7;
    double pmin = 0.1;
    double pmax = 0.7;
    double qmin = 0.1;
    double qmax = 0.7;
    double rmin = 0.1;
    double rmax = 0.7;

    double kg[] = new double[nxyztpqr*nxyztpqr];  
    double fg[] = new double[nxyztpqr];
    double ug[] = new double[nxyztpqr];
    double Ua[] = new double[nxyztpqr];

    double t_start, t_init, t_kf, t_simeq, t_end;

    System.out.println("boltzmann_7la.c running");
    System.out.println("Given Ux(x,y,z,t,p,q,r)+Uy(x,y,z,t,p,q,r)+Uz(x,y,z,t,p,q,r)+");
    System.out.println("Ut(x,y,z,t,p,q,r)+Up(x,y,z,t,p,q,r)+Uq(x,y,z,t,p,q,t)+");
    System.out.println("Ur(x,y,z,t,p,q,r,z) = Sf(x,y,z,t)");
    System.out.println("Sf(x,y,z,t,p,q,r) = RHS");
    System.out.println("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax ");
    System.out.println("tmin<=t<=tmax pmin<=p<=pmax qmin<=q<=qmax ");
    System.out.println("rmin<=r<=rmax  Boundaries");
    System.out.println("Analytic solution U(x,y,z,t,p,q,r) = x*y*z*t*p*q*r ");
    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("pmin="+pmin+", pmax="+pmax+", qmin="+qmin+", qmax="+qmax);
    System.out.println("rmin="+rmin+", rmax="+rmax);
    System.out.println("nx="+nx+", ny="+ny+", nz="+nz+", nt="+nt);
    System.out.println("np="+np+", nq="+nq+", nr="+nr);
    System.out.println(" ");
    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] = Ub(xg[i], ymin, zmin, tmin, pmin, qmin, rmin); // 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] = Ub(xmin, yg[ii], zmin, tmin, pmin, qmin, rmin); // 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] = Ub(xmin, ymin, zg[iii], tmin, pmin, qmin, rmin); // 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] = Ub(xmin, ymin, zmin, tg[iiii], pmin, qmin, rmin); // temp 
      System.out.println("iiii="+iiii+", Ua("+tg[iiii]+")="+Ua[iiii]);
    } // iiii
    System.out.println(" ");
    System.out.println("p grid and analytic solution at xmin,ymin,zmin ");
    hp = (pmax-pmin)/(double)(np-1);
    for(int iiiii=0; iiiii<np; iiiii++)
    {
      pg[iiiii] = pmin+iiiii*hp;
      Ua[iiiii] = Ub(xmin, ymin, zmin, tmin, pg[iiiii], qmin, rmin); // temp 
      System.out.println("iiiii="+iiiii+", Ua("+pg[iiiii]+")="+Ua[iiiii]);
    } // iiiii
    System.out.println(" ");
    System.out.println("q grid and analytic solution at xmin,ymin,zmin ");
    hq = (qmax-qmin)/(double)(nq-1);
    for(int iiiiii=0; iiiiii<nq; iiiiii++)
    {
      qg[iiiiii] = qmin+iiiiii*hq;
      Ua[iiiiii] = Ub(xmin, ymin, zmin, tmin, pmin, qg[iiiiii], rmin); // temp 
      System.out.println("iiiiii="+iiiiii+", Ua("+qg[iiiiii]+")="+Ua[iiiiii]);
    } // iiiiii
    System.out.println(" ");
    System.out.println("r grid and analytic solution at xmin,ymin,zmin ");
    hr = (rmax-rmin)/(double)(nr-1);
    for(int iiiiiii=0; iiiiiii<nr; iiiiiii++)
    {
      rg[iiiiiii] = rmin+iiiiiii*hr;
      Ua[iiiiiii] = Ub(xmin, ymin, zmin, tmin, pmin, qmin, rg[iiiiiii]); // temp 
      System.out.println("iiiiiii="+iiiiiii+", Ua("+rg[iiiiiii]+")="+Ua[iiiiiii]);
    } // iiiiiii
    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;
            for(int iiiii=0; iiiii<np; iiiii++)
            {
	      p = rmin+(double)iiiii*hp;
              for(int iiiiii=0; iiiiii<np; iiiiii++)
              {
	        q = rmin+(double)iiiiii*hq;
                for(int iiiiiii=0; iiiiiii<nq; iiiiiii++)
                {
	          r = rmin+(double)iiiiiii*hr;
          Ua[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)] = Ub(x,y,z,t,p,q,r);
		}
	      }
	    }
	  }
	}
      }
    }
    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++)
                  {
                    for(int iiiii=0; iiiii<np; iiiii++)
                    {
                      for(int jjjjj=0; jjjjj<np; jjjjj++)
                      {
                        for(int iiiiii=0; iiiiii<nq; iiiiii++)
                        {
                          for(int jjjjjj=0; jjjjjj<nq; jjjjjj++)
                          {
                            for(int iiiiiii=0; iiiiiii<nr; iiiiiii++)
                            {
                              for(int jjjjjjj=0; jjjjjjj<nr; jjjjjjj++)
                              {
		      kg[ss(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,
                            j,jj,jjj,jjjj,jjjjj,jjjjjj,jjjjjjj)] = 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++)
          {
            for(int iiiii=0; iiiii<np; iiiii++)
            {
              for(int iiiiii=0; iiiiii<nq; iiiiii++)
              {
                for(int iiiiiii=0; iiiiiii<nr; iiiiiii++)
                {
	          fg[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)] = 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];
            for(int iiiii=0; iiiii<np; iiiii++)
	    {
              p = pg[iiiii];
              for(int iiiiii=0; iiiiii<nq; iiiiii++)
	      {
                q = qg[iiiiii];
                for(int iiiiiii=0; iiiiiii<nr; iiiiiii++)
	        {
                  r = rg[iiiiiii];
            kg[ss(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,
                  i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)] = 1.0;
	    fg[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)] = Ub(x,y,z,t,p,q,r);
		} // iiiiiii
	      } // iiiiii
	    } // iiiii
	  } // 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 
    System.out.println("calling gauleg tmin="+tmin+", tmax="+tmax+
                       ", npt="+npt);
    new gaulegf(tmin, tmax, pp, wp, npp); // pp and wp set up for integrations 
    System.out.println("calling gauleg pmin="+pmin+", pmax="+pmax+
                       ", npp="+npp);
    new gaulegf(qmin, qmax, qq, wq, npq); // qq and wq set up for integrations 
    System.out.println("calling gauleg rmin="+rmin+", rmax="+rmax+
                       ", npr="+npr);
    new gaulegf(rmin, rmax, rr, wr, npr); // rr and wr 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]);
    System.out.println("pp[1]="+pp[1]+", pp[2]="+pp[2]+
                       ", wp[1]="+wp[1]+", wp[2]="+wp[2]);
    System.out.println("qq[1]="+qq[1]+", qq[2]="+qq[2]+
                       ", wq[1]="+wq[1]+", wq[2]="+wq[2]);
    System.out.println("rr[1]="+rr[1]+", rr[2]="+rr[2]+
                       ", wr[1]="+wr[1]+", wr[2]="+wr[2]);
    i=1; ii=1; iii=1; iiii=1; j=1; jj=1; jjj=1; jjjj=1; 
    iiii=1; iiiii=1; iiiiii=1; iiiiiii=1; jjjj=1; jjjjj=1; jjjjjj=1; jjjjjjj=1; 
    val = galk(xx[2],yy[2],zz[2],tt[2],pp[2],qq[2],rr[2]);
    System.out.println("galk(xx[2],yy[2],zz[2],tt[2],pp[2],qq[2],rr[2])="+val);
    val = galf(xx[2],yy[2],zz[2],tt[2],pp[2],qq[2],rr[2]);
    System.out.println("galf(xx[2],yy[2],zz[2],tt[2],pp[2],qq[2],rr[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  ");
    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++)
              {
                for(iiii=1; iiii<nt-1; iiii++)
                {
                  for(jjjj=0; jjjj<nt; jjjj++)
                  {
                    for(iiiii=1; iiiii<np-1; iiiii++)
                    {
                      for(jjjjj=0; jjjjj<np; jjjjj++)
                      {
                        for(iiiiii=1; iiiiii<nq-1; iiiiii++)
                        {
                          for(jjjjjj=0; jjjjjj<nq; jjjjjj++)
                          {
                            for(iiiiiii=1; iiiiiii<nr-1; iiiiiii++)
                            {
                              for(jjjjjjj=0; jjjjjjj<nr; jjjjjjj++)
                              {
              // compute integral for k(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)  
                                val = 0.0;
                                for(px=1; px<=npx; px++)
                                {
                                  for(py=1; py<=npy; py++)
                                  {
                                    for(pz=1; pz<=npz; pz++)
                                    {
                                      for(pt=1; pt<=npt; pt++)
                                      {
                                        for(pa=1; pa<=npp; pa++)
                                        {
                                          for(pq=1; pq<=npq; pq++)
                                          {
                                            for(pr=1; pr<=npr; pr++)
                                            {
      val = val + wx[px]*wy[py]*wz[pz]*wt[pt]*wp[pa]*wq[pq]*wr[pr]*
	     galk(xx[px],yy[py],zz[pz],tt[pt],pp[pa],qq[pq],rr[pr]);
			                    } // pr
		                          } // pq
					} // pa
				      } // pt
				    } // pz
		                  } // py
		                } // px
                    kg[ss(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,
                          j,jj,jjj,jjjj,jjjjj,jjjjjj,jjjjjjj)] = val;
		              } // end jjjjjjj loop
		            } // end iiiiiii loop
		          } // end jjjjjj loop
	         	} // end iiiiii loop
		      } // end jjjjj loop
	            } // end iiiii loop
		  } // 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="+ii+", jj="+jj);
        System.out.println("        , iii="+iii+", jjj="+jjj+
                           ", iiii="+iiii+", jjjj="+jjjj);
        System.out.println("        , iiiii="+iiiii+", jjjjj="+jjjjj+
                           ", iiiiii="+iiiiii+", jjjjjj="+jjjjjj);
        System.out.println("        , iiiiiii="+iiiiiii+", jjjjjjj="+jjjjjjj);
      } // end j loop 
    } // end i loop 

    // compute integral for f(i)  
    for(i=1; i<nx-1; i++)
    {
      for(ii=1; ii<ny-1; ii++)
      {
	for(iii=1; iii<nz-1; iii++)
	{
	  for(iiii=1; iiii<nt-1; iiii++)
	  {
	    for(iiiii=1; iiiii<np-1; iiiii++)
	    {
	      for(iiiiii=1; iiiiii<nq-1; iiiiii++)
	      {
	        for(iiiiiii=1; iiiiiii<nr-1; iiiiiii++)
	        {
                  val = 0.0;
                  for(px=1; px<=npx; px++)
                  {
                    for(py=1; py<=npy; py++)
                    {
	              for(pz=1; pz<=npz; pz++)
	              {
	                for(pt=1; pt<=npt; pt++)
	                {
	                  for(pa=1; pa<=npp; pa++)
	                  {
	                    for(pq=1; pq<=npq; pq++)
	                    {
	                      for(pr=1; pr<=npr; pr++)
	                      {
         val = val + wx[px]*wy[py]*wz[pz]*wt[pt]*wp[pa]*wq[pq]*wr[pr]*
	        galf(xx[px],yy[py],zz[pz],tt[pt],pp[pa],qq[pq],rr[pr]);
			      } // pr
		            } // pq
		          } // pa
		        } // pt
	              } // pz
	            } // py
                  } //px
                  fg[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)] = val;
		} // end iiiiiii loop
	      } // end iiiiii loop
	    } // end iiiii loop
	  } // end iiii loop
        } // end iii loop
      } // end ii loop
      System.out.println("Legendre integration="+val+", f at i="+i+
                         ", ii="+ii+", iii="+iii+", iiii="+iiii+
                         ", iiiii="+iiiii+", iiiiii="+iiiiii+
                         ", iiiiiii="+iiiiiii);
    } // end i loop
    System.out.println("k computed stiffness matrix, see above  ");
    System.out.println("f computed forcing function, see above  ");
    t_kf = System.currentTimeMillis();
    System.out.println("kf build took "+((t_kf-t_init)/1000.0)+" seconds");

    System.out.println("solve "+nxyztpqr+" equations in "+nxyztpqr+" unknowns");
    new simeq(nxyztpqr, 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++)
          {
            for(int iiiii=0; iiiii<np; iiiii++)
            {
              for(int iiiiii=0; iiiiii<nq; iiiiii++)
              {
                for(int iiiiiii=0; iiiiiii<nr; iiiiiii++)
                {
	err = Math.abs(ug[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)]-
                       Ua[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)]);
        if(err>maxerr) maxerr = err;
        avgerr = avgerr + err;
        System.out.println("ug["+i+","+ii+","+iii+","+iiii+
		           ","+iiiii+","+iiiiii+","+iiiiiii+"]="+
		           ug[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)]+
		           ", Ua="+Ua[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)]+
		           ", err="+(ug[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)]-
                           Ua[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)]));
		}
	      }
	    }
	  }
	}
      }
    }
    System.out.println("                         maxerr="+maxerr+", avgerr="+
                       avgerr/(double)(nxyztpqr));

    write_soln(ug);
    t_end = System.currentTimeMillis();
    System.out.println("total took "+((t_end-t_start)/1000.0)+" seconds");

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

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

  int ss(int i, int ii, int iii, int iiii, int iiiii, int iiiiii, int iiiiiii,
         int j, int jj, int jjj, int jjjj, int jjjjj, int jjjjjj, int jjjjjjj)
  {
    return s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)*nxyztpqr+
	   s(j,jj,jjj,jjjj,jjjjj,jjjjjj,jjjjjjj);
  } // end ss

  // PDE problem definition functions
    double Sf(double x, double y, double z, double t,
              double p, double q, double r)
                                           // forcing function, F(x,y,z,t)
  {
    return y*z*t*p*q*r +
           x*z*t*p*q*r +
           x*y*t*p*q*r +
           x*y*z*p*q*r +
           x*y*z*t*p*r +
	   x*y*z*t*p*q;
  }

    double Ub(double x, double y, double z, double t,
              double p, double q, double r)
                                      // analytic solution and boundaries
  {
    return x*y*z*p*q*r*t;
  }

    double galk(double x, double y, double z, double t,
                double p, double q, double r)
                       // Galerkin k stiffness function for this problem
  {
    return
        (L.phip(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(p,jjjjj,np-1,pg)*L.phi(q,jjjjjj,nq-1,qg)*
         L.phi(r,jjjjjjj,nr-1,rg)+
         L.phi(x,j,nx-1,xg)*L.phip(y,jj,ny-1,yg)*L.phi(z,jjj,nz-1,zg)*
         L.phi(t,jjjj,nt-1,tg)*L.phi(p,jjjjj,np-1,pg)*L.phi(q,jjjjjj,nq-1,qg)*
         L.phi(r,jjjjjjj,nr-1,rg)+
         L.phi(x,j,nx-1,xg)*L.phi(y,jj,ny-1,yg)*L.phip(z,jjj,nz-1,zg)*
	 L.phi(t,jjjj,nt-1,tg)*L.phi(p,jjjjj,np-1,pg)*L.phi(q,jjjjjj,nq-1,qg)*
         L.phi(r,jjjjjjj,nr-1,rg)+
         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)*L.phi(p,jjjjj,np-1,pg)*L.phi(q,jjjjjj,nq-1,qg)*
         L.phi(r,jjjjjjj,nr-1,rg)+
         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.phip(p,jjjjj,np-1,pg)*L.phi(q,jjjjjj,nq-1,qg)*
         L.phi(r,jjjjjjj,nr-1,rg)+
         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(p,jjjjj,np-1,pg)*L.phip(q,jjjjjj,nq-1,qg)*
         L.phi(r,jjjjjjj,nr-1,rg)+
         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(p,jjjjj,np-1,pg)*L.phi(q,jjjjjj,nq-1,qg)*
         L.phip(r,jjjjjjj,nr-1,rg))*
         (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)*L.phi(p,iiiii,np-1,pg)*L.phi(q,iiiiii,nq-1,qg)*
          L.phi(r,iiiiiii,nr-1,rg));
  } // end galk used by integration

  double galf(double x, double y, double z, double t,
              double p, double q, double r)
                                 // Galerkin f function for this problem
  {
    return Sf(x,y,z,t,p,q,r)*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)*
           L.phi(p,iiiii,np-1,pg)*L.phi(q,iiiiii,nq-1,qg)*
           L.phi(r,iiiiiii,nr-1,rg);
  } // end galf used by integration

  void write_soln(double Usol[])
  {
    String file_name = "boltzmann_7la.dat";
    try
    {
      FileWriter fout = new FileWriter(file_name);
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing "+file_name);
      fileout.write("7 "+nx+" "+ny+" "+nz+" "+nt+
                   " "+np+" "+nq+" "+nr+" \n");

      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++)
            {
              for(int iiiii=0; iiiii<np; iiiii++)
              {
                for(int iiiiii=0; iiiiii<nq; iiiiii++)
                {
                  for(int iiiiiii=0; iiiiiii<nr; iiiiiii++)
                  {
    fileout.write(xg[i]+" "+yg[ii]+" "+zg[iii]+" "+tg[iiii]+" "+
                  pg[iiiii]+" "+qg[iiiiii]+" "+rg[iiiiiii]+" "+
                  Usol[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)]+" \n");
                  } // end iiiiiii
                } // end iiiiii
              } // end iiiii
            } // end iiii
          } // end iii
        } // end ii
      } // end i

      fileout.close();
      System.out.println("finished writing nu2_laplace.dat");
    } // end try
    catch(FileNotFoundException exception)
    {
      System.out.println(file_name+"  not opened");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println(file_name+" output-file error");
    } // end try...catch
  } // end write_soln

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