dev_fdaPDE/Param_Functors.h at master · simonepanzeri/dev_fdaPDE · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
//
// Created by simonepanzeri on 26/11/2021.
//

#ifndef DEV_FDAPDE_PARAM_FUNCTORS_H
#define DEV_FDAPDE_PARAM_FUNCTORS_H

#include "Pde_Expression_Templates.h"

// Forward declaration!
template <UInt ORDER, UInt mydim, UInt ndim>
class FiniteElement;

// Convenience enum for options
enum class PDEParameterOptions{Constant, SpaceVarying};

template<PDEParameterOptions OPTION>
class Diffusion{
public:

    Diffusion(const Real* const K_ptr) :
            K_ptr_(K_ptr) {}

    //Diffusion(SEXP RGlobalVector) :
            //K_ptr_(REAL(RGlobalVector)) {}

    template<UInt ORDER, UInt mydim, UInt ndim>
    Real operator() (const FiniteElement<ORDER,mydim,ndim>& fe_, UInt iq, UInt i, UInt j) const;

private:
    const Real* const K_ptr_;
};

template<>
template<UInt ORDER, UInt mydim, UInt ndim>
Real Diffusion<PDEParameterOptions::Constant>::operator() (const FiniteElement<ORDER,mydim,ndim>& fe_, UInt iq, UInt i, UInt j) const {
    using EigenMap2Diff_matr = Eigen::Map<const Eigen::Matrix<Real,ndim,ndim> >;

    return fe_.stiff_impl(iq, i, j, EigenMap2Diff_matr(K_ptr_));
}

template<>
template<UInt ORDER, UInt mydim, UInt ndim>
Real Diffusion<PDEParameterOptions::SpaceVarying>::operator() (const FiniteElement<ORDER,mydim,ndim>& fe_, UInt iq, UInt i, UInt j) const {
    using EigenMap2Diff_matr = Eigen::Map<const Eigen::Matrix<Real,ndim,ndim> >;

    const UInt index = fe_.getGlobalIndex(iq) * EigenMap2Diff_matr::SizeAtCompileTime;
    return fe_.stiff_impl(iq, i, j, EigenMap2Diff_matr(&K_ptr_[index]));
}

template<PDEParameterOptions OPTION>
class Advection{
public:

    Advection(const Real* const b_ptr) :
            b_ptr_(b_ptr) {}

    //Advection(SEXP RGlobalVector) :
            //b_ptr_(REAL(RGlobalVector)) {}

    template<UInt ORDER, UInt mydim, UInt ndim>
    Real operator() (const FiniteElement<ORDER,mydim,ndim>& fe_, UInt iq, UInt i, UInt j) const;

    EOExpr<const Advection&> dot(const EOExpr<Grad>& grad) const {
        typedef EOExpr<const Advection&> ExprT;
        return ExprT(*this);
    }

private:
    const Real* const b_ptr_;
};

template<>
template<UInt ORDER, UInt mydim, UInt ndim>
Real Advection<PDEParameterOptions::Constant>::operator() (const FiniteElement<ORDER,mydim,ndim>& fe_, UInt iq, UInt i, UInt j) const {
    using EigenMap2Adv_vec = Eigen::Map<const Eigen::Matrix<Real,ndim,1> >;

    return fe_.grad_impl(iq, i, j, EigenMap2Adv_vec(b_ptr_));
}

template<>
template<UInt ORDER, UInt mydim, UInt ndim>
Real Advection<PDEParameterOptions::SpaceVarying>::operator() (const FiniteElement<ORDER,mydim,ndim>& fe_, UInt iq, UInt i, UInt j) const {
    using EigenMap2Adv_vec = Eigen::Map<const Eigen::Matrix<Real,ndim,1> >;

    const UInt index = fe_.getGlobalIndex(iq) * EigenMap2Adv_vec::SizeAtCompileTime;
    return fe_.grad_impl(iq, i, j, EigenMap2Adv_vec(&b_ptr_[index]));
}


class Reaction{
public:

    Reaction(const Real* const  c_ptr) :
            c_ptr_(c_ptr) {}

    //Reaction(SEXP RGlobalVector) :
            //c_ptr_(REAL(RGlobalVector)) {}

    template<UInt ORDER, UInt mydim, UInt ndim>
    Real operator() (const FiniteElement<ORDER, mydim, ndim>& fe_, UInt iq, UInt i, UInt j) const {
        const UInt index = fe_.getGlobalIndex(iq);
        return c_ptr_[index]*fe_.mass_impl(iq, i, j);
    }

    EOExpr<const Reaction&> operator* (const EOExpr<Mass>&  mass) const {
        typedef EOExpr<const Reaction&> ExprT;
        return ExprT(*this);
    }

private:
    const Real* const c_ptr_;
};


class ForcingTerm{
public:

    ForcingTerm() :
            u_ptr_(nullptr) {}
    ForcingTerm(const Real* const u_ptr) :
            u_ptr_(u_ptr) {}

    //ForcingTerm(SEXP RGlobalVector) :
            //u_ptr_(REAL(RGlobalVector)) {}

    template<UInt ORDER, UInt mydim, UInt ndim>
    Real integrate (const FiniteElement<ORDER, mydim, ndim>& fe_, UInt i) const {
        const UInt index = fe_.getGlobalIndex(0);
        return fe_.forcing_integrate(i, &u_ptr_[index]);
    }

private:
    const Real* const u_ptr_;
};

#endif //DEV_FDAPDE_PARAM_FUNCTORS_H