/*************************************************************************
 * FILE: flux_template.c                                                 *
 * --------------------------------------------------------------------- *
 * DESCRIPTION:                                                          *
 *   Implements flux routines to compute this, that, and the other       *
 *   thing.                                                              *
 * --------------------------------------------------------------------- *
 * EXTERNAL FUNCTIONS:                                                   *
 *   flux1: Compute this.                                                *
 *   flux2: Compute that.                                                *
 *   flux3: Compute the other thing.                                     *
 *************************************************************************/

#include <math.h>
#include <stdio.h>
#include "prophetc.h"
#include "grid.h"
#include "assemblist.h"
#include "pdeterms.h"
#include "mathpack.h"

/*************************************************************************
 * EXTERNAL FUNCTION: flux1                                              *
 * --------------------------------------------------------------------- *
 * Compute this. Here are the details:                                   *
 *                                                                       *
 *    this = alpha * beta * gamma * that * grad (the_other)              *
 *                                                                       *
 * Read variables: [that,the_other]                                      *
 * Write variables: [this]                                               *
 *                                                                       *
 * Database values used:                                                 *
 *    /library/physics/<material>/alpha                                  *
 *    /library/physics/<material>/this/beta                              *
 *    /library/physics/<material>/this/gamma.that                        *
 *************************************************************************/

flux1 (arglist)
    argdescrip
{

    int d, in, rhs = ((*imtx)%10 == 1), mtx = ((*imtx)/10 == 1);
    real **sol_, ***gradsol_, ***f_, ****df_, *****dgf_;
    real conc, grad, alpha, beta, gamma;
    property *tmpp; /* remove this? */

    if (*nsol != 1 || *ndep != 2) {
        ndberr ("flux1: expecting two inputs and one output\n");
        fluxwhine (nsol, msol, ndep, mdep);
        return (-1);
    }

    switch (*path) {

        case FT_CONFIG:
            df[0] = 1;        /* conc(a) */
            df[1] = 2;        /* grad(b) */
            return (0);
      
        case FT_DT:

            /* read values from the database (into static variables)
                that won't change during a time step */
            return (0);

        case -FT_DT:

            /* free any memory allocated in FT_DT */
            return (0);

        case FT_RUN:

            if (rhs || mtx) {
                sol_     = array2( *ndep, *nn, sol);
                gradsol_ = array3( *ndep, *dim, *nn, gradsol);
                if (rhs) { f_ = array3( *nsol, *dim, *nn, f); }
                if (mtx) {
                    df_  = array4( *ndep, *nsol, *dim, *nn, df);
                    dgf_ = array5( *ndep, *nsol, *dim, *dim, *nn, dgf);
                }

                /* Compute flux, Jacobian */
                alpha = rcoeff ("alpha", *ireg);
                if (alpha == 0) alpha = 1; /* Probably should complain */
                beta = rscoeff ("beta", mdep[0], *ireg);
                if (beta == 0) beta = 1; /* Probably should complain */
                alpha = rscoeff2 ("gamma", mdep[0], mdep[1], *ireg);
                if (gamma == 0) gamma = 1; /* Probably should complain */

                for (d = 0; d < *dim; d++) {
                    for (in = 0; in < *nn; in++) {
                        conc = sol_ [0][in]; grad = gradsol_ [1][d][in];
                        if (rhs) {
                            f_ [0][d][in] = alpha*beta*gamma * conc * grad;
                        }
                        if (mtx) {
                            df_ [0][0][d][in] = alpha*beta*gamma * grad;
                            dgf_ [1][0][d][d][in] = alpha*beta*gamma * conc;
                        }
                    }
                }
            }

            (void)array2free( *ndep, *nn, sol_);
            (void)array3free( *ndep, *dim, *nn, gradsol_);
            if (rhs) { (void)array3free( *nsol, *dim, *nn, f_); }
            if (mtx) {
                (void)array4free( *ndep, *nsol, *dim, *nn, df_);
                (void)array5free( *ndep, *nsol, *dim, *dim, *nn, dgf_);
            }
            break;

        default:
            /* Ignore all other calls */
            break;
    }

    return(0);
}

/* Similarly for flux2 and flux3 */