#define m 4
#define sstol 1.e-3

/*******************************************************************
 *  differential equation                                          *
 *  Q(x) du/dt + A(x,t) d2u/dx2 + B(x,t) du/dx + C(x,t) u = H(x,t) *
 *******************************************************************/
#define A(x,t) (cos(x))
#define B(x,t) (3*x*t*t)
#define C(x,t) (log(x+t))
#define Q(x)   (0.017*x)

/*************************************************
 *  exact solution (use for initial conditions)  *
 *************************************************/
#define    S(x,t) ((x*x*x + 1) * (1+exp(-t)))
#define  dSx(x,t) (3*x*x * (1+exp(-t)))
#define d2Sx(x,t) (6*x * (1+exp(-t)))
#define  dSt(x,t) (-(x*x*x + 1) * exp(-t))

/**********************
 *  forcing function  *
 **********************/
#define H(x,t) (Q(x)*dSt(x,t) + A(x,t)*d2Sx(x,t) + B(x,t)*dSx(x,t) + C(x,t)*S(x,t))

/**************************
 *  boundary conditions   *
 *  Type = 1 : Dirichlet  *
 *  Type = 2 : Neumann    *
 **************************/
#define BCTypeL    (2)
#define BCValuL(t) (dSx(1,t))
#define BCTypeR    (1)
#define BCValuR(t) (S(3,t))

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define hermite_fL(n)   (0.5  * (1+2*n)*(1+2*n)*(1-n))
#define hermite_fR(n)   (0.5  * (1-2*n)*(1-2*n)*(1+n))
#define hermite_gL(n,h) (h/8. * (2*n+1)*(2*n+1)*(2*n-1))
#define hermite_gR(n,h) (h/8. * (2*n-1)*(2*n-1)*(2*n+1))

#define hermite_dfL(n,h) (3./(2.*h) * (1+2*n) * (1-2*n))
#define hermite_dfR(n,h) (-hermite_dfL(n,h))
#define hermite_dgL(n)   (0.25 * (2*n+1) * (6*n-1))
#define hermite_dgR(n)   (0.25 * (2*n-1) * (6*n+1))

#define hermite_d2fL(n,h) (-12.*n/(h*h))
#define hermite_d2fR(n,h) (-hermite_d2fL(n,h))
#define hermite_d2gL(n,h) ((6*n+1) / h)
#define hermite_d2gR(n,h) ((6*n-1) / h)

#define L_fL(x,n,h,t) (A(x,t)*hermite_d2fL(n,h) + B(x,t)*hermite_dfL(n,h) + (C(x,t))*hermite_fL(n))
#define L_fR(x,n,h,t) (A(x,t)*hermite_d2fR(n,h) + B(x,t)*hermite_dfR(n,h) + (C(x,t))*hermite_fR(n))
#define L_gL(x,n,h,t) (A(x,t)*hermite_d2gL(n,h) + B(x,t)*hermite_dgL(n)   + (C(x,t))*hermite_gL(n,h))
#define L_gR(x,n,h,t) (A(x,t)*hermite_d2gR(n,h) + B(x,t)*hermite_dgR(n)   + (C(x,t))*hermite_gR(n,h))

#define J_fL(x,n)   (Q(x) * hermite_fL(n))
#define J_fR(x,n)   (Q(x) * hermite_fR(n))
#define J_gL(x,n,h) (Q(x) * hermite_gL(n,h))
#define J_gR(x,n,h) (Q(x) * hermite_gR(n,h))

#include "struct1d.h"
#include "fun1d.h"

