// LinElas3d_WG_HexaMesh.h
// James Liu, Graham Harper, ColoState; 2014/07--2017/11

#ifndef LinElas3d_WG_HexaMesh_H
#define LinElas3d_WG_HexaMesh_H

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

#include "GaussQuad.h"
#include "PtVec3d.h"
#include "HexaMesh.h"

////////////////////////////////////////////////////////////////////////////////
// Solving linear elasticity in 3d by WG(P0^3,P0^3;AT0^3,P0) on a hexa.mesh
// See ...
////////////////////////////////////////////////////////////////////////////////
/*
// JL20171102: TO BE FINISHED
// The Schur-complement approach: 4 functions

int LinElas3d_WG_HexaP03P03AT03P0_AsmBndryConds(Vector &GlbVecD,
                                                Vector &GlbVecN,
                                                PtVec3d (*fxnuD)(PtVec3d),
                                                PtVec3d (*fxntN)(PtVec3d),
                                                const HexaMesh &mesh,
                                                const GaussQuad &GQQ);

int LinElas3d_WG_HexaP03P03AT03P0_AsmGlbMats(BlockDiagMatrix &GlbMatEE,
                                             SparseBlockMatrix &GlbMatEC,
                                             SparseBlockMatrix &GlbMatCE,
                                             SparseBlockMatrix &GlbMatCC,
                                             const double lambda,
                                             const double mu,
                                             const HexaMesh &mesh,
                                             const GaussQuad &GQH,
                                             const GaussQuad &GQQ);

int LinElas3d_WG_HexaP03P03AT03P0_AsmSource(Vector &GlbVecE,
                                            PtVec3d (*fxnf)(PtVec3d),
                                            const HexaMesh &mesh,
                                            const GaussQuad &GQH);

int LinElas3d_WG_HexaP03P03AT03P0_ModiLinSys(SparseBlockMatrix &GlbMatEC,
                                             SparseBlockMatrix &GlbMatCE,
                                             SparseBlockMatrix &GlbMatCC,
                                             Vector &GlbVecE, Vector &GlbVecC,
                                             const HexaMesh &mesh,
                                             const Vector &GlbVecD,
                                             const Vector &GlbVecN);

// Common functions

int LinElas3d_WG_HexaP03P03AT03P0_DsplStrsDiv(FullMatrix &NumerDsplCofP03,
                                              FullMatrix &NumerStrsCofSym,
                                              Vector &NumerDivCofP0,
                                              const double &lambda, const double &mu,
                                              const HexaMesh &mesh, const Vector &sln,
                                              const GaussQuad &GQH, const GaussQuad &GQQ);

int LinElas3d_WG_HexaP03P03AT03P0_Err(double &L2ErrDspl, double &L2ErrStrs, double &L2ErrDiv,
                                      PtVec3d (*fxnu)(PtVec3d), Mat3 (*fxnsigma)(PtVec3d),
                                      double (*fxndivu)(PtVec3d), const HexaMesh &mesh,
                                      const FullMatrix &NumerDsplCofP03,
                                      const FullMatrix &NumerStrsCofSym,
                                      const Vector &NumerDivCofP0,
                                      const GaussQuad &GQH, const GaussQuad &GQQ);
*/
////////////////////////////////////////////////////////////////////////////////
// Solving linear elasticity in 3d by WG(Q0^3,Q0^3;RT[0]^3,Q0) on a hexa.mesh
// See Liu_PREP_2017
////////////////////////////////////////////////////////////////////////////////

// JL20170423: TO BE FINISHED BY James
// The single-matrix approach: 4 functions with suffix "1"

int LinElas3d_WG_HexaQ03Q03RT03Q0_AsmBndryConds1(
  Vector &GlbVecDirichlet, Vector &GlbVecNeumann, Vector &flag,
  PtVec3d (*fxnuD)(PtVec3d), PtVec3d (*fxntN)(PtVec3d), const Vector &ABC,
  const HexaMesh &mesh, const GaussQuad &GQQ);

int LinElas3d_WG_HexaQ03Q03RT03Q0_AsmGlbMat1(SparseBlockMatrix &GlbMat,
                                             const double &lambda,
                                             const double &mu,
                                             const HexaMesh &mesh,
                                             const GaussQuad &GQH,
                                             const GaussQuad &GQQ);

