// pde_blivet.java  reads  blivet.inp                                       
//               to get geometry and Dirichlet boundary values           
//
// uses pde_read_ucd.java lsfit.java nuderiv3d.java 
//
// The PDE to solve is Uxx(x,y,z)+Uyy(x,y,z)+Uzz(x,y,z)+(3/25)U(x,y,z)=0 
// the known solution for testing, is U(x,y,z)=sin((x+y+z)/5)            
// the solution is to be computed at:  14 points, 17 DOF                 
// (0.5, 0.5, 0.5)  (2.5, 0.5, 0.5)  (4.5, 0.5, 0.5)                     
// (0.5, 1.5, 0.5)  (2.5, 1.5, 0.5)  (4.5, 1.5, 0.5)                     
// (0.5, 2.5, 0.5)  (2.5, 2.5, 0.5)  (4.5, 2.5, 0.5)                     
// (0.5, 3.5, 0.5)  (2.5, 3.5, 0.5)  (4.5, 3.5, 0.5)                     
// (1.5, 4.5, 0.5)  (3.5, 4.5, 0.5)  refinement is allowed, if needed    
// for each x +/- 0.25, and for each y +/- 0.25, and for each z +/- 0.25 
// creating 14 + 28*28*28 = 21,966 DOF, 14 points above plus             
// 21,952 points (0.25, 0.25, 0.25) to (4.75, 3.75, 0.75)                  
//               
// fortunately, our algorithm solves for one DOF at a time,  
// rather than the typical solving of simultaneous equations.
// Given n boundary vertices, at DOF p, compute cxx[n+1] such that       
// Uxx = cxx[p]*U[p] + sum{i!=p, cxx[i]*Ubound[i]}                       
// similarly for Uyy and Uzz for this PDE                                
// Accumulate A = cxx[p]+cyy[p]+czz[p]+3/25   these are times the unknown U[p]
// Accumulate Y = sum{i!=p, cxx[i]*Ubound[ip[i]]} +                          
//                sum{i!=p, cyy[i]*Ubound[ip[i]]} +                          
//                sum{i!=p, czz[i]*Ubound[ip[i]]} +                          
// Now  A*U[p] + Y = 0  we have just summed the terms of the PDE         
// U[p] = -Y/A          then repeat for all p in DOF                     

// uses pde_read_ucd.java  the class defines all variables and arrays 
// all subscripts start with 1 
// vert struct {double x; double y; double z;} 
// cell struct {int i; int matl; int ctyp; int v[9];} 

