// pde44e_nuderiv4d.java
//                    forth order pde with fourth degree solution
//                    sufficient points must be available for each order,degree
//
// given a PDE, a scatering of points in 4D, boundary values,
//
// solve uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) +
//       utttt(x,y,z,t) + u(x,y,z,t) =  F(x,y,z,t)
// F(x,y,z,t) = 2*(exp(x)+exp(y)+exp(z)+exp(t))
//
// boundary conditions computed using u(x,y,z,t)
// analytic solution is uana(x,y,z,t) = exp(x)+exp(y)+exp(z)+exp(t)
//                                   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 cxx etc. at a subset of vertices
//
//  Uxx(xf,yf,zf,tf)= cxx(xf,yf,zf,tf)*U(xf,yf,zf,tf)+  unknown U m=70 points 
//                    cxx(xfi,yfi,zfi,tfi)*U(xfi,yfi,zfi.tfi)+...+ 
//                             unknown U(xfi,yfi,zfi,tfi) for each xf,yf,zf,tf
//                    cxx(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) - cxx(xbj,ybj,zbj,tbj)*U(xbj,ybj,zbj,tbj) -
//
//                                         ip[1]=j  any bound xbj,ybj,zbj,tbj
//
// same code for a21, a22, ... new N.nu4d4 on next point
// then rest of DOF points
// solve for U(xfi,yfi)  i=1..nfree

class pde44e_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;

  pde44e_nuderiv4d()
  {
    int order = 4;
    int m;
    final int nc=70;   // present limit of nuderiv4d
    double cxxxx[] = new double[nc];
    double cyyyy[] = new double[nc];
    double czzzz[] = new double[nc];
    double ctttt[] = new double[nc];
    int ip[] = new int[nc];
    double ap;
    double error, maxerror;

    t_start = System.currentTimeMillis();

    System.out.println("pde44e_nuderiv4d.java running");
    System.out.println("  uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) +");
    System.out.println("  utttt(x,y,z,t) + u(x,y,z,t) = F(x,y,z,t)");
    System.out.println("F(x,y,z,t) = 2*(exp(x)+exp(y)+exp(z)+exp(t))");
    System.out.println("uana(x,y)=exp(x)+exp(y)+exp(z)+exp(t)");
    System.out.println("                   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
    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]);
    }

    // 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 
    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]);
    }

    maxerror = 0.0;

    for(int point=0; point<nfree; point++) // for DOF 
    {
      m = nvert;
      N.nu4dxxxx(order, nvert, point, x, y, z, t, cxxxx);
      N.nu4dyyyy(order, nvert, point, x, y, z, t, cyyyy);
      N.nu4dzzzz(order, nvert, point, x, y, z, t, czzzz);
      N.nu4dtttt(order, nvert, point, x, y, z, t, ctttt);

      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(cxxxx[point], cyyyy[point], czzzz[point],
                           ctttt[point]);
      for(int i=0; i<m; i++)
      {
	if(i==point) continue;
        ap=aij(cxxxx[i], cyyyy[i], czzzz[i], ctttt[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 cxxxx, cyyyy,  ...
    } // 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 pde44e_nuderiv4d.java");
  }

  // PDE definition functions
  double f(double x, double y, double z, double t)
  {
      return 2.0*(Math.exp(x)+Math.exp(y)+Math.exp(z)+Math.exp(t));

  }

  // boundary and test function
    double uana(double x, double y, double z, double t)
  {
    return Math.exp(x)+Math.exp(y)+Math.exp(z)+Math.exp(t);
  }

  // PDE coefficients
    double aii(double axxxx, double ayyyy, double azzzz, double atttt)
  {
    // solve uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) +
    //       utttt(x,y,z,t) + u(x,y,z,t) =

    return axxxx + ayyyy + azzzz + atttt + 1.0;
  }
  // PDE derivatives
  double aij(double axxxx, double ayyyy, double azzzz, double atttt)
  {
    // solve uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) +
    //       utttt(x,y,z,t) +    no U term

    return    axxxx +  ayyyy + azzzz + atttt;
  }

  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 pde44e_nuderiv4d();
  }
} // end  class pde44e_nuderiv4d