int LinElas3d_WG_HexaQ03Q03RT03Q0_AsmSource1(Vector &GlbVecSource,
                                             PtVec3d (*fxnf)(PtVec3d),
                                             const HexaMesh &mesh,
                                             const GaussQuad &GQH);

int LinElas3d_WG_HexaQ03Q03RT03Q0_ModiLinSys1(
  SparseBlockMatrix &GlbMat, Vector &GlbRHS,
  const Vector &ABC, const HexaMesh &mesh, const Vector &flag,
  const Vector &GlbVecSource, const Vector &GlbVecDirichlet,
  const Vector &GlbVecNeumann);

// JL20170423: TO BE FINISHED BY Graham
// The Schur-complement approach: 4 functions

int LinElas3d_WG_HexaQ03Q03RT03Q0_AsmBndryConds(
  Vector &GlbVecD, Vector &GlbVecN, Vector &flag,
  PtVec3d (*fxnuD)(PtVec3d), PtVec3d (*fxntN)(PtVec3d),
  const Vector &ABC, const HexaMesh &mesh, const GaussQuad &GQQ);

int LinElas3d_WG_HexaQ03Q03RT03Q0_AsmBndryConds(Vector &GlbVecD,
                                                Vector &GlbVecN,
                                                Vector &flag,
                                                PtVec3d (*fxnuD)(PtVec3d),
                                                PtVec3d (*fxntN)(PtVec3d),
                                                const Vector &ABC,
                                                const HexaMesh &mesh,
                                                const GaussQuad &GQQ);

int LinElas3d_WG_HexaQ03Q03RT03Q0_AsmGlbMats(BlockDiagMatrix &GlbMatEE,
                                             SparseBlockMatrix &GlbMatEC,
                                             SparseBlockMatrix &GlbMatCE,
                                             SparseBlockMatrix &GlbMatCC,
                                             const double lambda,
                                             const double mu,
                                             const HexaMesh &mesh,
                                             const GaussQuad &GQH,
                                             const GaussQuad &GQQ);

int LinElas3d_WG_HexaQ03Q03RT03Q0_AsmSource(Vector &GlbVecE,
                                            PtVec3d (*fxnf)(PtVec3d),
                                            const HexaMesh &mesh,
                                            const GaussQuad &GQH);

int LinElas3d_WG_HexaQ03Q03RT03Q0_ModiLinSys(SparseBlockMatrix &GlbMatEC,
                                             SparseBlockMatrix &GlbMatCE,
                                             SparseBlockMatrix &GlbMatCC,
                                             Vector &GlbVecE, Vector &GlbVecC,
                                             const HexaMesh &mesh,
                                             const Vector &GlbVecD,
                                             const Vector &GlbVecN);

// Common functions

int LinElas3d_WG_HexaQ03Q03RT03Q0_DsplStrsDiv(FullMatrix &NumerDsplCofQ03,
                                              FullMatrix &NumerStrsCofSym,
                                              Vector &NumerDivCofQ0,
                                              const double &lambda, const double &mu,
                                              const HexaMesh &mesh, const Vector &sln,
                                              const GaussQuad &GQH, const GaussQuad &GQQ);

int LinElas3d_WG_HexaQ03Q03RT03Q0_Err(double &L2ErrDspl, double &L2ErrStrs, double &L2ErrDiv,
                                      PtVec3d (*fxnu)(PtVec3d), Mat3 (*fxnsigma)(PtVec3d),
                                      double (*fxndivu)(PtVec3d), const HexaMesh &mesh,
                                      const FullMatrix &NumerDsplCofQ03,
                                      const FullMatrix &NumerStrsCofSym,
                                      const Vector &NumerDivCofQ0,
                                      const GaussQuad &GQH, const GaussQuad &GQQ);

////////////////////////////////////////////////////////////////////////////////
// Solving linear elasticity in 3d by WG(P1^3,Prm;P0^{3x3},P0) on a hexa.mesh
// Stabilization needed
// See WangWangWangZhang_JCAM_2016
////////////////////////////////////////////////////////////////////////////////

