manufactured_solution.h

#ifndef __MANUFACTUREDSOLUTIONFUNCTION_H__
#define __MANUFACTUREDSOLUTIONFUNCTION_H__

#include <deal.II/lac/vector.h>

#include <deal.II/base/function.h>

//#include <Sacado.hpp>
//
//#include "physics/physics.h"
//#include "numerical_flux/numerical_flux.h"
#include "parameters/all_parameters.h"


namespace PHiLiP {



template <int dim, typename real>
class ManufacturedSolutionFunction : public dealii::Function<dim,real>
{
// We want the Point to be templated on the type,
// however, dealii does not template that part of the Function.
// Therefore, we end up overloading the functions and need to "import"
// those non-overloaded functions to avoid the warning -Woverloaded-virtual
// See: https://stackoverflow.com/questions/18515183/c-overloaded-virtual-function-warning-by-clang
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;
    using dealii::Function<dim,real>::vector_gradient;

public:
    const unsigned int nstate; 

    explicit ManufacturedSolutionFunction (const unsigned int nstate = 1);


    virtual real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const = 0;


    virtual dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const = 0;

    dealii::Tensor<1,dim,real> gradient_fd (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const;


    virtual dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const = 0;

    dealii::SymmetricTensor<2,dim,real> hessian_fd (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const;

    std::vector<real> stdvector_values (const dealii::Point<dim,real> &point) const;

  
    void vector_gradient (const dealii::Point<dim,real> &p,
                          std::vector<dealii::Tensor<1,dim, real> > &gradients) const;

    // Virtual functions inherited from dealii::Function
    //
    // virtual real value (const Point<dim,real> &p,
    //                               const unsigned int  component = 0) const;
  
    // virtual void vector_value (const Point<dim,real> &p,
    //                           Vector<real> &values) const;
  
    // virtual void value_list (const std::vector<Point<dim,real> > &points,
    //                         std::vector<real> &values,
    //                         const unsigned int              component = 0) const;
  
    // virtual void vector_value_list (const std::vector<Point<dim,real> > &points,
    //                                std::vector<Vector<real> > &values) const;
  
    // virtual void vector_values (const std::vector<Point<dim,real> > &points,
    //                            std::vector<std::vector<real> > &values) const;
  
    // virtual Tensor<1,dim, real> gradient (const Point<dim,real> &p,
    //                                                 const unsigned int  component = 0) const;
  
    // virtual void gradient_list (const std::vector<Point<dim,real> > &points,
    //                            std::vector<Tensor<1,dim, real> > &gradients,
    //                            const unsigned int              component = 0) const;
  
    // virtual void vector_gradients (const std::vector<Point<dim,real> > &points,
    //                               std::vector<std::vector<Tensor<1,dim, real> > > &gradients) const;
  
    // virtual void vector_gradient_list (const std::vector<Point<dim,real> > &points,
    //                                   std::vector<std::vector<Tensor<1,dim, real> > > &gradients) const;

protected:

    std::vector<double> base_values;
    std::vector<double> amplitudes;
    std::vector<dealii::Tensor<1,dim,real>> frequencies;
};

template <int dim, typename real>
class ManufacturedSolutionZero
        : public ManufacturedSolutionFunction<dim, real>
{
// We want the Point to be templated on the type,
// however, dealii does not template that part of the Function.
// Therefore, we end up overloading the functions and need to "import"
// those non-overloaded functions to avoid the warning -Woverloaded-virtual
// See: https://stackoverflow.com/questions/18515183/c-overloaded-virtual-function-warning-by-clang
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;

public:
    explicit ManufacturedSolutionZero(const unsigned int nstate = 1)
            :   ManufacturedSolutionFunction<dim,real>(nstate){}
    real value (const dealii::Point<dim,real> &/*point*/, const unsigned int /*istate*/ = 0) const override;
    dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &/*point*/, const unsigned int /*istate*/ = 0) const override;
    dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &/*point*/, const unsigned int /*istate*/ = 0) const override;
};

template <int dim, typename real>
class ManufacturedSolutionSine 
    : public ManufacturedSolutionFunction<dim, real>
{
// We want the Point to be templated on the type,
// however, dealii does not template that part of the Function.
// Therefore, we end up overloading the functions and need to "import"
// those non-overloaded functions to avoid the warning -Woverloaded-virtual
// See: https://stackoverflow.com/questions/18515183/c-overloaded-virtual-function-warning-by-clang
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;

public:
    explicit ManufacturedSolutionSine(const unsigned int nstate = 1)
        :   ManufacturedSolutionFunction<dim,real>(nstate){}
    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
};

template <int dim, typename real>
class ManufacturedSolutionCosine 
    : public ManufacturedSolutionFunction<dim, real>
{
// We want the Point to be templated on the type,
// however, dealii does not template that part of the Function.
// Therefore, we end up overloading the functions and need to "import"
// those non-overloaded functions to avoid the warning -Woverloaded-virtual
// See: https://stackoverflow.com/questions/18515183/c-overloaded-virtual-function-warning-by-clang
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;
public:
    explicit ManufacturedSolutionCosine(const unsigned int nstate = 1)
        :   ManufacturedSolutionFunction<dim,real>(nstate){}
    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
};

template <int dim, typename real>
class ManufacturedSolutionAdd 
    : public ManufacturedSolutionFunction<dim, real>
{
// We want the Point to be templated on the type,
// however, dealii does not template that part of the Function.
// Therefore, we end up overloading the functions and need to "import"
// those non-overloaded functions to avoid the warning -Woverloaded-virtual
// See: https://stackoverflow.com/questions/18515183/c-overloaded-virtual-function-warning-by-clang
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;
public:
    explicit ManufacturedSolutionAdd(const unsigned int nstate = 1)
        :   ManufacturedSolutionFunction<dim,real>(nstate){}
    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
};

template <int dim, typename real>
class ManufacturedSolutionExp
    : public ManufacturedSolutionFunction<dim, real>
{
// We want the Point to be templated on the type,
// however, dealii does not template that part of the Function.
// Therefore, we end up overloading the functions and need to "import"
// those non-overloaded functions to avoid the warning -Woverloaded-virtual
// See: https://stackoverflow.com/questions/18515183/c-overloaded-virtual-function-warning-by-clang
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;
public:
    explicit ManufacturedSolutionExp(const unsigned int nstate = 1)
        :   ManufacturedSolutionFunction<dim,real>(nstate){}
    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
};

template <int dim, typename real>
class ManufacturedSolutionPoly 
    : public ManufacturedSolutionFunction<dim, real>
{
// We want the Point to be templated on the type,
// however, dealii does not template that part of the Function.
// Therefore, we end up overloading the functions and need to "import"
// those non-overloaded functions to avoid the warning -Woverloaded-virtual
// See: https://stackoverflow.com/questions/18515183/c-overloaded-virtual-function-warning-by-clang
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;
public:
    explicit ManufacturedSolutionPoly(const unsigned int nstate = 1)
        :   ManufacturedSolutionFunction<dim,real>(nstate){}
    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
};

template <int dim, typename real>
class ManufacturedSolutionEvenPoly 
    : public ManufacturedSolutionFunction<dim, real>
{
// We want the Point to be templated on the type,
// however, dealii does not template that part of the Function.
// Therefore, we end up overloading the functions and need to "import"
// those non-overloaded functions to avoid the warning -Woverloaded-virtual
// See: https://stackoverflow.com/questions/18515183/c-overloaded-virtual-function-warning-by-clang
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;
public:
    explicit ManufacturedSolutionEvenPoly(const unsigned int nstate = 1)
        :   ManufacturedSolutionFunction<dim,real>(nstate){}
    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
};

template <int dim, typename real>
class ManufacturedSolutionAtan
    : public ManufacturedSolutionFunction<dim, real>
{
// We want the Point to be templated on the type,
// however, dealii does not template that part of the Function.
// Therefore, we end up overloading the functions and need to "import"
// those non-overloaded functions to avoid the warning -Woverloaded-virtual
// See: https://stackoverflow.com/questions/18515183/c-overloaded-virtual-function-warning-by-clang
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;
public:
    explicit ManufacturedSolutionAtan(const unsigned int nstate = 1)
        :   ManufacturedSolutionFunction<dim,real>(nstate)
    {
        n_shocks.resize(dim);
        S_j.resize(dim);
        x_j.resize(dim);

        for(unsigned int i = 0; i<dim; ++i){
            n_shocks[i] = 2;

            S_j[i].resize(n_shocks[i]);
            x_j[i].resize(n_shocks[i]);

            // S_j[i][0] =  10;
            // S_j[i][1] = -10;

            S_j[i][0] =  50;
            S_j[i][1] = -50;

            x_j[i][0] = -1/sqrt(2);
            x_j[i][1] =  1/sqrt(2);

            // x_j[i][0] = 1-1/sqrt(2);
            // x_j[i][1] = 1/sqrt(2);
        }
    }
    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;

private:
    std::vector<unsigned int> n_shocks; 
    std::vector<std::vector<real>> S_j; 
    std::vector<std::vector<real>> x_j; 
};

template <int dim, typename real>
class ManufacturedSolutionBoundaryLayer
    : public ManufacturedSolutionFunction<dim, real>
{
// We want the Point to be templated on the type,
// however, dealii does not template that part of the Function.
// Therefore, we end up overloading the functions and need to "import"
// those non-overloaded functions to avoid the warning -Woverloaded-virtual
// See: https://stackoverflow.com/questions/18515183/c-overloaded-virtual-function-warning-by-clang
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;
public:
    explicit ManufacturedSolutionBoundaryLayer(const unsigned int nstate = 1)
        : ManufacturedSolutionFunction<dim,real>(nstate)
        , epsilon(nstate)
    {
        for(int istate = 0; istate < (int)nstate; ++istate){
            for (int d=0; d<dim; d++){
                // epsilon[istate][d] = 0.1;   // smooth
                epsilon[istate][d] = 0.005; // strong
            }
        }
    }
    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;

private:
    std::vector<dealii::Tensor<1,dim,real>> epsilon;
};

template <int dim, typename real>
class ManufacturedSolutionSShock 
    : public ManufacturedSolutionFunction<dim, real>
{
// We want the Point to be templated on the type,
// however, dealii does not template that part of the Function.
// Therefore, we end up overloading the functions and need to "import"
// those non-overloaded functions to avoid the warning -Woverloaded-virtual
// See: https://stackoverflow.com/questions/18515183/c-overloaded-virtual-function-warning-by-clang
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;
public:
    explicit ManufacturedSolutionSShock(const unsigned int nstate = 1)
        :   ManufacturedSolutionFunction<dim,real>(nstate)
    {
        // setting constant for function
        // f(x,y) = a * tanh(b*sin(c*y + d) + e*x + f)

        // Ekelschot
        // Note: form given does not have brackets around b*(...)
        // a =  0.75;
        // b =  2.0;
        // c =  5.0;
        // d =  0.0;
        // e = -6.0;
        // f =  0.0;

        double scale_atan = 2.0;

        // shifted from [-1,1]^2 -> [0,1]
        a =   0.75;
        b =   2.0*scale_atan;
        c =  10.0;
        d =  -5.0;
        e = -12.0*scale_atan;
        f =   6.0*scale_atan;
    }
    real value(const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::Tensor<1,dim,real> gradient(const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::SymmetricTensor<2,dim,real> hessian(const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;

private:
    real a, b, c, d, e, f; 
};

template <int dim, typename real>
class ManufacturedSolutionQuadratic
    : public ManufacturedSolutionFunction<dim, real>
{
// We want the Point to be templated on the type,
// however, dealii does not template that part of the Function.
// Therefore, we end up overloading the functions and need to "import"
// those non-overloaded functions to avoid the warning -Woverloaded-virtual
// See: https://stackoverflow.com/questions/18515183/c-overloaded-virtual-function-warning-by-clang
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;
public:
    explicit ManufacturedSolutionQuadratic(const unsigned int nstate = 1)
        :   ManufacturedSolutionFunction<dim,real>(nstate)
    {
        // assigning the scaling coeffs for hessian diagonals
        for(unsigned int d = 0; d < dim; ++d){
            // diag(1, 4, 9, ...)
            alpha_diag[d] = (d+1)*(d+1);
        }
    }
    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;

private:
    std::array<real, dim> alpha_diag; 
};


template <int dim, typename real>
class ManufacturedSolutionExample
    : public ManufacturedSolutionFunction<dim, real>
{
// We want the Point to be templated on the type,
// however, dealii does not template that part of the Function.
// Therefore, we end up overloading the functions and need to "import"
// those non-overloaded functions to avoid the warning -Woverloaded-virtual
// See: https://stackoverflow.com/questions/18515183/c-overloaded-virtual-function-warning-by-clang
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;
public:
    explicit ManufacturedSolutionExample(const unsigned int nstate = 1)
        :   ManufacturedSolutionFunction<dim,real>(nstate){}
    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
};

template <int dim, typename real>
class ManufacturedSolutionNavahBase
    : public ManufacturedSolutionFunction<dim, real>
{
// We want the Point to be templated on the type,
// however, dealii does not template that part of the Function.
// Therefore, we end up overloading the functions and need to "import"
// those non-overloaded functions to avoid the warning -Woverloaded-virtual
// See: https://stackoverflow.com/questions/18515183/c-overloaded-virtual-function-warning-by-clang
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;
public:
    explicit ManufacturedSolutionNavahBase(const unsigned int nstate = 4)
        :   ManufacturedSolutionFunction<dim,real>(nstate)
    {
        // static_assert(dim==2, "ManufacturedSolutionNavahBase() should be created with dim=2");
        // static_assert(nstate==dim+2, "ManufacturedSolutionNavahBase() should be created with nstate=dim+2");
        
        const double pi = atan(1)*4.0;
        real L = 1.0;  
        c = pi/L; 
    }
    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
    dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
protected:
    std::array<dealii::Tensor<1,7,double>,5> ncm; 
    real c; 
    real primitive_value(const dealii::Point<dim,real> &point, const unsigned int istate = 0) const;
    dealii::Tensor<1,dim,real> primitive_gradient (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const;
    dealii::SymmetricTensor<2,dim,real> primitive_hessian (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const;
};

template <int dim, typename real>
class ManufacturedSolutionNavah_MS1
    : public ManufacturedSolutionNavahBase<dim, real>
{
public:
    explicit ManufacturedSolutionNavah_MS1(const unsigned int nstate = 4)
        :   ManufacturedSolutionNavahBase<dim,real>(nstate)
    {
        std::array<dealii::Tensor<1,7,double>,5> ncm; 
        /* MS-1 */
        ncm[0][0]= 1.0; ncm[0][1]=0.3; ncm[0][2]=-0.2; ncm[0][3]=0.3; ncm[0][4]=1.0; ncm[0][5]=1.0; ncm[0][6]=1.0;
        ncm[1][0]= 1.0; ncm[1][1]=0.3; ncm[1][2]= 0.3; ncm[1][3]=0.3; ncm[1][4]=3.0; ncm[1][5]=1.0; ncm[1][6]=1.0;
        ncm[2][0]= 1.0; ncm[2][1]=0.3; ncm[2][2]= 0.3; ncm[2][3]=0.3; ncm[2][4]=1.0; ncm[2][5]=1.0; ncm[2][6]=1.0;
        ncm[3][0]=18.0; ncm[3][1]=5.0; ncm[3][2]= 5.0; ncm[3][3]=0.5; ncm[3][4]=2.0; ncm[3][5]=1.0; ncm[3][6]=1.0;
        for(int j=0; j<7; j++) {
            ncm[4][j] = 0.0;
        }
        this->ncm=ncm; // done this way to minimize the use of keyword "this->"
    }
};

template <int dim, typename real>
class ManufacturedSolutionNavah_MS2
    : public ManufacturedSolutionNavahBase<dim, real>
{
public:
    explicit ManufacturedSolutionNavah_MS2(const unsigned int nstate = 4)
        :   ManufacturedSolutionNavahBase<dim,real>(nstate)
    {
        std::array<dealii::Tensor<1,7,double>,5> ncm; 
        /* MS-2 */
        ncm[0][0]=2.7; ncm[0][1]=0.9; ncm[0][2]=-0.9; ncm[0][3]=1.0; ncm[0][4]=1.5; ncm[0][5]=1.5; ncm[0][6]=1.5;
        ncm[1][0]=2.0; ncm[1][1]=0.7; ncm[1][2]= 0.7; ncm[1][3]=0.4; ncm[1][4]=1.0; ncm[1][5]=1.0; ncm[1][6]=1.0;
        ncm[2][0]=2.0; ncm[2][1]=0.4; ncm[2][2]= 0.4; ncm[2][3]=0.4; ncm[2][4]=1.0; ncm[2][5]=1.0; ncm[2][6]=1.0;
        ncm[3][0]=2.0; ncm[3][1]=1.0; ncm[3][2]= 1.0; ncm[3][3]=0.5; ncm[3][4]=1.0; ncm[3][5]=1.0; ncm[3][6]=1.5;
        for(int j=0; j<7; j++) {
            ncm[4][j] = 0.0;
        }
        this->ncm=ncm; // done this way to minimize the use of keyword "this->"
    }
};

template <int dim, typename real>
class ManufacturedSolutionNavah_MS3
    : public ManufacturedSolutionNavahBase<dim, real>
{
public:
    explicit ManufacturedSolutionNavah_MS3(const unsigned int nstate = 4)
        :   ManufacturedSolutionNavahBase<dim,real>(nstate)
    {
        std::array<dealii::Tensor<1,7,double>,5> ncm; 
        /* MS-3 */
        ncm[0][0]= 1.0; ncm[0][1]=0.1; ncm[0][2]=-0.2; ncm[0][3]=0.1; ncm[0][4]=1.0; ncm[0][5]=1.0; ncm[0][6]=1.0;
        ncm[1][0]= 2.0; ncm[1][1]=0.3; ncm[1][2]= 0.3; ncm[1][3]=0.3; ncm[1][4]=3.0; ncm[1][5]=1.0; ncm[1][6]=1.0;
        ncm[2][0]= 2.0; ncm[2][1]=0.3; ncm[2][2]= 0.3; ncm[2][3]=0.3; ncm[2][4]=1.0; ncm[2][5]=1.0; ncm[2][6]=1.0;
        ncm[3][0]=10.0; ncm[3][1]=1.0; ncm[3][2]= 1.0; ncm[3][3]=0.5; ncm[3][4]=2.0; ncm[3][5]=1.0; ncm[3][6]=1.0;
        for(int j=0; j<7; j++) {
            ncm[4][j] = 0.0;
        }
        this->ncm=ncm; // done this way to minimize the use of keyword "this->"
    }
};

template <int dim, typename real>
class ManufacturedSolutionNavah_MS4
    : public ManufacturedSolutionNavahBase<dim, real>
{
public:
    explicit ManufacturedSolutionNavah_MS4(const unsigned int nstate = 4)
        :   ManufacturedSolutionNavahBase<dim,real>(nstate)
    {
        std::array<dealii::Tensor<1,7,double>,5> ncm; 
        /* MS-4 */
        ncm[0][0]= 1.0; ncm[0][1]=  0.1; ncm[0][2]= -0.2; ncm[0][3]= 0.1; ncm[0][4]=1.0; ncm[0][5]=1.0; ncm[0][6]=1.0;
        ncm[1][0]= 2.0; ncm[1][1]=  0.3; ncm[1][2]=  0.3; ncm[1][3]= 0.3; ncm[1][4]=3.0; ncm[1][5]=1.0; ncm[1][6]=1.0;
        ncm[2][0]= 2.0; ncm[2][1]=  0.3; ncm[2][2]=  0.3; ncm[2][3]= 0.3; ncm[2][4]=1.0; ncm[2][5]=1.0; ncm[2][6]=1.0;
        ncm[3][0]=10.0; ncm[3][1]=  1.0; ncm[3][2]=  1.0; ncm[3][3]= 0.5; ncm[3][4]=2.0; ncm[3][5]=1.0; ncm[3][6]=1.0;
        ncm[4][0]= 0.6; ncm[4][1]=-0.03; ncm[4][2]=-0.02; ncm[4][3]=0.02; ncm[4][4]=2.0; ncm[4][5]=1.0; ncm[4][6]=3.0;
        this->ncm=ncm; // done this way to minimize the use of keyword "this->"
    }
};

template <int dim, typename real>
class ManufacturedSolutionNavah_MS5
    : public ManufacturedSolutionNavahBase<dim, real>
{
public:
    explicit ManufacturedSolutionNavah_MS5(const unsigned int nstate = 4)
        :   ManufacturedSolutionNavahBase<dim,real>(nstate)
    {
        std::array<dealii::Tensor<1,7,double>,5> ncm; 
        /* MS-5 */
        ncm[0][0]= 1.0; ncm[0][1]= 0.1; ncm[0][2]=-0.2; ncm[0][3]=0.1; ncm[0][4]=1.0; ncm[0][5]=1.0; ncm[0][6]=1.0;
        ncm[1][0]= 2.0; ncm[1][1]= 0.3; ncm[1][2]= 0.3; ncm[1][3]=0.3; ncm[1][4]=3.0; ncm[1][5]=1.0; ncm[1][6]=1.0;
        ncm[2][0]= 2.0; ncm[2][1]= 0.3; ncm[2][2]= 0.3; ncm[2][3]=0.3; ncm[2][4]=1.0; ncm[2][5]=1.0; ncm[2][6]=1.0;
        ncm[3][0]=10.0; ncm[3][1]= 1.0; ncm[3][2]= 1.0; ncm[3][3]=0.5; ncm[3][4]=2.0; ncm[3][5]=1.0; ncm[3][6]=1.0;
        ncm[4][0]=-6.0; ncm[4][1]=-0.3; ncm[4][2]=-0.2; ncm[4][3]=0.2; ncm[4][4]=2.0; ncm[4][5]=1.0; ncm[4][6]=3.0;
        this->ncm=ncm; // done this way to minimize the use of keyword "this->"
    }
};



template <int dim, typename real>
class ManufacturedSolutionFactory
{
    using ManufacturedSolutionEnum = Parameters::ManufacturedSolutionParam::ManufacturedSolutionType;
public:
    static std::shared_ptr< ManufacturedSolutionFunction<dim,real> > 
    create_ManufacturedSolution(
        Parameters::AllParameters const *const param, 
        int                                    nstate);

    static std::shared_ptr< ManufacturedSolutionFunction<dim,real> >
    create_ManufacturedSolution(
        ManufacturedSolutionEnum solution_type,
        int                      nstate);

};

} // namespace PHiLiP

#endif //__MANUFACTUREDSOLUTIONFUNCTION_H__