// Pde34_nuderiv4d.java
//                    third order pde with third degree solution, 4D
//                    sufficient points must be available for each order,degree
//
// given a PDE, a scatering of points in 4D, boundary values,
//
// solve uxxx(x,y,z,t) + 2 uyyy(x,y,z,t) + 3 uzzz(x,y,z,t) +
//     4 uttt(x,y,z,t) + 5 uxyt(x,y,z,t) + 6 uyyz(x,y,z,t) +
//     7 uztt(x,y,z,t) + 8 uyzt(x,y,z,t) + 9 uyyt(x,y,z,t) +
//     10  ut(x,y,z,t) + 11   u(x,y,z,t) =      F(x,y,z,t)
//     F(x,y,z,t) =  608 + 30*t*t   + 60*y*y   + 70*x*y   + 80*y*z   +
//              11*t*t*t + 22*z*z*z + 33*y*y*y + 44*x*x*x + 55*y*y*z +
//              66*y*y*t + 77*x*y*t + 88*y*z*t + 99*x*z   + 110*t
//
//
// boundary conditions computed using u(x,y,z,t)
// analytic solution is uana(x,y,z,t) = t*t*t + 2*z*z*z + 3*y*y*y + 4*x*x*x +
//                                    5*y*y*z + 6*y*y*t + 7*x*y*t + 8*y*z*t +
//                                    9*x*z   +10*t     +11; 
//                                   soultion, generally not known
//                                   used to compute boundary values
//                                   used for error checking
// U(xf,yf,zf,tf) unknown at  nfree  vertices, DOF
// U(xb,yb,zb,tb) known at  nbound  boundary vertices, nvert  total vertices
// nuderiv4d determines cxxx, cyyy, czzz, cttt, ...
//                                                   at a subset of vertices
//  Uxxx(xf,yf,zf,tf)= cxxx(xf,yf,zf,tf)*U(xf,yf,zf,tf)+  unknown U
//                     cxxx(xfi,yfi,zfi,tfi)*U(xfi,yfi,zfi.tfi)+...+ 
//                             unknown U(xfi,yfi,zfi,tfi) for each xf,yf,zf,tf
//                    cxxx(xbj,ybj,zbj,tbj)*U(xbj,ybj,zbj,tbj)+...
//                                                              known value
//  ...
//  f(xf,yf,zf,tf)    right hand side, known, computable
//
// substituting in the PDE, a set of linear simultaneous equations for U(xf,yf)
//  | a11  a12  ...  a1n |   |U(xf1,yf1,zf1,tf1)| = rhs(1)  n = nfree
//  | a21  a22  ...  a2n | * |U(xf2,yf2,zf2,tf2)| = rhs(2)
//              ...
//  | an1  an2  ...  ann |   |U(xfn,yfn,zfn,tfn)| = rhs(n)
//
//
// rhs(1) = f(xf1,yf1,zf1,tf1) - cxxx(xbj,ybj,zbj,tbj)*U(xbj,ybj,zbj,tbj) -
//
// same code for a21, a22, ... new N.nuderiv4d on next point
// then rest of DOF points
// solve for U(xfi,yfi)  i=1..nfree

