ConvectionDiffusionOperator.hpp 9.37 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#pragma once

#include <type_traits>

#include <amdis/LocalOperator.hpp>
#include <amdis/Output.hpp>
#include <amdis/common/Utility.hpp>
#include <amdis/common/ValueCategory.hpp>

namespace AMDiS
{
  /**
   * \addtogroup operators
   * @{
   **/

  /// convection-diffusion operator
  /// <A*grad(u),grad(v)> - <b*u, grad(v)> + <c*u, v> = <f, v> (conserving) or
  /// <A*grad(u),grad(v)> + <b*grad(u), v> + <c*u, v> = <f, v> (non conserving)
20
  template <class LocalContext, class GridFctA, class GridFctB, class GridFctC, class GridFctF, bool conserving = true>
21
  class ConvectionDiffusionOperator
22
23
      : public LocalOperator<ConvectionDiffusionOperator<LocalContext, GridFctA, GridFctB, GridFctC, GridFctF, conserving>,
                             LocalContext>
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
  {
    using A_range_type = typename GridFctA::Range;
    static_assert( Category::Scalar<A_range_type> || Category::Matrix<A_range_type>,
      "Expression A(x) must be of scalar or matrix type." );
    using b_range_type = typename GridFctB::Range;
    static_assert( Category::Scalar<b_range_type> || Category::Vector<b_range_type>,
      "Expression b(x) must be of scalar or vector type." );
    using c_range_type = typename GridFctC::Range;
    static_assert( Category::Scalar<c_range_type>,
      "Expression c(x) must be of scalar type." );
    using f_range_type = typename GridFctF::Range;
    static_assert( Category::Scalar<f_range_type>,
      "Expression f(x) must be of scalar type." );

  public:
    ConvectionDiffusionOperator(GridFctA const& gridFctA, GridFctB const& gridFctB,
40
                                GridFctC const& gridFctC, GridFctF const& gridFctF)
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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
      : gridFctA_(gridFctA)
      , gridFctB_(gridFctB)
      , gridFctC_(gridFctC)
      , gridFctF_(gridFctF)
    {}

    template <class Context, class RowNode, class ColNode, class ElementMatrix>
    void getElementMatrix(Context const& context,
                          RowNode const& rowNode, ColNode const& colNode,
                          ElementMatrix& elementMatrix)
    {
      static_assert(RowNode::isLeaf && ColNode::isLeaf,
        "Operator can be applied to Leaf-Nodes only." );

      static_assert(std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>{},
        "Galerkin operator requires equal finite elements for test and trial space." );

      using LocalBasisType = typename FiniteElementType_t<RowNode>::Traits::LocalBasisType;

      using RangeType = typename LocalBasisType::Traits::RangeType;
      using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
      using JacobianType = typename LocalBasisType::Traits::JacobianType;

      auto localFctA = localFunction(gridFctA_); localFctA.bind(context.element());
      auto localFctB = localFunction(gridFctB_); localFctB.bind(context.element());
      auto localFctC = localFunction(gridFctC_); localFctC.bind(context.element());

      auto const& localFE = rowNode.finiteElement();

      int quadDeg = std::max({this->getDegree(2,coeffOrder(localFctA),rowNode,colNode),
                              this->getDegree(1,coeffOrder(localFctB),rowNode,colNode),
                              this->getDegree(0,coeffOrder(localFctC),rowNode,colNode)});

      using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::dimension>;
      auto const& quad = QuadratureRules::rule(context.type(), quadDeg);

      for (auto const& qp : quad) {
        // Position of the current quadrature point in the reference element
        decltype(auto) local = context.local(qp.position());

        // The transposed inverse Jacobian of the map from the reference element to the element
        const auto jacobian = context.geometry().jacobianInverseTransposed(local);

        // The multiplicative factor in the integral transformation formula
        const auto factor = context.integrationElement(qp.position()) * qp.weight();

        // the values of the shape functions on the reference element at the quadrature point
        thread_local std::vector<RangeType> shapeValues;
        localFE.localBasis().evaluateFunction(local, shapeValues);

        // The gradients of the shape functions on the reference element
        thread_local std::vector<JacobianType> shapeGradients;
        localFE.localBasis().evaluateJacobian(local, shapeGradients);

        // Compute the shape function gradients on the real element
        thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients(shapeGradients.size());

        for (std::size_t i = 0; i < gradients.size(); ++i)
          jacobian.mv(shapeGradients[i][0], gradients[i]);

        const auto A = localFctA(local);
        const auto b = localFctB(local);
        const auto c = localFctC(local);

        IF_CONSTEXPR(conserving) {
          for (std::size_t i = 0; i < localFE.size(); ++i) {
            const int local_i = rowNode.localIndex(i);
            const auto gradAi = A * gradients[i];
            const auto gradBi = b * gradients[i];
            for (std::size_t j = 0; j < localFE.size(); ++j) {
              const int local_j = colNode.localIndex(j);
              elementMatrix[local_i][local_j] += (dot(gradAi, gradients[j]) + (c * shapeValues[i] - gradBi) * shapeValues[j]) * factor;
            }
          }
        } else {
          for (std::size_t i = 0; i < localFE.size(); ++i) {
            const int local_i = rowNode.localIndex(i);
            const auto grad_i = A * gradients[i];
            grad_i += b * shapeValues[i];
            for (std::size_t j = 0; j < localFE.size(); ++j) {
              const int local_j = colNode.localIndex(j);
              elementMatrix[local_i][local_j] += (dot(grad_i, gradients[j]) + c * shapeValues[i] * shapeValues[j]) * factor;
            }
          }
        }
      }

      localFctA.unbind();
      localFctB.unbind();
      localFctC.unbind();
    }


    template <class Context, class Node, class ElementVector>
    void getElementVector(Context const& context,
                          Node const& node,
                          ElementVector& elementVector)
    {
      static_assert(Node::isLeaf
        "Operator can be applied to Leaf-Nodes only." );

      using RangeType = typename FiniteElementType_t<RowNode>::Traits::LocalBasisType::Traits::RangeType;

      auto localFctF = localFunction(gridFctF_); localFctF.bind(context.element());

      auto const& localFE = node.finiteElement();

      int quad_order = this->getDegree(0,coeffOrder(localFctF),node);

      using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::dimension>;
      auto const& quad = QuadratureRules::rule(context.type(), quad_order);

      for (auto const& qp : quad) {
        // Position of the current quadrature point in the reference element
        decltype(auto) local = context.local(qp.position());

        // The multiplicative factor in the integral transformation formula
        const auto factor = context.integrationElement(qp.position()) * qp.weight();

        // the values of the shape functions on the reference element at the quadrature point
        std::vector<RangeType> shapeValues;
        localFE.localBasis().evaluateFunction(local, shapeValues);

        const auto f = localFctF(local);

        for (std::size_t i = 0; i < localFE.size(); ++i) {
          const int local_i = rowNode.localIndex(i);
          elementVector[local_i][local_j] += f * shapeValues[i] * factor;
        }
      }

      localFctF.unbind();
    }

  private:

    template <class LF>
    using HasLocalFunctionOrder = decltype( order(std::declval<LF>()) );

    template <class LocalFct>
    int coeffOrder(LocalFct const& localFct)
    {
      using Concept = Dune::Std::is_detected<HasLocalFunctionOrder, LocalFunction>;
      return Dune::Hybrid::ifElse(Concept{},
        [&](auto id) { return order(id(localFct)); },
        [] (auto)    { return 0; });
    }

  private:
    GridFctA gridFctA;
    GridFctB gridFctB;
    GridFctC gridFctC;
    GridFctF gridFctF;
  };

196
197
198
199
200
201
202
203
  template <class PreGridFctA, class PreGridFctB, class PreGridFctC, class PreGridFctF, bool conserving>
  struct PreConvectionDiffusionOperator
  {
    PreGridFctA gridFctA;
    PreGridFctB gridFctB;
    PreGridFctC gridFctC;
    PreGridFctF gridFctF;
  };
204

205
206
207
  template <class PreGridFctA, class PreGridFctB, class PreGridFctC, class PreGridFctF, bool conserving = true>
  auto convectionDiffusion(PreGridFctA const& gridFctA, PreGridFctB const& gridFctB,
                           PreGridFctC const& gridFctC, PreGridFctF const& gridFctF,
208
209
                           bool_t<conserving> = {})
  {
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
    using Pre = PreConvectionDiffusionOperator<PreGridFctA, PreGridFctB, PreGridFctC, PreGridFctF, conserving>;
    return Pre{gridFctA, gridFctB, gridFctC, gridFctF};
  }

  template <class LocalContext, class... GrdFcts, bool conserving, class GridView>
  auto makeLocalOperator(PreConvectionDiffusionOperator<GridFcts..., conserving> const& pre, GridView const& gridView)
  {
    auto gridFctA = makeGridFunction(pre.gridFctA, gridView);
    auto gridFctB = makeGridFunction(pre.gridFctB, gridView);
    auto gridFctC = makeGridFunction(pre.gridFctC, gridView);
    auto gridFctF = makeGridFunction(pre.gridFctF, gridView);

    using GridFctOp = ConvectionDiffusionOperator<LocalContext,
      decltype(gridFctA), decltype(gridFctB), decltype(gridFctC), decltype(gridFctF), conserving>;

    GridFctOp localOperator{gridFctA, gridFctB, gridFctC, gridFctF};
    return localOperator;
227
228
229
230
231
  }

  /** @} **/

} // end namespace AMDiS