/* nag_pde_parab_1d_cd_ode_remesh (d03psc) Example Program. 
 *
 * Copyright 2001 Numerical Algorithms Group.
 *
 * Mark 7, 2001. 
 */

#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd03.h>
#include <nagx01.h>

int ex1(void), ex2(void);

static void uvin1(Integer, Integer, Integer, const double[],
                  const double[], double[],Integer, double[], 
                  Nag_Comm *);

static void uvin2(Integer, Integer, Integer, const double[],
                  const double[], double[],Integer, double[], 
                  Nag_Comm *);

static void pdef1(Integer, double, double, const double[],
                  const double[], Integer, const double[],
                  const double[], double[], double[], double[],
                  double[], Integer *, Nag_Comm *);

static void pdef2(Integer, double, double, const double[],
                  const double[], Integer, const double[],
                  const double[], double[], double[], double[],
                  double[], Integer *, Nag_Comm *);

static void bndry1(Integer, Integer, double, const double[],
                   const double[], Integer, const double[],
                   const double[], Integer, double[],
                   Integer *, Nag_Comm *);

static void bndry2(Integer, Integer, double, const double[],
                   const double[], Integer, const double[],
                   const double[], Integer, double[],
                   Integer *, Nag_Comm *);

static void monit1(double, Integer, Integer, const double[], 
                   const double[], double[], Nag_Comm *);

static void monit2(double, Integer, Integer, const double[], 
                   const double[], double[], Nag_Comm *);


static void nmflx1(Integer, double, double, Integer, 
                   const double[], const double[],
                   const double[], double[], Integer *,
                   Nag_Comm *, Nag_D03_Save *);

static void nmflx2(Integer, double, double, Integer, 
                   const double[], const double[],
                   const double[], double[], Integer *,
                   Nag_Comm *, Nag_D03_Save *);

static void exact(double, double *, const double *, Integer, 
                  Integer);

#define P(I,J) p[npde*((J)-1)+(I)-1]
#define UE(I,J) ue[npde*((J)-1)+(I)-1]
#define U(I,J) u[npde*((J)-1)+(I)-1]
#define UOUT(I,J,K) uout[npde*(intpts*((K)-1)+(J)-1)+(I)-1]

int main(void)
{
  Vprintf("d03psc Example Program Results\n");
  ex1();
  ex2();
  return 0;
}