class pde_blivet
{
  pde_blivet(String file_name)
  {
      int k, m;
    boolean last = false;
    boolean debug = true;
    int order = 2;
    int npwr;
    double u, ue, err, maxerr;
    int ctypn[]  = { 0, 1, 2, 3, 4, 4, 5, 6, 8}; // UCD constants 
    int ndof = 14; // could be dynamically changed 
    double unk[][] = {
      {0.5, 0.5, 0.5}, {2.5, 0.5, 0.5}, {4.5, 0.5, 0.5},
      {0.5, 1.5, 0.5}, {2.5, 1.5, 0.5}, {4.5, 1.5, 0.5},
      {0.5, 2.5, 0.5}, {2.5, 2.5, 0.5}, {4.5, 2.5, 0.5},
      {0.5, 3.5, 0.5}, {2.5, 3.5, 0.5}, {4.5, 3.5, 0.5},
      {1.5, 4.5, 0.5}, {3.5, 4.5, 0.5}};  
    double xg[]; // allocated later when size is known 
    double yg[];
    double zg[];
    double A;  // later solve for U, A*U=Y 
    double Y;
    double U[] = new double[14]; // DOF that must be computed 
    double ub[];                // boundaries zero based
    double cxx[], cyy[], czz[]; // discrete derivatives coefficients 
    int ip[];

    System.out.println("pde_blivet.java running");
    pde_read_ucd B = new pde_read_ucd(file_name);

    System.out.println("pde_read_ucd returned:");
    if(debug) System.out.println("1 based indexing");

    System.out.println("n_vertices (boundary)= "+B.n_vertices);
    System.out.println("n_cells (boundary faces)= "+B.n_cells);
    System.out.println("n_ndata (Dirichlet at vertices)= "+B.n_ndata);
    System.out.println(" ");

    ub = new double[B.n_vertices];
    for(int i=0; i<B.n_vertices; i++) ub[i] = 0.0;
    if(debug)
    {
      System.out.println("boundary vertices:");
      for(int i=1; i<=B.n_vertices; i++)
      {
        System.out.println(i+" x="+B.vert[i].x+
                           ", y="+B.vert[i].y+
                           ", z="+B.vert[i].z);
      } 
      System.out.println("");

      System.out.println("boundary cells:");
      for(int i=1; i<=B.n_cells; i++)
      {
        System.out.println(i+"  cell="+B.cell[i].i+
                           ", matl="+B.cell[i].matl+
                           ", type="+B.cell[i].ctyp);
        k = ctypn[B.cell[i].ctyp];
        for(int j=1; j<=k; j++) System.out.println(j+" "+B.cell[i].v[j]);
        System.out.println("");
      } 
      System.out.println("");

      if(B.n_ndata==1)
      {
        System.out.println("Dirichlet boundary:");
        System.out.println(B.nodes1);
        for(int i=1; i<=B.n_vertices; i++)
          System.out.println(i+"  "+ B.dirichlet_bound[i]);
      }
      System.out.println("");
    } // end debug print 
    for(int i=1; i<=B.n_vertices; i++)
      ub[i-1] = Math.sin((B.vert[i].x+B.vert[i].y+B.vert[i].z)/5.0);


    System.out.println("check input by fitting");
    lsfit LS = new lsfit(3,3,0); // 3rd power, 3D 
    last = false;
    for(int i=1; i<=B.n_vertices; i++)
    {
      if(i==B.n_vertices) last = true;
      LS.fit3_data(B.vert[i].x, B.vert[i].y, B.vert[i].z,
                   B.dirichlet_bound[i], last);
    }
    // check 
    maxerr = 0.0;
    for(int i=1; i<=B.n_vertices; i++)
    {
      ue = B.dirichlet_bound[i];
      u = LS.eval3_poly(B.vert[i].x, B.vert[i].y, B.vert[i].z);
      err = Math.abs(ue-u);
      maxerr = Math.max(maxerr, err);
    }
    System.out.println("fit3 maxerr="+maxerr);

    maxerr = 0.0;
    for(int i=0; i<ndof; i++)
    {
      ue = Math.sin((unk[i][0]+unk[i][1]+unk[i][2])/5.0); // exact solution 
      u = LS.eval3_poly(unk[i][0],unk[i][1],unk[i][2]); // just raw fit 
      err = Math.abs(ue-u);
      maxerr = Math.max(maxerr, err);
    }
    System.out.println("DOF fit check maxerr="+maxerr);


    // build PDE solver 
    System.out.println("build PDE  A, Y and solve ");
    xg = new double[B.n_vertices+1];
    yg = new double[B.n_vertices+1];
    zg = new double[B.n_vertices+1];
    cxx = new double[B.n_vertices+1];
    cyy = new double[B.n_vertices+1];
    czz = new double[B.n_vertices+1];
    ip = new int[B.n_vertices+1];

    nuderiv3dg N = new nuderiv3dg();
    N.self_test(); // better late than never

    for(int i=0; i<B.n_vertices; i++)
    {
      xg[i] = B.vert[i+1].x; // vert is 1 based indexing 
      yg[i] = B.vert[i+1].y;
      zg[i] = B.vert[i+1].z;
    } // last entry becomes each unk later 

    // build A and Y, compute U[p] = -Y/A 
    for(int i=0; i<ndof; i++)
    {
      xg[B.n_vertices] = unk[i][0];
      yg[B.n_vertices] = unk[i][1];
      zg[B.n_vertices] = unk[i][2];
      System.out.println("unk="+B.n_vertices+", xg="+xg[B.n_vertices]+
                         ", yg="+yg[B.n_vertices]+
                         ", zg="+zg[B.n_vertices]);
      // each case unknown is last
      A = 0.0;
      Y = 0.0;
      if(i==0 && debug)
      {
        for(int j=0; j<=B.n_vertices; j++)
          System.out.println("j="+j+", xg="+xg[j]+
                             ", yg="+yg[j]+", zg="+zg[j]);
        System.out.println("input to nu3dxx ");
        System.out.println("nu3dxx(order="+order+
	                   ", npoint="+(B.n_vertices+1)+
                           ", point="+B.n_vertices);
      }
      // nu3dxx(int order, int npoint, int point,
      //        double x[], double y[], double z[], double c[]); 
      m = N.nu3dxx(order, B.n_vertices+1, B.n_vertices,
	       xg, yg, zg, cxx, ip);    
      for(int j=0; j<=B.n_vertices; j++)
        System.out.println("after xx j="+j+", xg="+xg[j]+
                           ", yg="+yg[j]+", zg="+zg[j]);
      System.out.println("m points are");
      for(int j=0; j<m; j++)
	System.out.println("ip["+j+"]="+ip[j]);
      A += cxx[0]; // unk coef U[ip[0]]
      for(int j=1; j<m; j++) Y += cxx[j]*ub[ip[j]];

      m = N.nu3dyy(order, B.n_vertices+1, B.n_vertices,
	       xg, yg, zg, cyy, ip);    
      for(int j=0; j<=B.n_vertices; j++)
        System.out.println("after yy j="+j+", xg="+xg[j]+
                           ", yg="+yg[j]+", zg="+zg[j]);
      System.out.println("m points are");
      for(int j=0; j<m; j++)
	System.out.println("ip["+j+"]="+ip[j]);
      A += cyy[0]; // unk coef U[ip[0]]
      for(int j=1; j<m; j++) Y += cyy[j]*ub[ip[j]];

      m = N.nu3dzz(order, B.n_vertices+1, B.n_vertices,
	       xg, yg, zg, czz, ip);    
      for(int j=0; j<=B.n_vertices; j++)
        System.out.println("after zz j="+j+", xg="+xg[j]+
                           ", yg="+yg[j]+", zg="+zg[j]);
      System.out.println("m points are");
      for(int j=0; j<m; j++)
	System.out.println("ip["+j+"]="+ip[j]);
      A += czz[0]; // unk coef U[ip[0]]
      for(int j=1; j<m; j++) Y += czz[j]*ub[ip[j]];

      // collect terms of PDE  n_vertices is p
      A += 3.0/25.0;
      U[i] = -Y/A;
    }

    // print solution with checking 
    maxerr = 0.0;
    for(int i=0; i<ndof; i++)
    {
      ue = Math.sin((unk[i][0]+unk[i][1]+unk[i][2])/5.0); // exact solution 
      err = Math.abs(ue-U[i]);
      maxerr = Math.max(maxerr, err);
      System.out.println(i+" x="+unk[i][0]+
                         ", y="+unk[i][1]+
                         ", z="+unk[i][2]+
                         ", U="+U[i]+
                         ", ue="+ue+
                         ", err="+err);
    }
    System.out.println("computed values of DOF maxerr="+maxerr);

    System.out.println("pde_blivet.java ends");
  } // end pde_blivet constructor

  public static void main(String args[])
  {
    // read in name of .inp                                             
    String file_name = "blivet.inp";
    if(args.length>0)
      file_name = args[0];
    System.out.println("input UCD file_name is "+file_name);
    new pde_blivet(file_name);
  }
} // end class pde_blivet