// pde_spike.c  reads  spike.inp                                         
//              to get geometry and Dirichlet boundary values           
// The PDE to solve is 
// Uxxxx(x,y,z) + Uyyyy(x,y,z) + Uzzzz(x,y,z)+
// 2*Uxx(x,y,z) + 2*Uyy(x,y,z) + 2*Uzz(x,y,z) + 3*U(x,y,z) = 0 
// the known solution for testing, is U(x,y,z)=sin(x+y+z+1)            
// the solution is to be computed at points inside               
// for each x step 0.25, and for each y step 0.25, and for each z step 0.25
//                                                                       
// fortunately, our algorithm solves for one DOF at a time,              
// rather than the typical solving of simultaneous equation for all points   
// 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 = cxxxx[p], cxx[p], ...   these are times the unknown U[p]
// Accumulate Y = sum{i!=p, cxxxx[i]*Ubound[i]} +                          
//                sum{i!=p, cxx[i]*Ubound[i]} +                          
//                sum{i!=p, ...[i]*Ubound[i]} +                          
//           
// Now  A*U[p] + Y = 0  we have just summed the terms of the PDE         
// solve for U[p] A * U[p] = -Y                     
// see code below for details of indexing                                

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "lsfit.h"       // check boundary
#include "nuderiv3dg.h"  // for discretization number of known
                         // points importamnt 

#define _MAIN_
#include "pde_read_ucd.h" // 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];} 

#undef abs
#define abs(a) ((a)<0.0?(-(a)):(a))
#undef max
#define max(a,b) ((a)>(b)?(a):(b))
#undef min
#define min(a,b) ((a)<(b)?(a):(b))