int ex1(void)
{
  const Integer npde=1, npts=61, ncode=0, nxi=0, nxfix=0, itype=1, 
    neqn=npde*npts+ncode, intpts=7, lisave=25+nxfix+neqn, 
    nwkres=npde*(3*npts+3*npde+32)+7*npts+3, lenode=11*neqn+50,
    mlu=3*npde-1, lrsave=(3*mlu+1)*neqn+nwkres+lenode;
  static double xout[7] = { .2,.3,.4,.5,.6,.7,.8 };
  double con, dxmesh, tout, trmesh, ts, xratio;
  Integer exit_status, i, ind, ipminf, it, itask, itol, 
    itrace, m, nrmesh;
  Boolean remesh;
  double *algopt=0, *atol=0, *rsave=0, *rtol=0, 
    *u=0, *uout=0, *x=0, *xfix=0, *xi=0;
  Integer *isave=0;
  NagError fail;
  Nag_Comm comm;
  Nag_D03_Save saved;

  INIT_FAIL(fail);
  exit_status = 0;

  /* Allocate memory */

  if ( !(algopt = NAG_ALLOC(30, double)) ||
       !(atol = NAG_ALLOC(1, double)) ||
       !(rsave = NAG_ALLOC(lrsave, double)) ||
       !(rtol = NAG_ALLOC(1, double)) ||
       !(u = NAG_ALLOC(npde*npts, double)) ||
       !(uout = NAG_ALLOC(npde*intpts*itype, double)) ||
       !(x = NAG_ALLOC(npts, double)) ||
       !(xfix = NAG_ALLOC(1, double)) ||
       !(xi = NAG_ALLOC(1, double)) ||
       !(isave = NAG_ALLOC(lisave, Integer)) )
    {
      Vprintf("Allocation failure\n");
      exit_status = 1;
      goto END;
    }
    
  Vprintf("\n\nExample 1\n\n");

  itrace = 0;
  itol = 1;
  atol[0] = 1.0e-4;
  rtol[0] = 1.0e-4;
    
  Vprintf("  npts = %4ld", npts);
  Vprintf(" atol = %10.3e", atol[0]);
  Vprintf(" rtol = %10.3e\n\n", rtol[0]);

  /* Initialise mesh */

  for (i = 0; i < npts; ++i) x[i] = i/(npts-1.0);
  xfix[0] = 0.0;

  /* Set remesh parameters */

  remesh = 1;
  nrmesh = 3;
  dxmesh = 0.0;
  trmesh = 0.0;
  con = 2.0/(npts-1.0);
  xratio = 1.5;
  ipminf = 0;

  xi[0] = 0.0;
  ind = 0;
  itask = 1;

  for (i = 0; i < 30; ++i) algopt[i] = 0.0;

  /* b.d.f. integration */

  algopt[0] = 1.0;
  algopt[12] = 0.005;

  /* Loop over output value of t */

  ts = 0.0;
  tout = 0.0;
  for (it = 0; it < 3; ++it) 
    {
      tout = 0.1*(it+1);

      d03psc(npde, &ts, tout, pdef1, nmflx1, bndry1, uvin1, u, 
             npts, x, ncode, d03pek, nxi, xi, neqn, rtol, atol, 
             itol, Nag_OneNorm, Nag_LinAlgBand, algopt, remesh, 
             nxfix, xfix, nrmesh, dxmesh, trmesh, ipminf, xratio, 
             con, monit1, rsave, lrsave, isave, lisave, itask, 
             itrace, 0, &ind, &comm, &saved, &fail);

      if (fail.code != NE_NOERROR)
        {
          Vprintf("Error from d03psc.\n%s\n", fail.message);
          exit_status = 1;
          goto END;
        }

      Vprintf(" t = %6.3f\n", ts);
      Vprintf(" x        ");

      for (i = 1; i <= intpts; ++i) 
        {
          Vprintf("%9.4f", xout[i-1]);
          Vprintf(i%7 == 0 || i == 7 ?"\n":"");
        }

      /* Interpolate at output points */

      m = 0;
      d03pzc(npde, m, u, npts, x, xout, intpts, itype, uout, &fail);

      if (fail.code != NE_NOERROR)
        {
          Vprintf("Error from d03pzc.\n%s\n", fail.message);
          exit_status = 1;
          goto END;
        }

      Vprintf(" Approx u ");

      for (i = 1; i <= intpts; ++i)
        {
          Vprintf("%9.4f", UOUT(1,i,1));
          Vprintf(i%7 == 0 || i == 7 ?"\n":"");
        }       
      Vprintf("\n");
    }

  Vprintf(" Number of integration steps in time = %6ld\n", isave[0]);
  Vprintf(" Number of function evaluations = %6ld\n", isave[1]);
  Vprintf(" Number of Jacobian evaluations =%6ld\n", isave[2]);
  Vprintf(" Number of iterations = %6ld\n\n", isave[4]);
 END:

  if (algopt) NAG_FREE(algopt);
  if (atol) NAG_FREE(atol);
  if (rsave) NAG_FREE(rsave);
  if (rtol) NAG_FREE(rtol);
  if (u) NAG_FREE(u);
  if (uout) NAG_FREE(uout);
  if (x) NAG_FREE(x);
  if (xfix) NAG_FREE(xfix);
  if (xi) NAG_FREE(xi);
  if (isave) NAG_FREE(isave);

  return exit_status;

}


static void uvin1(Integer npde, Integer npts, Integer nxi, const double x[],
                  const double xi[], double u[],Integer ncode, double v[],
                  Nag_Comm *comm)
{
  Integer i;

  for (i = 1; i <= npts; ++i) 
    {
      if (x[i-1] > 0.2 && x[i-1] <= 0.4) 
        {
          U(1, i) = sin(nag_pi*(5.0*x[i-1]-1.0));
        } else {
          U(1, i) = 0.0;
        }
    }
  return;
}


static void pdef1(Integer npde, double t, double x, const double u[],
                  const double ux[], Integer ncode, const double v[],
                  const double vdot[], double p[], double c[], double d[],
                  double s[], Integer *ires, Nag_Comm *comm)
{
  P(1, 1) = 1.0;
  c[0] = 0.002;
  d[0] = ux[0];
  s[0] = 0.0;

  return;
}


static void bndry1(Integer npde, Integer npts, double t, const double x[],
                   const double u[], Integer ncode, const double v[],
                   const double vdot[], Integer ibnd, double g[],
                   Integer *ires, Nag_Comm *comm)
{
  /* Zero solution at both boundaries */

  if (ibnd == 0) 
    {
      g[0] = U(1, 1);
    } else {
      g[0] = U(1, npts);
    }
  return;
}

static void monit1(double t, Integer npts, Integer npde, const double x[], 
                   const double u[], double fmon[], Nag_Comm *comm)
{
  double h1, h2, h3;
  Integer i;

  for (i = 2; i <= npts-1; ++i) 
    {
      h1 = x[i - 1] - x[i - 2];
      h2 = x[i] - x[i - 1];
      h3 = 0.5*(x[i] - x[i - 2]);

      /* Second derivatives */

      fmon[i-1] = fabs(((U(1,i+1) - U(1,i))/h2 - (U(1,i) - U(1,i-1))/h1)/h3);
    }
  fmon[0] = fmon[1];
  fmon[npts-1] = fmon[npts-2];

  return;
}


