/* pde_blivet.c  reads  blivet.inp                                       */
/*               to get geometry and Dirichlet boundary values           */
/* 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]   these are times the unknown U[p]*/
/* Accumulate Y = sum{i!=p, cxx[i]*Ubound[i]} +                          */
/*                sum{i!=p, cyy[i]*Ubound[i]} +                          */
/*                sum{i!=p, czz[i]*Ubound[i]} +                          */
/*                2/25  the coefficient of U(x,y,z) in the PDE           */
/* 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                     */
/* see code below for details of indexing                                */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#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];} */
#include "lsfit.h"      /* onlu used for testing */
#include "nuderiv3dg.h"  /* for discretization */

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


int main(int argc, char * argv[])
{
  int i, j, k, m, which, last;
  int debug = 1;
  int order = 2; /* not enough independent points for 3 in blivet.inp */
                 /* thus second case of blivet_nr.inp                 */
  int npwr;
  double u, ue, err, maxerr;
  double x, y, z;
  int  ctypn[9] = { 0, 1, 2, 3, 4, 4, 5, 6, 8}; /* UCD constants */
  int ndof = 14; /* could be dynamically changed */
  double unk[14][3] = {
    {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[14]; /* DOF that must be computed */
  double *cxx, *cyy, *czz; /* discrete derivatives coefficients */
  int *ip;
  double *ub; /* boundary, zero index based, recomputed accurately */

  printf("pde_blivet.c running\n");
  which = 0; /* use data as read */
  if(argc>1)
  {
    printf("reading %s\n", argv[1]);
    pde_read_ucd(argv[1]);
    which = 1; /* overwrite boundary */
  }
  else
  {
    printf("reading %s\n", "blivet.inp");
    pde_read_ucd("blivet.inp");
    which = 0; /* built in solution sin((x+y+z)/5) */
  }
  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];
      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++)
  {
    if(which!=0) /* due to builtin checking below */
    {
      ub[i-1] = sin((vert[i].x + vert[i].y + vert[i].z)/5.0);
      printf("ub[%d]=%9.6f\n", i-1, ub[i-1]);
    } 
    else
      ub[i-1] = dirichlet_bound[i];
  } 

  printf("check input by fitting\n");
  fit_init(3,3); /* 3rd power, 3D */
  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-1], last);
  }
  /* check */
  maxerr = 0.0;
  for(i=1; i<=n_vertices; i++)
  {
    ue = ub[i-1];
    u = eval3_poly(vert[i].x, vert[i].y, vert[i].z);
    err = abs(ue-u);
    maxerr = max(maxerr, err);
  }
  printf("fit3 maxerr=%g\n", maxerr);

  maxerr = 0.0;
  for(i=0; i<ndof; i++)
  {
    ue = sin((unk[i][0]+unk[i][1]+unk[i][2])/5.0); /* exact solution */
    u = eval3_poly(unk[i][0],unk[i][1],unk[i][2]); /* just raw fit */
    err = abs(ue-u);
    maxerr = max(maxerr, err);
  }
  printf("DOF fit check maxerr=%g\n", maxerr);


  /* build PDE solver */
  printf("build PDE  A, Y and solve \n");
  xg =  (double *)malloc((n_vertices+1)*sizeof(double));
  yg =  (double *)malloc((n_vertices+1)*sizeof(double));
  zg =  (double *)malloc((n_vertices+1)*sizeof(double));
  cxx = (double *)malloc((n_vertices+1)*sizeof(double));
  cyy = (double *)malloc((n_vertices+1)*sizeof(double));
  czz = (double *)malloc((n_vertices+1)*sizeof(double));
  ip =     (int *)malloc((n_vertices+1)*sizeof(int));
  for(i=0; i<n_vertices; i++)
  {
    xg[i] = vert[i+1].x; /* vert is 1 based indexing */
    yg[i] = vert[i+1].y;
    zg[i] = vert[i+1].z;
  } /* last entry becomes each unk later */

  self_test(); /* nuderiv3dg.c */

  /* 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[n_vertices] = unk[i][0];
    yg[n_vertices] = unk[i][1];
    zg[n_vertices] = unk[i][2];
    if(i==0 && debug)
    {
      printf("input to nu3dxx for DOF %d\n",i);
      printf("m=nu3dxx(order=%d,npoints=%d,point=%d,\n",
	     order, n_vertices+1, n_vertices);
      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 = nu3dxx(int order, int npoint, int point,
       double x[], double y[], double z[], double c[], ip[]); */
    m = nu3dxx(order, n_vertices+1, n_vertices,
	       xg, yg, zg, cxx, ip);    
    A += cxx[0];
    for(j=1; j<m; j++)
    {
      Y += cxx[j]*ub[ip[j]];
    }
    m = nu3dyy(order, n_vertices+1, n_vertices,
           xg, yg, zg, cyy, ip);    
    A += cyy[0];
    for(j=1; j<m; j++)
    {
      Y += cyy[j]*ub[ip[j]];
    }
    m = nu3dzz(order, n_vertices+1, n_vertices,
           xg, yg, zg, czz, ip);    
    A += czz[0];
    for(j=1; j<m; j++)
    {
      Y += czz[j]*ub[ip[j]];
    }
    A += 3.0/25.0;
    U[i] = -Y/A;
  }

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


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