// fem_check47h_la.java  Galerkin FEM
//                    Biharmonic 4th order PDE in seven dimensions
// solve Uxxxx(x,y,z,t,u,v,w) + Uyyyy(x,y,z,t,u,v,w) + Uzzzz(x,y,z,t,u,v,w) + 
//       Utttt(x,y,z,t,u,v,w) + Uuuuu(x,y,z,t,u,v,w) + Uvvvv(x,y,z,t,u,v,w) + 
//       Uwwww(x,y,z,t,u,v,w) + 
//       2 Uxx(x,y,z,t,u,v,w) + 2 Uyy(x,y,z,t,u,v,w) + 2 Uzz(x,y,z,t,u,v,w) +
//       2 Utt(x,y,z,t,u,v,w) + 2 Uuu(x,y,z,t,u,v,w) + 2 Uvv(x,y,z,t,u,v,w) + 
//       2 Uww(x,y,z,t,u,v,w) + 
//       8 U(x,y,z,t,u,v,w)   = F(x,y,z,t,u,v,w)
//
//       F(x,y,z,t,u,v,w) = sin(x+y+z+t+u+v+w)
//
// boundary conditions computed using uana(x,y,z,t,u,v,w)
// analytic solution is U(x,y,z,t,u,v,w) = sin(x+y+z+t+u+v+w)
//
// subscripts in code use x[i], y[i2], z[i3], t[i4], u[i5], v[i6], w[i7]
//
// using Lagrange Phi functions and Gauss-Legendre integration
// solve for U numerical approximation U(x_i,y_i2,z_i3,t_i4,u_i5,v_i6,w_i7)
// kg * uu = fg, solve simultaneous equations for U
//

class fem_check47h_la
{
  final int nx = 3;
  final int ny = 3;
  final int nz = 3;
  final int nt = 3;
  final int nu = 3;
  final int nv = 3;
  final int nw = 3;
  final int nxyztuvw = nx*ny*nz*nt*nu*nv*nt;
  double xg[] = new double[nx];
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double tg[] = new double[nt];
  double ug[] = new double[nu];
  double vg[] = new double[nv];
  double wg[] = new double[nw];
  int i, j, i2, j2, i3, j3, i4, j4, i5, j5, i6, j6, i7, j7;
  laphi L = new laphi();

