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.72 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
      : 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 RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;

      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();
66
      std::size_t numLocalFe = localFE.size();
67
68
69
70
71

      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)});

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

75
76
77
78
79
      LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
      auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(context, quad);
      auto const& shapeValuesCache = cache.evaluateFunctionAtQp(context, quad);

      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
80
        // Position of the current quadrature point in the reference element
81
        decltype(auto) local = context.local(quad[iq].position());
82
83
84
85
86

        // 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
87
        const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
88
89

        // the values of the shape functions on the reference element at the quadrature point
90
        auto const& shapeValues = shapeValuesCache[iq];
91
92

        // The gradients of the shape functions on the reference element
93
        auto const& shapeGradients = shapeGradientsCache[iq];
94
95

        // 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
          for (std::size_t i = 0; i < numLocalFe; ++i) {
110
            const int local_i = rowNode.localIndex(i);
111
112
            gradAi = A * gradients[i];
            gradBi = b * gradients[i];
113
            for (std::size_t j = 0; j < numLocalFe; ++j) {
114
115
116
117
118
              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
          for (std::size_t i = 0; i < numLocalFe; ++i) {
121
            const int local_i = rowNode.localIndex(i);
122
123
            grad_i = A * gradients[i];
            grad_i+= b * shapeValues[i];
124
            for (std::size_t j = 0; j < numLocalFe; ++j) {
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
              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 LocalBasisType = typename FiniteElementType_t<Node>::Traits::LocalBasisType;
147
148
149
150

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

      auto const& localFE = node.finiteElement();
151
      std::size_t numLocalFe = localFE.size();
152
153
154

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

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

158
159
      LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
      auto const& shapeValuesCache = cache.evaluateFunctionAtQp(context, quad);
160

161
162
163
      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
        // Position of the current quadrature point in the reference element
        decltype(auto) local = context.local(quad[iq].position());
164
165

        // the values of the shape functions on the reference element at the quadrature point
166
        auto const& shapeValues = shapeValuesCache[iq];
167

168
169
        // The multiplicative factor in the integral transformation formula
        const auto factor = localFctF(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
170

171
        for (std::size_t i = 0; i < numLocalFe; ++i) {
172
          const int local_i = node.localIndex(i);
173
          elementVector[local_i] += shapeValues[i] * factor;
174
175
176
177
178
179
180
181
182
183
184
185
186
187
        }
      }

      localFctF.unbind();
    }

  private:

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

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

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

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

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

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

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

  /** @} **/

} // end namespace AMDiS