#ifndef PIPE_H
#define PIPE_H

#define GRAV                9.806
#define PIPE_CHARACTERISTIC 10000000.0

#include <petsc.h>

typedef struct {
  PetscScalar q; /* flow rate */
  PetscScalar h; /* pressure */
} PipeField;

typedef struct {
  PetscScalar Q0, H0; /* boundary values in upstream */
  PetscScalar QL, HL; /* boundary values in downstream */
} PipeBoundary;

/* pipe                 */
/*----------------------*/
struct _p_Pipe {
  /* identification variables */
  PetscInt id;
  PetscInt networkid; /* which network this pipe belongs */

  /* solver objects */
  Vec        x;
  PipeField *xold;
  PetscReal  dt;
  DM         da;
  PetscInt   nnodes; /* number of nodes in da discretization */
  Mat       *jacobian;

  /* physics */
  PetscReal    length; /* pipe length */
  PetscReal    a;      /* natural flow speed */
  PetscReal    fric;   /* friction */
  PetscReal    D;      /* diameter */
  PetscReal    A;      /* area of cross section */
  PetscReal    R;
  PetscReal    rad;
  PipeBoundary boundary; /* boundary conditions for H and Q */
} PETSC_ATTRIBUTEALIGNED(PetscMax(sizeof(double), sizeof(PetscScalar)));

typedef struct _p_Pipe *Pipe;

extern PetscErrorCode PipeCreate(MPI_Comm, Pipe *);
extern PetscErrorCode PipeDestroy(Pipe *);
extern PetscErrorCode PipeSetParameters(Pipe, PetscReal, PetscReal, PetscReal, PetscReal);
extern PetscErrorCode PipeSetUp(Pipe);
extern PetscErrorCode PipeCreateJacobian(Pipe, Mat *, Mat *[]);
extern PetscErrorCode PipeDestroyJacobian(Pipe);

extern PetscErrorCode PipeComputeSteadyState(Pipe, PetscScalar, PetscScalar);
extern PetscErrorCode PipeIFunctionLocal(DMDALocalInfo *, PetscReal, PipeField *, PipeField *, PipeField *, Pipe);
extern PetscErrorCode PipeIFunctionLocal_Lax(DMDALocalInfo *, PetscReal, PipeField *, PipeField *, PetscScalar *, Pipe);
extern PetscErrorCode PipeRHSFunctionLocal(DMDALocalInfo *, PetscReal, PipeField *, PetscScalar *, Pipe);
extern PetscErrorCode PipeMonitor(TS, PetscInt, PetscReal, Vec, void *);

extern PetscErrorCode PipeCreateJacobian(Pipe, Mat *, Mat *[]);
extern PetscErrorCode PipeDestroyJacobian(Pipe);
#endif