main()
{
struct triblock22 G, T, M;
struct kvec2 b, bcombo, u, t1, cvec, cbar;

int     i, j, k, m2m2, m2m1, mm1, tsteps, ti;
double  x[m+1], h[m], c[2*m], cc[2*m], xi[m], midpt, 
        diff[m+1], oldU[m+1], U[m+1], dU[m+1], BCTL[2], BCTR[2], BCML[2], BCMR[2], 
        t, init_t, dt, theta, tau, taubar;

/******************
 *  define nodes  *
 ******************/
x[0] = 1.0;
x[1] = 1.3;
x[2] = 1.7;
x[3] = 2.3;
x[4] = 3.0;

/**********************
 *  define xi values  *
 **********************/
for (i = 0; i < m; i++)
    xi[i] = 1. / sqrt(12);

/******************
 *  time control  *
 ******************/
init_t = 0;
dt = 0.1;
tsteps = 200;
theta = 0.6;
tau = theta * dt;
taubar = (1-theta) * dt;

/**********************************************************
 *  define finite element lengths and collocation points  *
 **********************************************************/
j = -1;
for (i = 0; i < m; i++)
    {
    h[i] = x[i+1] - x[i];
    midpt = 0.5 * (x[i] + x[i+1]);

    c[++j] = -xi[i];
    cc[j]  = midpt - xi[i] * h[i];

    c[++j] =  xi[i];
    cc[j]  = midpt + xi[i] * h[i];
    }

/************************
 *  initial conditions  *
 ************************/
t = init_t;

/*  not BC dependent  */
for (i = 1; i <= m-1; i++)
    {
    u.el[i-1].el[1] =   S(x[i], t);
    u.el[i]  .el[0] = dSx(x[i], t);
    }

/*  BC dependent  */
if (BCTypeL == 1)  /* Dirichlet on L */
    u.el[0].el[0] = dSx(x[0], t);
else  /* Neumann on L */
    u.el[0].el[0] = S(x[0], t);

if (BCTypeR == 1)  /* Dirichlet on R */
    u.el[m-1].el[1] = dSx(x[m], t);
else  /* Neumann on R */
    u.el[m-1].el[1] = S(x[m], t);

/**************************************
 *  time-independent and initial BCs  *
 **************************************/
if (BCTypeL == 1)  /* Dirichlet on L */
    {
    BCTL[0] = J_fR(cc[0], c[0]);
    BCTL[1] = J_fR(cc[1], c[1]);
    BCML[0] = L_fR(cc[0], c[0], h[0], t);
    BCML[1] = L_fR(cc[1], c[1], h[0], t);
    }
else  /* Neumann on L */
    {
    BCTL[0] = J_gR(cc[0], c[0], h[0]);
    BCTL[1] = J_gR(cc[1], c[1], h[0]);
    BCML[0] = L_gR(cc[0], c[0], h[0], t);
    BCML[1] = L_gR(cc[1], c[1], h[0], t);
    }

m2m2 = 2*m-2;
m2m1 = m2m2 + 1;
mm1  = m-1;
if (BCTypeR == 1)  /* Dirichlet on R */
    {
    BCTR[0] = J_fL(cc[m2m2], c[m2m2]);
    BCTR[1] = J_fL(cc[m2m1], c[m2m1]);
    BCMR[0] = L_fL(cc[m2m2], c[m2m2], h[mm1], t);
    BCMR[1] = L_fL(cc[m2m1], c[m2m1], h[mm1], t);
    }
else  /* Neumann on R */
    {
    BCTR[0] = J_gL(cc[m2m2], c[m2m2], h[mm1]);
    BCTR[1] = J_gL(cc[m2m1], c[m2m1], h[mm1]);
    BCMR[0] = L_gL(cc[m2m2], c[m2m2], h[mm1], t);
    BCMR[1] = L_gL(cc[m2m1], c[m2m1], h[mm1], t);
    }

/******************************
 *  matrices T and initial M  *
 ******************************/
make_T(&T, cc, c, h);
make_M(&M, cc, c, h, t);

/**********************************************
 *  initial RHS vector at collocation points  *
 **********************************************/
k = -1;
for (i = 0; i < m; i++)
    for (j = 0; j < 2; j++)
        {
	k++;
        b.el[i].el[j] = H(cc[k], t);
        }

/******************
 *  time advance  *
 ******************/
for (ti = 1; ti <= tsteps; ti++)
    {    
    /****************************************
     *  matrix * solution at old time step  *
     ****************************************/
    saxpy_triblock22(&G, -taubar, &M, &T);
    mult_tri2xu (&t1, &G, &u);

    /***************************************
     *  forcing function at old time step  *
     ***************************************/
    for (i = 0; i < m; i++)
        for (j = 0; j < 2; j++)
	    cbar.el[i].el[j] = taubar * b.el[i].el[j];

    cbar.el[0]  .el[0] += (-taubar * BCML[0] + BCTL[0]) * BCValuL(t);
    cbar.el[0]  .el[1] += (-taubar * BCML[1] + BCTL[1]) * BCValuL(t);

    cbar.el[mm1].el[0] += (-taubar * BCMR[0] + BCTR[0]) * BCValuR(t);
    cbar.el[mm1].el[1] += (-taubar * BCMR[1] + BCTR[1]) * BCValuR(t);

    /*******************
     *  new time step  *
     *******************/
    t += dt; 

    /*******************
     *  make matrix G  *
     *******************/
    make_M(&M, cc, c, h, t);
    saxpy_triblock22(&G, tau, &M, &T);

    /**************************
     *  BCs at new time step  *
     **************************/
    if (BCTypeL == 1)  /* Dirichlet on L */
        {
	BCML[0] = L_fR(cc[0], c[0], h[0], t);
	BCML[1] = L_fR(cc[1], c[1], h[0], t);
        }
    else  /* Neumann on L */
        {
	BCML[0] = L_gR(cc[0], c[0], h[0], t);
	BCML[1] = L_gR(cc[1], c[1], h[0], t);
        }

    if (BCTypeR == 1)  /* Dirichlet on R */
        {
	BCMR[0] = L_fL(cc[m2m2], c[m2m2], h[mm1], t);
	BCMR[1] = L_fL(cc[m2m1], c[m2m1], h[mm1], t);
        }
    else  /* Neumann on R */
        {
	BCMR[0] = L_gL(cc[m2m2], c[m2m2], h[mm1], t);
	BCMR[1] = L_gL(cc[m2m1], c[m2m1], h[mm1], t);
        }

    /***********************************************
     *  evaluate RHS vector at collocation points  *
     ***********************************************/
    k = -1;
    for (i = 0; i < m; i++)
        for (j = 0; j < 2; j++)
	    {
	    k++;
	    b.el[i].el[j] = H(cc[k], t);
	    cvec.el[i].el[j] = tau * b.el[i].el[j];
	    }

    cvec.el[0]  .el[0] -= (tau * BCML[0] + BCTL[0]) * BCValuL(t);
    cvec.el[0]  .el[1] -= (tau * BCML[1] + BCTL[1]) * BCValuL(t);

    cvec.el[mm1].el[0] -= (tau * BCMR[0] + BCTR[0]) * BCValuR(t);
    cvec.el[mm1].el[1] -= (tau * BCMR[1] + BCTR[1]) * BCValuR(t);

    /***************************
     *  form final RHS vector  *
     ***************************/
    addkv2(&bcombo, &cvec, &cbar);
    addkv2(&bcombo, &t1, &bcombo);

    /************************
     *  solve: Gu = bcombo  *
     ************************/
    tridiag2(&u, &G, &bcombo);

    /*********************
     *  compute results  *
     *********************/
    /*  not BC dependent  */
    for (i = 1; i <= m-1; i++)
        {
	oldU[i] = U[i];
	 U[i] = u.el[i-1].el[1];
	dU[i] = u.el[i]  .el[0];
        }

    /*  BC dependent  */
    if (BCTypeL == 1)  /* Dirichlet on L */
        {    
	oldU[0] = U[0];
	U[0] = BCValuL(t);
	dU[0] = u.el[0]  .el[0];
        }
    else  /* Neumann on L */
        {
	oldU[0] = U[0];
	U[0] = u.el[0]  .el[0];
	dU[0] = BCValuL(t);
        }

    if (BCTypeR == 1)  /* Dirichlet on R */
        {
	oldU[m] = U[m];
	U[m] = BCValuR(t);
	dU[m] = u.el[m-1].el[1];
        }
    else  /* Neumann on R */
        {
	oldU[m] = U[m];
	U[m] = u.el[m-1].el[1];
	dU[m] = BCValuR(t);
        }

    subvec(diff, U, oldU);

    /***********************
     *  results to screen  *
     ***********************/
    blank(1);
    printf("time step = %d\ntime = %lf\ndt = %lf\n", ti, t, dt);

    for (i = 0; i <= m; i++)
        printf("%7.4lf    %8.4lf     %8.4lf\n", x[i], U[i], dU[i]);

    /*
    for (i = 0; i <= m; i++)
        printf("           %8.4lf     %8.4lf\n", S(x[i], t), dSx(x[i], t));
    blank(1);
    */

    printf("max = %lf\n", maxvec(diff));
    if (maxvec(diff) < sstol)
        {
        printf("steady state reached\n");
	ti = tsteps;
        }
    blank(1);
    }

} /* end of main */