int main(int argc, char * argv[])
{
  int i, j, k, m, last, c;
  int debug = 0;
  int order = 4; // need enough independent points for 3 in spike.inp
  int npwr = 4;
  double u, ue, err, maxerr, avgerr, rmserr;
  double x, y, z;
  int  ctypn[9] = {0, 1, 2, 3, 4, 4, 5, 6, 8}; // UCD constants 
  int ndof = 0; // dynamically changed 
  double unk[500][3];
  double *xg; // allocated later when size is known 
  double *yg;
  double *zg;
  int point = 0; // xg index, unknown vertex for solving
  double A;  // later solve for U, A*U=Y 
  double Y;
  double *U; // DOF that must be computed 
  double xmin, xmax, hx, ymin, ymax, hy, zmin, zmax, hz;
  double *cxxxx, *cxx, *cyyyy, *cyy, *czzzz, *czz; // discrete derivatives
                                                   // coefficients 
  int *ip;
  double *ub; // boundary, zero index based 

  printf("pde_spike.c running\n");
  if(argc>1)
  {
    printf("reading %s\n", argv[1]);
    pde_read_ucd(argv[1]);
  }
  else
  {
    printf("reading %s\n", "spike.inp");
    pde_read_ucd("spike.inp");
  }
  printf("pde_read_ucd returned:\n");
  if(debug) printf("1 based indexing\n");

  printf("n_vertices (boundary)= %d \n",n_vertices);
  printf("n_cells (boundary faces)= %d \n",n_cells);
  printf("n_ndata (Dirichlet at vertices)= %d \n",n_ndata);
  printf("\n");

  if(debug)
  {
    printf("boundary vertices:\n");
    for(i=1; i<=n_vertices; i++)
    {
      printf("%d  x=%f, y=%f, z=%f \n",i, vert[i].x, vert[i].y, vert[i].z);
    } 
    printf("\n");

    printf("boundary cells:\n");
    for(i=1; i<=n_cells; i++)
    {
      printf("%d  matl=%d, type=%d ", cell[i].i, cell[i].matl, cell[i].ctyp);
      k = ctypn[cell[i].ctyp];
      if(debug)
      {
        for(j=1; j<=k; j++) 
          printf("%d ", cell[i].v[j]);
        printf("\n");
      }
    } 
    printf("\n");

    if(n_ndata==1)
    {
      printf("Dirichlet boundary:\n");
      //  printf("%s", nodes1);
      for(i=1; i<=n_vertices; i++)
      {
        printf("%d  %f\n", i, dirichlet_bound[i]);
      } 
      printf("\n");
    }
    printf("\n");
  } // end debug print 

  ub =  (double *)malloc((n_vertices+1)*sizeof(double));
  for(i=1; i<=n_vertices; i++)
  {
    ub[i] = dirichlet_bound[i];
  } 

  printf("check input by fitting 3 variables with 4th power\n");
  fit_init(3,4); // 3D, 4th power not enough points for 5Th 
  last = 0;
  for(i=1; i<=n_vertices; i++)
  {
    if(i==n_vertices) last = 1;
    fit3_data(vert[i].x, vert[i].y, vert[i].z, ub[i], last);
  }
  // check accuracy of fit (bad means maybe non smooth
  printf("check boundary fit\n");
  maxerr = 0.0;
  for(i=1; i<=n_vertices; i++)
  {
    ue = ub[i];
    u = eval3_poly(vert[i].x, vert[i].y, vert[i].z);
    err = abs(ue-u);
    maxerr = max(maxerr, err);
  }
  printf("boundary fit3 maxerr=%g\n", maxerr);


  // build PDE solver 
  printf("build PDE  A, Y and solve \n");
  xg =    (double *)malloc((n_vertices+2)*sizeof(double));
  yg =    (double *)malloc((n_vertices+2)*sizeof(double));
  zg =    (double *)malloc((n_vertices+2)*sizeof(double));
  cxxxx = (double *)malloc((n_vertices+2)*sizeof(double));
  cyyyy = (double *)malloc((n_vertices+2)*sizeof(double));
  czzzz = (double *)malloc((n_vertices+2)*sizeof(double));
  cxx =   (double *)malloc((n_vertices+2)*sizeof(double));
  cyy =   (double *)malloc((n_vertices+2)*sizeof(double));
  czz =   (double *)malloc((n_vertices+2)*sizeof(double));
  ip =    (int *)   malloc((n_vertices+2)*sizeof(int));

  for(i=1; i<=n_vertices; i++)
  {
    xg[i] = vert[i].x; // vert is 1 based indexing 
    yg[i] = vert[i].y; // [0] will be unknown
    zg[i] = vert[i].z;
  } // [0] entry becomes each unk later 
  printf("built xg, yg, zg \n");

  if(pde_closed_ucd()==0)
  {
    printf("UCD polyhedra not closed or repeated points, bye.\n");
    return 0;
  }
  // find points inside volume
  // first get limits
  xmin = vert[1].x;
  ymin = vert[1].y;
  zmin = vert[1].z;
  xmax = vert[1].x;
  ymax = vert[1].y;
  zmax = vert[1].z;
  for(i=2; i<=n_vertices; i++)
  {
    xmin = min(xmin, vert[i].x);
    ymin = min(ymin, vert[i].y);
    zmin = min(zmin, vert[i].z);
    xmax = max(xmax, vert[i].x);
    ymax = max(ymax, vert[i].y);
    zmax = max(zmax, vert[i].z);
  }
  hx = (xmax-xmin)/(double)8;
  hy = (ymax-ymin)/(double)8;
  hz = (zmax-zmin)/(double)8;
  printf("xmin=%f, xmax=%f, hx=%f \n", xmin, xmax, hx);
  printf("ymin=%f, xmax=%f, hy=%f \n", ymin, ymax, hy);
  printf("zmin=%f, xmax=%f, hz=%f \n", zmin, zmax, hz);

  for(x=xmin+hx/2.0; x<xmax-hx/2.0; x+=hx)
  {
    for(y=ymin+hy/2.0; y<ymax-hy/2.0; y+=hy)
    {
      for(z=zmin+hz/2.0; z<zmax-hz/2.0; z+=hz)
      {
        c = pde_inside3_ucd(x, y, z);
        if(debug) printf("inside c=%d, at %f %f %f\n", c, x, y, z);
        if(c==1)
	{
          unk[ndof][0] = x;
          unk[ndof][1] = y;
          unk[ndof][2] = z;
          ndof++;
	}
      } 
    } 
  }
  // solve at these internal points
  U =   (double *)malloc((ndof+2)*sizeof(double));
  printf("ndof=%d, vertices for solution \n", ndof);         
  self_test(); // nuderiv3dg.c 
  printf("passed nuderiv3dg.c self_test \n");

  // ndof = 1; // for testing
  // unk[0][0] = 0.0;
  // unk[0][1] = 0.0;
  // unk[0][2] = 0.0;
  // printf("test point = %f, %f, %f \n", unk[0][0], unk[0][1], unk[0][2]);

  // build A and Y, compute U[p] = -Y/A 
  for(i=0; i<ndof; i++)
  {
    A = 0.0; // do one DOF at a time, rather than simeq 
    Y = 0.0;
    xg[point] = unk[i][0];
    yg[point] = unk[i][1];
    zg[point] = unk[i][2];

    order = 4;

    if(i==0 && debug)
    {
      // m = nu3dxxxx(int order, int npoint, int point,
      //   double xg[], double yg[], double zg[], double c[], ip[]); 
      printf("m=nu3dxxxx(order=%d,npoints=%d,point=%d,\n",
	     order, n_vertices+1, point);
      for(j=0; j<=n_vertices; j++)
        printf("j=%d, xg=%f, yg=%f, zg=%f \n",j,xg[j],yg[j],zg[j]);
    }
    m = nu3dxxxx(order, n_vertices+1, point, xg, yg, zg, cxxxx, ip);    
    A += cxxxx[point]; // the unknown
    for(j=1; j<m; j++)
    {
      if(i==0 && debug) printf("ip[j=%d]=%d,cxxxx[j]=%f, ub[ip[j]]=%f \n",
                               j, ip[j], cxxxx[j], ub[ip[j]]);
      Y += cxxxx[j]*ub[ip[j]];
    }
    if(i==0 && debug) printf("just cxxxxx A=%f, Y=%f \n", A, Y);
   
    if(i==0 && debug)
    {
      printf("m=nu3dyyyy(order=%d,npoints=%d,point=%d,\n",
	     order, n_vertices+1, point);
    }
    m = nu3dyyyy(order, n_vertices+1, point, xg, yg, zg, cyyyy, ip);    
    A += cyyyy[point];
    for(j=1; j<m; j++)
    {
      Y += cyyyy[j]*ub[ip[j]];
    }

    if(i==0 && debug)
    {
      printf("m=nu3dzzzz(order=%d,npoints=%d,point=%d,\n",
	     order, n_vertices+1, point);
    }
    m = nu3dzzzz(order, n_vertices+1, point, xg, yg, zg, czzzz, ip);    
    A += czzzz[point];
    for(j=1; j<m; j++)
    {
      Y += czzzz[j]*ub[ip[j]];
    }

    if(i==0 && debug)
    {
      printf("m=nu3dxx(order=%d,npoints=%d,point=%d,\n",
	     order, n_vertices+1, point);
    }
    m = nu3dxx(order, n_vertices+1, point, xg, yg, zg, cxx, ip);    
    A += 2.0*cxx[point];
    for(j=1; j<m; j++)
    {
      Y += 2.0*cxx[j]*ub[ip[j]];
    }

    if(i==0 && debug)
    {
      printf("m=nu3dyy(order=%d,npoints=%d,point=%d,\n",
	     order, n_vertices+1, point);
    }
    m = nu3dyy(order, n_vertices+1, point, xg, yg, zg, cyy, ip);    
    A += 2.0*cyy[point];
    for(j=1; j<m; j++)
    {
      Y += 2.0*cyy[j]*ub[ip[j]];
    }

    if(i==0 && debug)
    {
      printf("m=nu3dzz(order=%d,npoints=%d,point=%d,\n",
	     order, n_vertices+1, point);
    }
    m = nu3dzz(order, n_vertices+1, point, xg, yg, zg, czz, ip);    
    A += 2.0*czz[point];
    for(j=1; j<m; j++)
    {
      Y += 2.0*czz[j]*ub[ip[j]];
    }
    A += 3.0; //  coef of U() term in PDE
    if(debug) printf("ready to solve A=%e. Y=%e \n", A, Y);
    U[i] = -Y/A;
  } // end i on ndof
  printf("finished ndof loop \n");

  // print solution with checking 
  printf("checking solution at %d points\n", ndof);
  maxerr = 0.0;
  avgerr = 0.0;
  rmserr = 0.0;
  for(i=0; i<ndof; i++)
  {
    ue = sin(unk[i][0]+unk[i][1]+unk[i][2]+1.0); // exact solution 
    err = abs(ue-U[i]);
    maxerr = max(maxerr, err);
    avgerr = avgerr + err;
    rmserr = rmserr + err*err;
    printf("%d x=%7.4f, y=%7.4f, z=%7.4f , U=%7.4f, ue=%f, err=%g\n",
           i, unk[i][0], unk[i][1], unk[i][2], U[i], ue, err);
  }
  avgerr = avgerr/(double)ndof;
  rmserr = sqrt(rmserr/(double)ndof);
  printf("computed values of DOF maxerr=%g\n", maxerr);
  printf("average error=%g, RMS error=%g\n", avgerr, rmserr);

  printf("pde_spike.c ends\n");
  return 0;
} // end pde_spike.c