// JL20170423: TO BE FINISHED BY Graham
// The single-matrix approach: 4 functions with suffix "1"

int LinElas3d_WG_HexaP13PrmP033P0_AsmBndryConds1(Vector &GlbVecDirichlet,
                                                 Vector &GlbVecNeumann,
                                                 PtVec3d (*fxnuD)(PtVec3d),
                                                 PtVec3d (*fxntN)(PtVec3d),
                                                 const HexaMesh &mesh,
                                                 const GaussQuad &GQQ);

int LinElas3d_WG_HexaP13PrmP033P0_AsmGlbMat1(SparseMatrix &GlbMat,
                                             const double &lambda,
                                             const double &mu,
                                             const HexaMesh &mesh,
                                             const GaussQuad &GQH,
                                             const GaussQuad &GQQ);

int LinElas3d_WG_HexaP13PrmP033P0_AsmSource1(Vector &GlbVecSource,
                                             PtVec3d (*fxnf)(PtVec3d),
                                             const HexaMesh &mesh,
                                             const GaussQuad &GQH);

int LinElas3d_WG_HexaP13PrmP033P0_ModiLinSys1(SparseMatrix &GlbMat,
                                              Vector &GlbRHS,
                                              const HexaMesh &mesh,
                                              const Vector &GlbVecSource,
                                              const Vector &GlbVecDirichlet,
                                              const Vector &GlbVecNeumann);

// JL20170423: TO BE FINISHED BY James
// The Schur-complement approach: 4 functions

int LinElas3d_WG_HexaP13PrmP033P0_AsmBndryConds(Vector &GlbVecD,
                                                Vector &GlbVecN,
                                                PtVec3d (*fxnuD)(PtVec3d),
                                                PtVec3d (*fxntN)(PtVec3d),
                                                const HexaMesh &mesh,
                                                const GaussQuad &GQQ);

int LinElas3d_WG_HexaP13PrmP033P0_AsmGlbMats(BlockDiagMatrix &GlbMatEE,
                                             SparseBlockMatrix &GlbMatEC,
                                             SparseBlockMatrix &GlbMatCE,
                                             SparseBlockMatrix &GlbMatCC,
                                             const double lambda,
                                             const double mu,
                                             const HexaMesh &mesh,
                                             const double &rho,
                                             const GaussQuad &GQH,
                                             const GaussQuad &GQQ);

int LinElas3d_WG_HexaP13PrmP033P0_AsmSource(Vector &GlbVecE,
                                            PtVec3d (*fxnf)(PtVec3d),
                                            const HexaMesh &mesh,
                                            const GaussQuad &GQH);

int LinElas3d_WG_HexaP13PrmP033P0_ModiLinSys(SparseBlockMatrix &GlbMatEC,
                                             SparseBlockMatrix &GlbMatCE,
                                             SparseBlockMatrix &GlbMatCC,
                                             Vector &GlbVecE, Vector &GlbVecC,
                                             const HexaMesh &mesh,
                                             const Vector &GlbVecD,
                                             const Vector &GlbVecN);

// Common functions 

int LinElas3d_WG_HexaP13PrmP033P0_DsplStrsDiv(FullMatrix &NumerDsplCofP13,
                                              FullMatrix &NumerStrsCofSym, Vector &NumerDivCofP0,
                                              const double &lambda, const double &mu,
                                              const HexaMesh &mesh, const Vector &sln,
                                              const GaussQuad &GQH, const GaussQuad &GQQ);

int LinElas3d_WG_HexaP13PrmP033P0_Err(double &L2ErrDspl, double &L2ErrStrs, double &L2ErrDiv,
                                      PtVec3d (*fxnu)(PtVec3d), Mat3 (*fxnsigma)(PtVec3d),
                                      double (*fxndivu)(PtVec3d), const HexaMesh &mesh,
                                      const FullMatrix &NumerDsplCofP13,
                                      const FullMatrix &NumerStrsCofSym,
                                      const Vector &NumerDivCofP0,
                                      const GaussQuad &GQH, const GaussQuad &GQQ);

#endif  // LinElas3d_WG_HexaMesh_H
// LinElas3d_WG_HexaMesh.h