static void nmflx1(Integer npde, double t, double x, Integer ncode, 
                   const double v[], const double uleft[],
                   const double uright[], double flux[], Integer *ires,
                   Nag_Comm *comm, Nag_D03_Save *saved)
{
  flux[0] = uleft[0];

  return;
}


int ex2(void)
{
  const Integer npde=1, npts=61, ncode=0, nxi=0, nxfix=0, itype=1, 
    neqn=npde*npts+ncode, intpts=7, lisave=25+nxfix+neqn, 
    nwkres=npde*(3*npts+3*npde+32)+7*npts+3, lenode=11*neqn+50,
    mlu=3*npde-1, lrsave=(3*mlu+1)*neqn+nwkres+lenode;
  static double xout[7] = { 0.,.3,.4,.5,.6,.7,1. };
  double con, dxmesh, tout, trmesh, ts, xratio;
  Integer exit_status, i, ind, ipminf, it, itask, itol, itrace, m, 
    nrmesh;
  Boolean remesh;
  double *algopt=0, *atol=0, *rsave=0, *rtol=0, 
    *u=0, *ue=0, *uout=0, *x=0, *xfix=0, *xi=0;
  Integer *isave=0;
  NagError fail;
  Nag_Comm comm;
  Nag_D03_Save saved;

  INIT_FAIL(fail);
  exit_status = 0;

  /* Allocate memory */

  if ( !(algopt = NAG_ALLOC(30, double)) ||
       !(atol = NAG_ALLOC(1, double)) ||
       !(rsave = NAG_ALLOC(lrsave, double)) ||
       !(rtol = NAG_ALLOC(1, double)) ||
       !(u = NAG_ALLOC(npts, double)) ||
       !(ue = NAG_ALLOC(npde*intpts, double)) ||
       !(uout = NAG_ALLOC(npde*intpts*itype, double)) ||
       !(x = NAG_ALLOC(npts, double)) ||
       !(xfix = NAG_ALLOC(1, double)) ||
       !(xi = NAG_ALLOC(1, double)) ||
       !(isave = NAG_ALLOC(lisave, Integer)) )
    {
      Vprintf("Allocation failure\n");
      exit_status = 1;
      goto END;
    }
    
  Vprintf("\n\nExample 2\n\n");

  itrace = 0;
  itol = 1;
  atol[0] = 5e-4;
  rtol[0] = 0.05;
    
  Vprintf(" npts = %4ld", npts);
  Vprintf(" atol = %10.3e", atol[0]);
  Vprintf(" rtol = %10.3e\n\n", rtol[0]);

  /* Initialise mesh */

  for (i = 0; i < npts; ++i) x[i] = i/(npts-1.0);
  xfix[0] = 0.0;

  /* Set remesh parameters */

  remesh = TRUE;
  nrmesh = 5;
  dxmesh = 0.0;
  trmesh = 0.0;
  con = 1.0/(npts-1.0);
  xratio = 1.5;
  ipminf = 0;

  xi[0] = 0.0;
  ind = 0;
  itask = 1;

  for (i = 0; i < 30; ++i) algopt[i] = 0.0;

  /* Theta integration */

  algopt[0] = 2.0;
  algopt[5] = 2.0;
  algopt[6] = 2.0;

  /* Max. time step */

  algopt[12] = 0.0025;

  ts = 0.0;
  tout = 0.0;
  for (it = 0; it < 2; ++it) 
    {
      tout = 0.2*(it+1);

      d03psc(npde, &ts, tout, pdef2, nmflx2, bndry2, uvin2, u, 
             npts, x, ncode, d03pek, nxi, xi, neqn, rtol, atol, 
             itol, Nag_OneNorm, Nag_LinAlgBand, algopt, remesh, 
             nxfix, xfix, nrmesh, dxmesh, trmesh, ipminf, xratio, 
             con, monit2, rsave, lrsave, isave, lisave, itask, 
             itrace, 0, &ind, &comm, &saved, &fail);

      if (fail.code != NE_NOERROR)
        {
          Vprintf("Error from d03psc.\n%s\n", fail.message);
          exit_status = 1;
          goto END;
        }

      Vprintf(" t = %6.3f\n", ts);
      Vprintf("        x        Approx u    Exact u\n\n");

      /* Interpolate at output points */

      m = 0;
      d03pzc(npde, m, u, npts, x, xout, intpts, itype, uout, &fail);

      if (fail.code != NE_NOERROR)
        {
          Vprintf("Error from d03pzc.\n%s\n", fail.message);
          exit_status = 1;
          goto END;
        }

      /* Check against exact solution */

      exact(tout, ue, xout, npde, intpts);
      for (i = 1; i <= intpts; ++i) 
        {
          Vprintf("   %9.4f", xout[i-1]);
          Vprintf("   %9.4f", UOUT(1,i,1));
          Vprintf("   %9.4f\n", UE(1,i));
        }
    }

  Vprintf(" Number of integration steps in time = %6ld\n", isave[0]);
  Vprintf(" Number of function evaluations = %6ld\n", isave[1]);
  Vprintf(" Number of Jacobian evaluations =%6ld\n", isave[2]);
  Vprintf(" Number of iterations = %6ld\n\n", isave[4]);
 END:

  if (algopt) NAG_FREE(algopt);
  if (atol) NAG_FREE(atol);
  if (rsave) NAG_FREE(rsave);
  if (rtol) NAG_FREE(rtol);
  if (u) NAG_FREE(u);
  if (ue) NAG_FREE(ue);
  if (uout) NAG_FREE(uout);
  if (x) NAG_FREE(x);
  if (xfix) NAG_FREE(xfix);
  if (xi) NAG_FREE(xi);
  if (isave) NAG_FREE(isave);

  return exit_status;
}