public strictfp class pde34_nuderiv4d
{
  int debug = 1;
  int nvert = 70;                 // number of unique vertices 
  double x[] = new double[nvert];
  double y[] = new double[nvert];
  double z[] = new double[nvert];
  double t[] = new double[nvert];
  int uniqueb[] = new int[nvert];  // boundary vertices 
  int nfree = 1;                // degrees of freedom, must be first vertices 
  int nbound = nvert-nfree;        // number of boundary vertices 
  double t_start, t_init, t_kf, t_simeq, t_end;

  pde34_nuderiv4d()
  {
    int order = 3;
    int m;
    final int nc=70;
    double ct[] = new double[nc];
    double cxxx[] = new double[nc];
    double cxyt[] = new double[nc];
    double cyyy[] = new double[nc];
    double cyyz[] = new double[nc];
    double cyyt[] = new double[nc];
    double cyzt[] = new double[nc];
    double czzz[] = new double[nc];
    double cztt[] = new double[nc];
    double cttt[] = new double[nc];
    int ip[] = new int[nc];
    double ap;
    double error, maxerror;
    double xmin, xmax, ymin, ymax, zmin, zmax, tmin, tmax, umin, umax;

    t_start = System.currentTimeMillis();

    System.out.println("pde34_nuderiv4d.java running");
    System.out.println("  uxxx(x,y,z,t) + 2 uyyy(x,y,z,t) + 3 uzzz(x,y,z,t) +");
    System.out.println("4 uttt(x,y,z,t) + 5 uxyt(x,y,z,t) + 6 uyyz(x,y,z,t) +");
    System.out.println("7 uztt(x,y,z,t) + 8 uyzt(x,y,z,t) + 9  uyyt(x,y,z,t) +");
    System.out.println("10  ut(x,y,z,t) + 11   u(x,y,z,t) = F(x,y,z,t)");
    System.out.println("F(x,y,z,t) = ");
    System.out.println(" 608 + 30*t*t   + 60*y*y   + 70*x*y   + 80*y*z   + ");
    System.out.println("11*t*t*t + 22*z*z*z + 33*y*y*y + 44*x*x*x + 55*y*y*z+");
    System.out.println("66*y*y*t + 77*x*y*t + 88*y*z*t + 99*x*z   + 110*t ");
    System.out.println("             ");
    System.out.println("uana(x,y)=t^3 + 2 z^3 + 3 y^3 + 4 x^3 +");
    System.out.println("          5 y^2 z + 6 y y t + 7 x y t + 8 y z t +");
    System.out.println("          9 x z 10 t + 11   boundary and solution");

    nuderiv4d N = new nuderiv4d();

    double U[] = new double[nvert];          // known boundary 
    double A[][] = new double[nfree][nfree]; // matrix 
    double rhs[] = new double[nfree];        // right hand side, forcing function 
    double ug[] = new double[nfree];         // solution computed 

    // read in or generate data points, DOF (free) first, then boundary
    gen_data4(1);

    // initialize vertices
    xmin = x[0]; xmax = x[0];
    ymin = y[0]; ymax = y[0];
    zmin = z[0]; zmax = z[0];
    tmin = t[0]; tmax = t[0];
    System.out.println("vertex, x, y");
    for(int i=0; i<nvert; i++)
    {
      U[i] = 0.0;  // known at boundaries
      System.out.println(i+"  x="+x[i]+", y="+y[i]+", z="+z[i]+", t="+t[i]);
      xmin = Math.min(xmin,x[i]);
      xmax = Math.max(xmax,x[i]);
      ymin = Math.min(ymin,y[i]);
      ymax = Math.max(ymax,y[i]);
      zmin = Math.min(zmin,z[i]);
      zmax = Math.max(zmax,z[i]);
      tmin = Math.min(tmin,t[i]);
      tmax = Math.max(tmax,t[i]);
    }

    // clear A and rhs
    for(int i=0; i<nfree; i++)
    {
      rhs[i] = 0.0;
      for(int k=0; k<nfree; k++) A[i][k] = 0.0;
    }

    // set boundary values 
    umin = uana(x[0],y[0],z[0],t[0]);
    umax = umin;
    for(int i=nfree; i<nvert; i++)
    {
      U[i]=uana(x[i],y[i],z[i],t[i]);
      System.out.println("bound "+i+"  U[i]="+U[i]);
      umin = Math.min(umin,U[i]);
      umax = Math.max(umax,U[i]);
    }

    System.out.println("xmin="+xmin+",  xmax="+xmax);
    System.out.println("ymin="+ymin+",  ymax="+ymax);
    System.out.println("zmin="+zmin+",  zmax="+zmax);
    System.out.println("tmin="+tmin+",  tmax="+tmax);
    System.out.println("umin="+umin+",  umax="+umax);
    maxerror = 0.0;

    for(int point=0; point<nfree; point++) // for DOF 
    {
      m = nvert;
      N.nu4dxxx(order, nvert, point, x, y, z, t, cxxx);
      if(point==0 && debug>0)
      {
	for(int k=0; k<nvert; k++)
	  System.out.println("cxxx["+k+"]="+cxxx[k]);
      }
      N.nu4dyyy(order, nvert, point, x, y, z, t, cyyy);
      N.nu4dzzz(order, nvert, point, x, y, z, t, czzz);
      N.nu4dttt(order, nvert, point, x, y, z, t, cttt);
      N.nu4dxyt(order, nvert, point, x, y, z, t, cxyt);
      N.nu4dyyz(order, nvert, point, x, y, z, t, cyyz);
      N.nu4dztt(order, nvert, point, x, y, z, t, cztt);
      N.nu4dyzt(order, nvert, point, x, y, z, t, cyzt);
      N.nu4dyyt(order, nvert, point, x, y, z, t, cyyt);
      N.nu4dt  (order, nvert, point, x, y, z, t, ct);

      System.out.println("order="+order+", m="+m+", point="+point);
      // build row for point
      rhs[point] = f(x[point],y[point],z[point],t[point]);
      A[point][point]+=aii(cxxx[point], cyyy[point], czzz[point],
                           cttt[point], cxyt[point], cyyz[point],
                           cztt[point], cyzt[point], cyyt[point],
                             ct[point]);
      for(int i=0; i<m; i++)
      {
	if(i==point) continue;
        ap=aij(cxxx[i], cyyy[i], czzz[i], cttt[i],
               cxyt[i], cyyz[i], cztt[i], cyzt[i],
               cyyt[i],   ct[i]);

        if(i<nfree)
	{
	  A[point][i] += ap;
          if(debug>0) 
            System.out.println("point="+point+", i="+i+", ap="+ap);
        }
	else // boundary
	{
	  rhs[point] -= ap*U[i];
          if(debug>0)
            System.out.println("point="+point+", i="+i+", ap*U="+
	 		       (-ap*U[i]));
	}
      } // end i across cxxx, cyyy,  ...
    } // end point
    if(debug>0)
    {
      System.out.println("A matrix and rhs");
      for(int i=0; i<nfree; i++)
      {
        for(int k=0; k<nfree; k++)
	  System.out.println("A["+i+"]["+k+"]="+A[i][k]);
        System.out.println("rhs["+i+"]="+rhs[i]);
      }
      System.out.println(" ");
    }
    // solve linear equations
    new simeq(A, rhs, ug);
    for(int i=0; i<nfree; i++)
    {
      error = Math.abs(ug[i]-uana(x[i],y[i],z[i],t[i]));
      maxerror = Math.max(maxerror, error);
      System.out.println("vertex "+i+" x="+x[i]+", y="+y[i]+
                         ", z="+z[i]+", t="+t[i]);
      System.out.println("    computed="+ug[i]+", actual="+
                         uana(x[i],y[i],z[i],t[i])+", error = "+error);
    }
    System.out.println("maxerror = "+maxerror);
    System.out.println(" ");
    System.out.println("end pde34_nuderiv4d.java");
  }

  // PDE definition functions
  double f(double x, double y, double z, double t)
  {
    return 0.0;

  }

  // boundary and test function
    double uana(double x, double y, double z, double t)
  {
    return t*t*t + 2.0*z*z*z + 3.0*y*y*y + 4.0*x*x*x* + 5.0*y*y*z +
       6.0*y*y*t + 7.0*x*y*t + 8.0*y*z*t + 9.0*x*z    + 10.0*t    + 11.0;
  }

  // PDE coefficients
    double aii(double uxxx, double uyyy, double uzzz,
               double uttt, double uxyt, double uyyz,
               double uztt, double uyzt, double uyyt,
               double ut)
  {
    // solve uxxx(x,y,z,t) + 2 uyyy(x,y,z,t) + 3 uzzz(x,y,z,t) +
    //     4 uttt(x,y,z,t) + 5 uxyt(x,y,z,t) + 6 uyyz(x,y,z,t) +
    //     7 uztt(x,y,z,t) + 8 uyzt(x,y,z,t) + 9 uyyt(x,y,z,t) +
    //    10   ut(x,y,z,t) + 11   u(x,y,z,t) = f(x,y,z,t)

    return    uxxx +  2.0*uyyy +  3.0*uzzz + 
          4.0*uttt +  5.0*uxyt +  6.0*uyyz +
          7.0*uztt +  8.0*uyzt +  9.0*uyyt  +
         10.0*ut   + 11.0;
  }
  // PDE derivatives
  double aij(double uxxx, double uyyy, double uzzz,
             double uttt, double uxyt, double uyyz,
             double uztt, double uyzt, double uyyt,
             double ut)
  {
    // solve uxxx(x,y,z,t) + 2 uyyy(x,y,z,t) + 3 uzzz(x,y,z,t) +
    //     4 uttt(x,y,z,t) + 5 uxyt(x,y,z,t) + 6 uyyz(x,y,z,t) +
    //     7 uztt(x,y,z,t) + 8 uyzt(x,y,z,t) + 9 uyyt(x,y,z,t) +
    //    10   ut(x,y,z,t)      no U term

    return    uxxx +  2.0*uyyy + 3.0*uzzz + 
          4.0*uttt +  5.0*uxyt + 6.0*uyyz +
          7.0*uztt +  8.0*uyzt + 9.0*uyyt  +
         10.0*ut;
  }

  void gen_data4(int which)
  {
    double rx = 0.25;
    double ry = 0.28;
    double rz = 0.33;
    double rt = 0.37;
    double drx = 0.0096;
    double dry = 0.0094;
    double drz = 0.0091;
    double drt = 0.0089;
    double thx = 0.2;
    double thy = 1.7;
    double thz = 2.3;
    double tht = 3.32;
    double dthx = 0.61;
    double dthy = 0.72;
    double dthz = 0.49;
    double dtht = 0.53;
    double d5th = 0.628318531;
    double d3rd = 1.047197551;
    double th5th = 1.256637061; 
    double r = 0.2;
    double theta1 = 0.05;
    double phi1 = 0.1;
    double psi1 = 0.15;
    int    nn = 0;
    double theta, phi, psi;

    if(which==0)
    { 
      for(int i=0; i<nvert; i++)
      {
        x[i] = rx*Math.sin(thx);
        y[i] = ry*Math.cos(thy);
        z[i] = rz*Math.sin(thz);
        t[i] = rt*Math.cos(tht);
        rx = rx + drx;
        ry = ry + dry;
        rz = rz + drz;
        rt = rt + drt;
        thx = thx + dthx;
        thy = thy + dthy;
        thz = thz + dthz;
        tht = tht + dtht;
      }
    }
    else
    {
      nn = 0;
      psi = psi1;
      for(int ipsi=0; ipsi<3; ipsi++)
      {
	phi = phi1;
	for(int iphi=0; iphi<5; iphi++)
	{
	  theta = theta1;
	  for(int itheta=0; itheta<5; itheta++)
	  {
	    x[nn] = r*Math.sin(psi)*Math.sin(phi)*Math.cos(theta); 
            y[nn] = r*Math.sin(psi)*Math.sin(phi)*Math.sin(theta);
            z[nn] = r*Math.sin(psi)*Math.cos(phi);
            t[nn] = r*Math.cos(psi);
            nn++;
            if(nn==nvert) return;
            r += 0.011;
	    theta += th5th;
          }
          theta1 += 0.008;
          phi += d5th;
        }
        phi1 += 0.009;
        psi += d3rd;
      }
    } // end if 
  } // end gen_data4

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