Am Montag, 13. Mai 2022, finden Wartungsarbeiten am Gitlab-Server (Update auf neue Version statt). Der Dienst wird daher am Montag für einige Zeit nicht verfügbar sein.
On Monday, May 13th 2022, the Gitlab server will be updated. The service will therefore not be accessible for some time on Monday.

ConvectionDiffusionOperator.hpp 9.57 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#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
   * @{
   **/

17
  /// convection-diffusion operator, see \ref convectionDiffusion
18
19
  /// <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
      : 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)});

74
      using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::LocalContext::mydimension>;
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
      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
96
        using WorldVector = FieldVector<RangeFieldType,Context::dow>;
97
98
        thread_local std::vector<WorldVector> gradients;
        gradients.resize(shapeGradients.size());
99
100
101
102
103

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

        const auto A = localFctA(local);
104
        WorldVector b; b = localFctB(local);
105
106
107
        const auto c = localFctC(local);

        IF_CONSTEXPR(conserving) {
108
          WorldVector gradAi, gradBi;
109
110
          for (std::size_t i = 0; i < localFE.size(); ++i) {
            const int local_i = rowNode.localIndex(i);
111
112
            gradAi = A * gradients[i];
            gradBi = b * gradients[i];
113
114
115
116
117
118
            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 {
119
          WorldVector grad_i;
120
121
          for (std::size_t i = 0; i < localFE.size(); ++i) {
            const int local_i = rowNode.localIndex(i);
122
123
            grad_i = A * gradients[i];
            grad_i+= b * shapeValues[i];
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
            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)
    {
143
      static_assert(Node::isLeaf,
144
145
        "Operator can be applied to Leaf-Nodes only." );

146
      using RangeType = typename FiniteElementType_t<Node>::Traits::LocalBasisType::Traits::RangeType;
147
148
149
150
151
152
153

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

      auto const& localFE = node.finiteElement();

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

154
      using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::LocalContext::dimension>;
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
      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) {
171
172
          const int local_i = node.localIndex(i);
          elementVector[local_i] += f * shapeValues[i] * factor;
173
174
175
176
177
178
179
180
181
182
183
184
185
186
        }
      }

      localFctF.unbind();
    }

  private:

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

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

  private:
194
195
196
197
    GridFctA gridFctA_;
    GridFctB gridFctB_;
    GridFctC gridFctC_;
    GridFctF gridFctF_;
198
199
  };

200
  template <class PreGridFctA, class PreGridFctB, class PreGridFctC, class PreGridFctF, class c>
201
202
  struct PreConvectionDiffusionOperator
  {
203
    static constexpr bool conserving = c::value;
204
205
206
207
208
    PreGridFctA gridFctA;
    PreGridFctB gridFctB;
    PreGridFctC gridFctC;
    PreGridFctF gridFctF;
  };
209

210
211
212
  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,
213
214
                           bool_t<conserving> = {})
  {
215
    using Pre = PreConvectionDiffusionOperator<PreGridFctA, PreGridFctB, PreGridFctC, PreGridFctF, bool_t<conserving>>;
216
217
218
    return Pre{gridFctA, gridFctB, gridFctC, gridFctF};
  }

219
220
  template <class LocalContext, class... T, class GridView>
  auto makeLocalOperator(PreConvectionDiffusionOperator<T...> const& pre, GridView const& gridView)
221
  {
222
    using Pre = PreConvectionDiffusionOperator<T...>;
223
224
225
226
227
228
    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,
229
      decltype(gridFctA), decltype(gridFctB), decltype(gridFctC), decltype(gridFctF), Pre::conserving>;
230
231
232

    GridFctOp localOperator{gridFctA, gridFctB, gridFctC, gridFctF};
    return localOperator;
233
234
235
236
237
  }

  /** @} **/

} // end namespace AMDiS