  fem_check47h_la()
  {
    double x, y, z, t, u, v, w, hx, hy, hz, ht, hu, hv, hw;
    final int npx = 8; // Gauss Legendre integration order
    final int npy = 8;
    final int npz = 8;
    final int npt = 8;
    final int npu = 8;
    final int npv = 8;
    final int npw = 8;
    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 uuu[] = new double[npu+1];
    double wu[] = new double[npu+1];
    double vv[] = new double[npv+1];
    double wv[] = new double[npv+1];
    double www[] = new double[npw+1];
    double ww[] = new double[npw+1];
    double val, err, avgerr, maxerr;
    int px, py, pz, pt, pu, pv, pw;

    double xmin = 0.0; // problem parameters 
    double xmax = 0.5; 
    double ymin = 0.0;
    double ymax = 0.5; 
    double zmin = 0.0;
    double zmax = 0.5; 
    double tmin = 0.0;
    double tmax = 0.5; 
    double umin = 0.0;
    double umax = 0.5; 
    double vmin = 0.0;
    double vmax = 0.5; 
    double wmin = 0.0;
    double wmax = 0.5; 

    double kg[][] = new double[nxyztuvw][nxyztuvw];  
    double fg[] = new double[nxyztuvw];
    double uu[] = new double[nxyztuvw];
    double Ua[] = new double[nxyztuvw];
    double tstart, tnow;

    tstart = System.currentTimeMillis();
    System.out.println("fem_check47h_la.java running ");
    System.out.println("Given Uxxxx(x,y,z,t,u,v,w) + Uyyyy(x,y,z,t,u,v,w) +");
    System.out.println("      Uzzzz(x,y,z,t,u,v,w) + Utttt(x,y,z,t,u,v,w) +");
    System.out.println("      Uuuuu(x,y,z,t,u,v,w) + Uvvvv(x,y,z,t,u,v,w) +");
    System.out.println("      Uwwww(x,y,z,t,u,v,w) + ");
    System.out.println("      2 Uxx(x,y,z,t,u,v,w) + 2 Uyy(x,y,z,t,u,v,w) +");
    System.out.println("      2 Uzz(x,y,z,t,u,v,w) + 2 Utt(x,y,z,t,u,v,w) +");
    System.out.println("      2 Uuu(x,y,z,t,u,v,w) + 2 Uvv(x,y,z,t,u,v,w) +");
    System.out.println("      2 Uww(x,y,z,t,u,v,w) + ");
    System.out.println("      8 U(x,y,z,t,u,v,w)   = F(x,y,z,t,u,v,w)");
    System.out.println("F(x,y,z,t,u,v,w) = sin(x+y+z+t+u+v+w)");
    System.out.println(" ");
    System.out.println("Biharmonic 4th order PDE in seven dimensions");
    System.out.println("Solve by  Galerkin Finite Element Method");
    System.out.println("using Lagrange Phi functions and Gauss-Legendre integration");

    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("umin="+umin+", umax="+umax+", vmin="+vmin+
                       ", vmax="+vmax);
    System.out.println("wmin="+wmin+", wmax="+wmax);
    System.out.println("nx="+nx+", ny="+ny+", nz="+nz+", nt="+nt);
    System.out.println("nu="+nu+", nv="+nv+", nw="+nw);
    System.out.println("x grid and analytic solution ");
    hx = (xmax-xmin)/(double)(nx-1);
    for(i=0; i<nx; i++)
    {
      xg[i] = xmin+i*hx;
      val = uana(xg[i],ymin,zmin,tmin,umin,vmin,wmin);
      System.out.println("i="+i+", U("+xg[i]+
      ","+ymin+","+zmin+","+tmin+","+umin+","+vmin+","+wmin+")="+val);
    } // end i  
    System.out.println(" ");
    System.out.println("y grid and analytic solution ");
    hy = (ymax-ymin)/(double)(ny-1);
    for(i2=0; i2<ny; i2++)
    {
      yg[i2] = ymin+i2*hy;
      val = uana(xmax, yg[i2], zmin, tmin, umin, vmin, wmin);
      System.out.println("i2="+i2+", U("+xmax+","+yg[i2]+
      ","+zmin+","+tmin+","+umin+","+vmin+","+wmin+")="+val);
    } // end i2
    System.out.println(" ");
    System.out.println("z grid and analytic solution ");
    hz = (zmax-zmin)/(double)(nz-1);
    for(i3=0; i3<nz; i3++)
    {
      zg[i3] = zmin+i3*hz;
      val = uana(xmax, ymax, zg[i3], tmin, umin, vmin, wmin);
      System.out.println("i3="+i3+", U("+xmax+","+ymax+","+zg[i3]+
      ","+tmin+","+umin+","+vmin+","+wmin+")="+val);
    } // end i3  
    System.out.println(" ");
    System.out.println("t grid and analytic solution ");
    ht = (tmax-tmin)/(double)(nt-1);
    for(i4=0; i4<nt; i4++)
    {
      tg[i4] = tmin+i4*ht;
      val = uana(xmax, ymax, zmax, tg[i4], umin, vmin, wmin);
      System.out.println("i4="+i4+", U("+xmax+","+ymax+","+zmax+
      ","+tg[i4]+","+umin+","+vmin+","+wmin+")="+val);
    } // end i4  
    System.out.println(" ");
    System.out.println("u grid and analytic solution ");
    hu = (umax-umin)/(double)(nu-1);
    for(i5=0; i5<nu; i5++)
    {
      ug[i5] = umin+i5*hu;
      val = uana(xmax, ymax, zmax, tmax, ug[i5], vmin, wmin);
      System.out.println("i5="+i5+", U("+xmax+","+ymax+","+zmax+
      ","+tmax+","+ug[i5]+","+vmin+","+wmin+")="+val);
    } // end i5  
    System.out.println(" ");
    System.out.println("v grid and analytic solution ");
    hv = (vmax-vmin)/(double)(nv-1);
    for(i6=0; i6<nv; i6++)
    {
      vg[i6] = vmin+i6*hv;
      val = uana(xmax, ymax, zmax, tmax, umax, vg[i6], wmin);
      System.out.println("i6="+i6+", U("+xmax+","+ymax+","+zmax+
      ","+tmax+","+umax+","+vg[i6]+","+wmin+")="+val);
    } // end i6  
    System.out.println(" ");
    System.out.println("w grid and analytic solution ");
    hw = (wmax-wmin)/(double)(nw-1);
    for(i7=0; i7<nw; i7++)
    {
      wg[i7] = wmin+i7*hw;
      val = uana(xmax, ymax, zmax, tmax, umax, vmax, wg[i7]);
      System.out.println("i7="+i7+", U("+xmax+","+ymax+","+zmax+
      ","+tmax+","+umax+","+vmax+","+wg[i7]+")="+val);
    } // end i7  
    System.out.println(" ");

    // put solution in Ua vector 
    for(i=0; i<nx; i++)
    {
      for(i2=0; i2<ny; i2++)
      {
        for(i3=0; i3<nz; i3++)
        {
          for(i4=0; i4<nt; i4++)
          {
            for(i5=0; i5<nu; i5++)
            {
              for(i6=0; i6<nv; i6++)
              {
                for(i7=0; i7<nw; i7++)
                {
                  Ua[s(i,i2,i3,i4,i5,i6,i7)] = 
                      uana(xg[i],yg[i2],zg[i3],tg[i4],ug[i5],vg[i6],wg[i7]);
		}
	      }
	    }
	  }
	}
      }
    }
    System.out.println(" ");

    System.out.println("initialize Galerkin k stiffness matrix and f RHS");
    // initialize k and f 
    for(i=0; i<nx; i++)
    {
      for(j=0; j<nx; j++)
      {
        for(i2=0; i2<ny; i2++)
        {
          for(j2=0; j2<ny; j2++)
          {
            for(i3=0; i3<nz; i3++)
            {
              for(j3=0; j3<nz; j3++)
              {
                for(i4=0; i4<nt; i4++)
                {
                  for(j4=0; j4<nt; j4++)
                  {
                    for(i5=0; i5<nu; i5++)
                    {
                      for(j5=0; j5<nu; j5++)
                      {
                        for(i6=0; i6<nv; i6++)
                        {
                          for(j6=0; j6<nv; j6++)
                          {
                            for(i7=0; i7<nw; i7++)
                            {
                              for(j7=0; j7<nw; j7++)
                              {
		    kg[s(i,i2,i3,i4,i5,i6,i7)][s(j,j2,j3,j4,j5,j6,j7)] = 0.0;
			      }
			    }
			  }
			}
		      }
		    }
		  }
		}
	      }
	    }
	  }
        }
      }
    }
    for(i=0; i<nx; i++)
    {
      for(i2=0; i2<ny; i2++)
      {
	for(i3=0; i3<nz; i3++)
	{
	  for(i4=0; i4<nt; i4++)
	  {
	    for(i5=0; i5<nu; i5++)
	    {
	      for(i6=0; i6<nv; i6++)
	      {
		for(i7=0; i7<nw; i7++)
		{
		  fg[s(i,i2,i3,i4,i5,i6,i7)] = 0.0;
		}
	      }
	    }
	  }
	}
      }
    }

    System.out.println("incorporate boundary conditions");
    // set boundary conditions, all faces 
    System.out.println("set kg diagonal and fg");
    for(i=0; i<nx; i++)
    {
      x = xg[i];
      for(i2=0; i2<ny; i2++)
      {
	y = yg[i2];
	for(i3=0; i3<nz; i3++)
	{
	  z = zg[i3];
	  for(i4=0; i4<nt; i4++)
	  {
	    t = tg[i4];
	    for(i5=0; i5<nu; i5++)
	    {
	      u = ug[i5];
	      for(i6=0; i6<nv; i6++)
	      {
		v = vg[i6];
		for(i7=0; i7<nw; i7++)
		{
		  w = wg[i7];
		  kg[s(i,i2,i3,i4,i5,i6,i7)][s(i,i2,i3,i4,i5,i6,i7)] = 1.0;
	          val = uana(x,y,z,t,u,v,w);
	          fg[s(i,i2,i3,i4,i5,i6,i7)] = val;
		}
	      }
	    }
	  }
	}
      }
    }

    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]);
    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("zz[1]="+zz[1]+", zz[2]="+zz[2]+
                       ", wz[1]="+wz[1]+", wz[2]="+wz[2]);
    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("tt[1]="+tt[1]+", tt[2]="+tt[2]+
                       ", wt[1]="+wt[1]+", wt[2]="+wt[2]);
    System.out.println("calling gauleg umin="+umin+", umax="+umax+
                       ", npu="+npu);
    new gaulegf(umin, umax, uuu, wu, npu); // uuu and wu set up for integrations 
    System.out.println("uuu[1]="+uuu[1]+", uuu[2]="+uuu[2]+
                       ", wu[1]="+wu[1]+", wu[2]="+wu[2]);
    System.out.println("calling gauleg vmin="+vmin+", vmax="+tmax+
                       ", npv="+npv);
    new gaulegf(vmin, vmax, vv, wv, npv); // tt and wt set up for integrations 
    System.out.println("vv[1]="+vv[1]+", vv[2]="+vv[2]+
                       ", wv[1]="+wv[1]+", wv[2]="+wv[2]);
    System.out.println("calling gauleg wmin="+wmin+", wmax="+wmax+
                       ", npw="+npw);
    new gaulegf(wmin, wmax, www, ww, npw); // www and ww set up for integrations 
    System.out.println("www[1]="+www[1]+", www[2]="+www[2]+
                       ", ww[1]="+ww[1]+", ww[2]="+ww[2]);

    i=1; i2=1; i3=1; i4=1; i5=1; i6=1; i7=1; // for print test check
    j=1; j2=1; j3=1; j4=1; j5=1; j6=1; j7=1;
    val = galk(xx[2],yy[2],zz[2],tt[2],uuu[2],vv[2],www[2]);
    System.out.println("galk(xx[2],yy[2],zz[2],tt[2],uuu[2],vv[2],www[2])="+val);
    val = galf(xx[2],yy[2],zz[2],tt[2],uuu[2],vv[2],www[2]);
    System.out.println("galf(xx[2],yy[2],zz[2],tt[2],uuu[2],vv[2],www[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(i2=1; i2<ny-1; i2++)
        {
          for(j2=0; j2<ny; j2++)
          {
            for(i3=1; i3<nz-1; i3++)
            {
              for(j3=0; j3<nz; j3++)
              {
                for(i4=1; i4<nt-1; i4++)
                {
                  for(j4=0; j4<nt; j4++)
                  {
                    for(i5=1; i5<nu-1; i5++)
                    {
                      for(j5=0; j5<nu; j5++)
                      {
                        for(i6=1; i6<nv-1; i6++)
                        {
                          for(j6=0; j6<nv; j6++)
                          {
                            for(i7=1; i7<nw-1; i7++)
                            {
                              for(j7=0; j7<nw; j7++)
                              {
                                // compute integral for k(i,j)  
                                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(pu=1; pu<=npu; pu++)
                                        {
                                          for(pv=1; pv<=npv; pv++)
                                          {
                                            for(pw=1; pw<=npw; pw++)
                                            {
                val = val + wx[px]*wy[py]*wz[pz]*wt[pt]*wu[pu]*wv[pv]*ww[pw]*
		    galk(xx[px],yy[py],zz[pz],tt[pt],uuu[pu],vv[pv],www[pw]);
					    }
					  }
					}
				      }
				    }
				  }
				}
                kg[s(i,i2,i3,i4,i5,i6,i7)][s(j,j2,j3,j4,j5,j6,j7)] = val;
			      }
			    }
			  }
			}
		      }
		    }
		  }
		}
	      }
            }
          }
        } 
      } 
    } 

    // compute integral for f(i)  
    for(i=1; i<nx-1; i++)
    {
      for(i2=1; i2<ny-1; i2++)
      {
        for(i3=1; i3<nz-1; i3++)
        {
          for(i4=1; i4<nt-1; i4++)
          {
            for(i5=1; i5<nu-1; i5++)
            {
              for(i6=1; i6<nv-1; i6++)
              {
                for(i7=1; i7<nw-1; i7++)
                {
                  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(pu=1; pu<=npu; pu++)
                          {
                            for(pv=1; pv<=npv; pv++)
                            {
                              for(pw=1; pw<=npw; pw++)
                              {
	val = val + wx[px]*wy[py]*wz[pz]*wt[pt]*wu[pu]*wv[pv]*ww[pw]*
	      galf(xx[px],yy[py],zz[pz],tt[pt],uuu[pu],vv[pv],www[pw]);
		              }
		            }
		          }
	                }
	              }
                    }
                  }
                  fg[s(i,i2,i3,i4,i5,i6,i7)] = val;
		}
	      }
	    }
	  }
	}
      }
    }
    tnow = System.currentTimeMillis();
    System.out.println("kg computed stiffness matrix, see above ");
    System.out.println("fg computed forcing function, RHS,see above, in"+
                       ((tnow-tstart)/1000.0)+" seconds");
    System.out.println("solve kg * uu = fg ");
    new simeq(kg, fg, uu);

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("uu computed Galerkin, Ua analytic, error  ");
    for(i=0; i<nx; i++)
    {
      for(i2=0; i2<ny; i2++)
      {
        for(i3=0; i3<nz; i3++)
        {
          for(i4=0; i4<nt; i4++)
          {
            for(i5=0; i5<nu; i5++)
            {
              for(i6=0; i6<nv; i6++)
              {
                for(i7=0; i7<nw; i7++)
                {
        err = Math.abs(uu[s(i,i2,i3,i4,i5,i6,i7)]-Ua[s(i,i2,i3,i4,i5,i6,i7)]);
        if(err>maxerr) maxerr = err;
        avgerr = avgerr + err;
        System.out.println("uu["+i+","+i2+","+i3+","+i4+","+i5+","+i6+","+i7+
			   "]=");
        System.out.println(+uu[s(i,i2,i3,i4,i5,i6,i7)]+
			   ", Ua="+Ua[s(i,i2,i3,i4,i5,i6,i7)]+
                           ", abserr="+err);
		}
	      }
	    }
	  }
	}
      }
    }
    System.out.println("nx="+nx+", ny="+ny+", nz="+nz+", nt="+nt);
    System.out.println("nu="+nu+", nv="+nv+", nw="+nw);
    System.out.println("npx="+npx+", npy="+npy+", npz="+npz+", npt="+npt);
    System.out.println("npu="+npu+", npv="+npv+", npw="+npw);
    System.out.println("                  maxerr="+maxerr+", avgerr="+
                       avgerr/(double)(nxyztuvw));
    tnow = System.currentTimeMillis();
    System.out.println("total time "+((tnow-tstart)/3600000.0)+ " hours");
  } // end fem_check47h_la constructor

  // indexing function  returne 0 to nxyztuvw-1
  int s(int i, int i2, int i3, int i4, int i5, int i6, int i7)
  {
	// row or column index to kg, index to uu and fg
    return i*ny*nz*nt*nu*nv*nw + i2*nz*nt*nu*nv*nw +
           i3*nt*nu*nv*nw + i4*nu*nv*nw + i5*nv*nw + i6*nw + i7;
  } // end s

  // PDE problem definition functions
    double F(double x, double y, double z, double t, 
             double u, double v, double w)
  {
    return Math.sin(x+y+z+t+u+v+w);
  }

  double uana(double x, double y, double z, double t, 
              double u, double v, double w)
  {
    return Math.sin(x+y+z+t+u+v+w);
  }

  double galk(double x, double y, double z, double t,
                double u, double v, double w)
  {                        // Galerkin k stiffness function for this problem 
    return
    (    L.phipppp(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(u,j5,nu-1,ug)*     L.phi(v,j6,nv-1,vg)*
	 L.phi(w,j7,nw-1,wg)+
	 L.phi(x,j,nx-1,xg)*  L.phipppp(y,j2,ny-1,yg)*     L.phi(z,j3,nz-1,zg)*
	 L.phi(t,j4,nt-1,tg)*     L.phi(u,j5,nu-1,ug)*     L.phi(v,j6,nv-1,vg)*
	 L.phi(w,j7,nw-1,wg)+
	 L.phi(x,j,nx-1,xg)*      L.phi(y,j2,ny-1,yg)* L.phipppp(z,j3,nz-1,zg)*
	 L.phi(t,j4,nt-1,tg)*     L.phi(u,j5,nu-1,ug)*     L.phi(v,j6,nv-1,vg)*
	 L.phi(w,j7,nw-1,wg)+
	 L.phi(x,j,nx-1,xg)*      L.phi(y,j2,ny-1,yg)*     L.phi(z,j3,nz-1,zg)*
	 L.phipppp(t,j4,nt-1,tg)* L.phi(u,j5,nu-1,ug)*     L.phi(v,j6,nv-1,vg)*
	 L.phi(w,j7,nw-1,wg)+
	 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.phipppp(u,j5,nu-1,ug)*     L.phi(v,j6,nv-1,vg)*
	 L.phi(w,j7,nw-1,wg)+
	 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(u,j5,nu-1,ug)* L.phipppp(v,j6,nv-1,vg)*
	 L.phi(w,j7,nw-1,wg)+
	 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(u,j5,nu-1,ug)*     L.phi(v,j6,nv-1,vg)*
	 L.phipppp(w,j7,nw-1,wg)+

     2.0*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)*     L.phi(u,j5,nu-1,ug)*     L.phi(v,j6,nv-1,vg)*
	 L.phi(w,j7,nw-1,wg)+
     2.0*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)*     L.phi(u,j5,nu-1,ug)*     L.phi(v,j6,nv-1,vg)*
	 L.phi(w,j7,nw-1,wg)+
     2.0*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)*     L.phi(u,j5,nu-1,ug)*     L.phi(v,j6,nv-1,vg)*
         L.phi(w,j7,nw-1,wg)+
     2.0*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)*   L.phi(u,j5,nu-1,ug)*     L.phi(v,j6,nv-1,vg)*
         L.phi(w,j7,nw-1,wg)+
     2.0*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.phipp(u,j5,nu-1,ug)*     L.phi(v,j6,nv-1,vg)*
	 L.phi(w,j7,nw-1,wg)+
     2.0*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(u,j5,nu-1,ug)*   L.phipp(v,j6,nv-1,vg)*
	 L.phi(w,j7,nw-1,wg)+
     2.0*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(u,j5,nu-1,ug)*     L.phi(v,j6,nv-1,vg)*
	 L.phipppp(w,j7,nw-1,wg)+

     8.0*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(u,j5,nu-1,ug)*      L.phi(v,j6,nv-1,vg)*
	 L.phi(w,j7,nw-1,wg)                                               )*
        (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)*     L.phi(u,i5,nu-1,ug)*      L.phi(v,i6,nv-1,vg)*
	 L.phi(w,i7,nw-1,wg) );
  } // end galk used by integration 

  double galf(double x, double y, double z, double t,
              double u, double v, double w)
  {                                // Galerkin f function for this problem
    return F(x,y,z,t,u,v,w)*
           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)*L.phi(u,i5,nu-1,ug)*L.phi(v,i6,nv-1,vg)*
           L.phi(w,i7,nw-1,wg);
  } // end galf used by integration 

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