// Darcy3d_WG_TetraMesh.h
// James Liu, Graham Harper, ColoState; 2014/07--2018/05

// Solving Darcy in 3d on a tetrahedral mesh:
// -- Using WG(P0,P0;RT0) elements:     Good
// -- Using WG(P1,P0;P0^3)+stabilizer:  Theoretically good

// #ifndef Darcy3d_WG_TetraMesh_H
// #define Darcy3d_WG_TetraMesh_H

#include <vector>

#include "matrix.h"
#include "vector.h"

#include "GaussQuad.h"
#include "mat3.h"
#include "PtVec3d.h"
#include "TetraMesh.h"

////////////////////////////////////////////////////////////////////////////////
// For solving Darcy in 3d by the lowest order WG(P0,P0;RT0)Tetra elements
////////////////////////////////////////////////////////////////////////////////

// The single-matrix approach: 4 functions with suffix "1"

int Darcy3d_WG_TetraP0P0RT0_AsmBndryConds1(Vector &GlbVecDirichlet,
                                           Vector &GlbVecNeumann,
                                           double (*fxnpD)(PtVec3d),
                                           double (*fxnuN)(PtVec3d),
                                           const TetraMesh &mesh,
                                           const GaussQuad &GQT);

int Darcy3d_WG_TetraP0P0RT0_AsmGlbMat1(SparseMatrix &GlbMat,
                                       const TetraMesh &mesh, Mat3 *PermK,
                                       const GaussQuad &GQTe);

int Darcy3d_WG_TetraP0P0RT0_AsmSource1(Vector &GlbVecSource,
                                       double (*fxnf)(PtVec3d),
                                       const TetraMesh &mesh,
                                       const GaussQuad &GQTe);

int Darcy3d_WG_TetraP0P0RT0_ModiLinSys1(SparseMatrix &GlbMat, Vector &GlbRHS,
                                        const Vector &GlbVecSource,
                                        const Vector &GlbVecDirichlet,
                                        const Vector &GlbVecNeumann,
                                        const TetraMesh &mesh);

// The Schur-complement approach: 4 functions

int Darcy3d_WG_TetraP0P0RT0_AsmGlbMats(BlockDiagMatrix &GlbMatEE,
                                       SparseBlockMatrix &GlbMatEC,
                                       SparseBlockMatrix &GlbMatCE,
                                       SparseBlockMatrix &GlbMatCC,
                                       const TetraMesh &mesh, Mat3 *PermK,
                                       const GaussQuad &GQTe,
                                       const GaussQuad &GQT);

int Darcy3d_WG_TetraP0P0RT0_AsmBndryConds(Vector &GlbVecD, Vector &GlbVecN,
                                          double (*fxnpD)(PtVec3d),
                                          double (*fxnuN)(PtVec3d),
                                          const TetraMesh &mesh,
                                          const GaussQuad &GQT);

int Darcy3d_WG_TetraP0P0RT0_AsmSource(Vector &GlbVecE, double (*fxnf)(PtVec3d),
                                      const TetraMesh &mesh, const GaussQuad &GQTe);

int Darcy3d_WG_TetraP0P0RT0_ModiLinSys(SparseBlockMatrix &GlbMatEF,
                                       SparseBlockMatrix &GlbMatFE,
                                       SparseBlockMatrix &GlbMatFF,
                                       Vector &GlbVecE, Vector &GlbVecF,
                                       Vector &GlbVecD, Vector &GlbVecN,
                                       const TetraMesh &mesh);

// Common functions

int Darcy3d_WG_TetraP0P0RT0_Err(double &L2ErrPres, double &L2ErrVel,
                                const TetraMesh &mesh,
                                const Vector &NumerPresEm,
                                const FullMatrix &NumerVelCofRT0,
                                double (*fxnp)(PtVec3d),
                                PtVec3d (*fxnu)(PtVec3d));

int Darcy3d_WG_TetraP0P0RT0_PresVelFlux(Vector &NumerPresEm,
                                        FullMatrix &NumerVelCofRT0,
                                        FullMatrix &NumerFlux,
                                        const TetraMesh &mesh,
                                        Mat3 *PermK, Vector &sln,
                                        const GaussQuad &GQTe,
                                        const GaussQuad &GQT);

int Darcy3d_WG_TetraP0P0RT0_ProjPres(Vector &ProjPresEm, Vector &ProjPresFc,
                                     double (*fxnp)(PtVec3d),
                                     double (*fxnpD)(PtVec3d),
                                     const TetraMesh &mesh,
                                     const GaussQuad &GQTe,
                                     const GaussQuad &GQT);

////////////////////////////////////////////////////////////////////////////////
// For solving Darcy in 3d by WG(P1,P0;P0^3)Tetra elements with a stabilizer
////////////////////////////////////////////////////////////////////////////////

// The single-matrix approach: 4 functions with suffix "1"
/*
int Darcy3d_WG_TetraP1P0P03_AdjLinSys1(SparseMatrix &GlbMat, Vector &GlbRHS,
                                       const Vector &GlbRhsSource,
                                       const Vector &GlbRhsDirichlet,
                                       const Vector &GlbRhsNeumann,
                                       const TetraMesh &mesh);

int Darcy3d_WG_TetraP1P0P03_AsmBndryConds1(Vector &GlbVecDirichlet,
                                           Vector &GlbVecNeumann,
                                           double (*fxnpD)(PtVec3d),
                                           double (*fxnuN)(PtVec3d),
                                           const TetraMesh &mesh,
                                           const GaussQuad &GQTe);

int Darcy3d_WG_TetraP1P0P03_AsmGlbMat1(SparseMatrix &GlbMat,
                                       const TetraMesh &mesh, Mat3 *PermK,
                                       double rho,
                                       const GaussQuad &GQTe,
                                       const GaussQuad &GQT);

int Darcy3d_WG_TetraP1P0P03_AsmSource1(Vector &GlbRhsSource,
                                       double (*fxnf)(PtVec3d),
                                       const TetraMesh &mesh,
                                       const GaussQuad &GQTe);
*/

// Common functions
/*
int Darcy3d_WG_TetraP1P0P03_ProjPresVel(Vector &CofProjPresEm,
                                        Vector &CofProjPresFc,
                                        double (*fxnp)(PtVec3d),
                                        PtVec3d (*fxnu)(PtVec3d),
                                        const TetraMesh &mesh,
                                        const GaussQuad &GQTe,
                                        const GaussQuad &GQT);
*/

// #endif  // Darcy3d_WG_TetraMesh_H