static void uvin2(Integer npde, Integer npts, Integer nxi, const double x[],
                  const double xi[], double u[],Integer ncode, double v[],
                  Nag_Comm *comm)
{
  double t;

  t = 0.0;
  exact(t, u, x, npde, npts);

  return;
}


static void pdef2(Integer npde, double t, double x, const double u[],
                  const double ux[], Integer ncode, const double v[],
                  const double vdot[], double p[], double c[], double d[],
                  double s[], Integer *ires, Nag_Comm *comm)
{
  P(1, 1) = 1.0;
  c[0] = 0.0;
  d[0] = 0.0;
  s[0] = -100.0*u[0]*(u[0]-1.0)*(u[0]-0.5);

  return;
}


static void bndry2(Integer npde, Integer npts, double t, const double x[],
                   const double u[], Integer ncode, const double v[],
                   const double vdot[], Integer ibnd, double g[],
                   Integer *ires, Nag_Comm *comm)
{
  /* Solution known to be constant at both boundaries */

  double ue[1];

  if (ibnd == 0) 
    {
      exact(t, ue, &x[0], npde, 1);
      g[0] = UE(1, 1) - U(1, 1);
    } else {
      exact(t, ue, &x[npts-1], npde, 1);
      g[0] = UE(1, 1) - U(1,npts);
    }

  return;
}



static void nmflx2(Integer npde, double t, double x, Integer ncode, 
                   const double v[], const double uleft[],
                   const double uright[], double flux[], Integer *ires,
                   Nag_Comm *comm, Nag_D03_Save *saved)
{
  flux[0] = uleft[0];

  return;
}


static void monit2(double t, Integer npts, Integer npde, const double x[], 
                   const double u[], double fmon[], Nag_Comm *comm)
{
  static double xa = 0.0;
  static Integer icount = 0;
  double h1, ux, uxmax, xl, xleft, xmax, xr, xright;
  Integer i;

  /* Locate shock */

  uxmax = 0.0;
  xmax = 0.0;
  for (i = 2; i <= npts-1; ++i) 
    {
      h1 = x[i - 1] - x[i - 2];
      ux = fabs((U(1, i) - U(1, i-1))/h1);
      if (ux > uxmax) 
        {
          uxmax = ux;
          xmax = x[i - 1];
        }
    }

  /* Assign width (on first call only) */

  if (icount == 0) 
    {
      icount = 1;
      xleft = xmax - x[0];
      xright = x[npts-1] - xmax;
      if (xleft > xright) 
        {
          xa = xright;
        } else {
          xa = xleft;
        }
    }
  xl = xmax - xa;
  xr = xmax + xa;

  /* Assign monitor function */

  for (i = 0; i < npts; ++i) 
    {
      if (x[i] > xl && x[i] < xr) 
        {
          fmon[i] = 1.0 + cos(nag_pi*(x[i] - xmax)/xa);
        } else {
          fmon[i] = 0.0;
        }
    }
  return;
}


static void exact(double t, double *u, const double *x, Integer npde, 
                  Integer npts)
{
  /* Exact solution (for comparison and b.c. purposes) */

  double del, psi, rm, rn, s;
  Integer i;

  s = 0.1;
  del = 0.01;
  rm = -1.0/del;
  rn = s/del + 1.0;

  for (i = 1; i <= npts; ++i) 
    {
      psi = x[i - 1] - t;
      if (psi < s) 
        {
          U(1, i) = 1.0;
        } else if (psi > del + s) {
          U(1, i) = 0.0;
        } else {
          U(1, i) = rm*psi + rn;
        }
    }
  return;
}