GridFunctionOperator.hpp 9.07 KB
Newer Older
1
2
3
#pragma once

#include <cassert>
4
#include <type_traits>
5

6
#include <amdis/GridFunctions.hpp>
7
#include <amdis/LocalOperator.hpp>
8
#include <amdis/common/Transposed.hpp>
9
10
11
#include <amdis/common/TypeTraits.hpp>
#include <amdis/typetree/FiniteElementType.hpp>
#include <amdis/utility/QuadratureFactory.hpp>
12
13
14

namespace AMDiS
{
15
16
17
18
19
  /**
   * \addtogroup operators
   * @{
   **/

20
  /// \brief The main implementation of an CRTP-base class for operators using a grid-function
21
  /// coefficient to be used in an \ref Assembler.
22
23
24
25
26
27
28
  /**
   * An Operator that takes a GridFunction as coefficient.
   * Provides quadrature rules and handles the evaluation of the GridFunction at
   * local coordinates.
   *
   * The class is specialized, by deriving from it, in \ref GridFunctionOperator.
   *
29
30
31
32
   * \tparam Derived  The class derived from GridFunctionOperatorBase
   * \tparam LC       The Element or Intersection type
   * \tparam GF       The GridFunction, a LocalFunction is created from, and
   *                  that is evaluated at quadrature points.
33
34
   *
   * **Requirements:**
35
36
   * - `LC` models either Entity (of codim 0) or Intersection
   * - `GF` models the \ref Concepts::GridFunction
37
   **/
38
  template <class Derived, class LC, class GF>
39
  class GridFunctionOperatorBase
40
      : public LocalOperator<Derived, LC>
41
  {
42
    template <class, class>
43
44
    friend class LocalOperator;

45
46
    using ContextType = Impl::ContextTypes<LC>;
    using Super = LocalOperator<Derived, LC>;
47
48

  private:
49
50
    using GridFunction = GF;

51
    /// The type of the localFunction associated with the GridFunction
52
    using LocalFunction = decltype(localFunction(std::declval<GF>()));
53

54
    /// The Codim=0 entity of the grid, the localFunction can be bound to
55
56
    using Element = typename ContextType::Entity;

57
    /// The geometry-type of the grid element
58
59
    using Geometry = typename Element::Geometry;

60
    /// The type of the local coordinates in the \ref Element
61
    using LocalCoordinate = typename GF::EntitySet::LocalCoordinate;
62

63
    /// A factory for QuadratureRules that incooperate the order of the LocalFunction
64
    using QuadFactory = QuadratureFactory<typename Geometry::ctype, LC::mydimension, LocalFunction>;
65

66
  public:
67
    /// \brief Constructor. Stores a copy of `gridFct`.
68
    /**
69
70
71
     * A GridFunctionOperator takes a gridFunction and the
     * differentiation order of the operator, to calculate the
     * quadrature degree in \ref getDegree.
72
     **/
73
74
    template <class GridFct>
    GridFunctionOperatorBase(GridFct&& gridFct, int termOrder)
75
      : gridFct_(FWD(gridFct))
76
      , termOrder_(termOrder)
77
78
    {}

79
    /// Create a quadrature factory from a PreQuadratureFactory, e.g. class derived from \ref QuadratureFactory
80
81
82
    /// \tparam  PQF  A PreQuadratureFactory
    template <class PQF>
    void setQuadFactory(PQF&& pre)
83
    {
84
      using ctype = typename Geometry::ctype;
85
      quadFactory_.emplace(
86
        makeQuadratureFactory<ctype, LC::mydimension, LocalFunction>(FWD(pre)));
87
88
89
    }

  protected:
90
    /// Return expression value at LocalCoordinates
91
    auto coefficient(LocalCoordinate const& local) const
92
    {
93
      assert( this->bound_ );
94
      return (*localFct_)(local);
95
96
    }

97
    /// Create a quadrature rule using the \ref QuadratureFactory by calculating the
98
99
    /// quadrature order necessary to integrate the (bi)linear-form accurately.
    template <class... Nodes>
100
    auto const& getQuadratureRule(Dune::GeometryType type, Nodes const&... nodes) const
101
    {
102
103
104
      assert( bool(quadFactory_) );
      int quadDegree = this->getDegree(termOrder_, quadFactory_->order(), nodes...);
      return quadFactory_->rule(type, quadDegree);
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
  private:
    /// \brief Binds operator to `element` and `geometry`.
    /**
     * Binding an operator to the currently visited element in grid traversal.
     * Since all operators need geometry information, the `Element::Geometry`
     * object `geometry` is created once, during grid traversal, and passed in.
     *
     * By default, it binds the \ref localFct_ to the `element` and the Quadrature
     * factory to the localFunction.
     **/
    void bind_impl(Element const& element, Geometry const& geometry)
    {
      assert( bool(quadFactory_) );
      localFct_.emplace(localFunction(gridFct_));
      localFct_->bind(element);
      quadFactory_->bind(localFct_.value());
    }

    /// Unbinds operator from element.
    void unbind_impl()
    {
      localFct_->unbind();
      localFct_.reset();
    }

132
  private:
133
    /// The gridFunction to be used within the operator
134
    GridFunction gridFct_;
135
136

    /// localFunction associated with gridFunction. Mus be updated whenever gridFunction changes.
137
    Dune::Std::optional<LocalFunction> localFct_;
138

139
    /// Assign each element type a quadrature rule
140
    Dune::Std::optional<QuadFactory> quadFactory_;
141

142
143
    /// the derivative order of this operator (in {0, 1, 2})
    int termOrder_ = 0;
144
145
146
  };


147
148
  /// \brief The transposed operator, implemented in term of its transposed by
  /// calling \ref getElementMatrix with inverted arguments.
149
150
  template <class Derived, class Transposed>
  class GridFunctionOperatorTransposed
151
      : public LocalOperator<Derived, typename Transposed::LocalContext>
152
  {
153
    template <class, class>
154
155
    friend class LocalOperator;

156
157
158
    template <class T, class... Args>
    using Constructable = decltype( new T(std::declval<Args>()...) );

159
  public:
160
161
162
    template <class... Args,
      std::enable_if_t<Dune::Std::is_detected<Constructable, Transposed, Args...>::value, int> = 0>
    GridFunctionOperatorTransposed(Args&&... args)
163
      : transposedOp_(FWD(args)...)
164
165
    {}

166
    /// Redirects the setQuadFactory call top the transposed operator
167
168
169
    /// \tparam  PQF  A PreQuadratureFactory
    template <class PQF>
    void setQuadFactory(PQF&& pre)
170
    {
171
      transposedOp_.setQuadFactory(FWD(pre));
172
173
174
    }

  private:
175
    /// Redirects the bind call top the transposed operator
176
177
    template <class Element, class Geometry>
    void bind_impl(Element const& element, Geometry const& geometry)
178
    {
179
      transposedOp_.bind(element, geometry);
180
181
    }

182
    /// Redirects the unbind call top the transposed operator
183
    void unbind_impl()
184
    {
185
      transposedOp_.unbind();
186
187
    }

188
    /// Apply the assembling to the transposed elementMatrix with interchanged row-/colNode
189
190
191
192
193
194
195
196
    /**
     * \tparam CG  ContextGeometry
     * \tparam RN  RowNode
     * \tparam CN  ColNode
     * \tparam Mat ElementMatrix
     **/
    template <class CG, class RN, class CN, class Mat>
    void getElementMatrix(CG const& contextGeometry, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
197
    {
198
      auto elementMatrixTransposed = transposed(elementMatrix);
199
      transposedOp_.getElementMatrix(
200
        contextGeometry, colNode, rowNode, elementMatrixTransposed);
201
    }
202

203
204
205
  private:
    Transposed transposedOp_;
  };
206

207

208
#ifndef DOXYGEN
209
  template <class Tag, class PreGridFct, class PQF>
210
  struct PreGridFunctionOperator
211
212
  {
    Tag tag;
213
    PreGridFct expr;
214
    PQF quadFactory;
215
  };
216
#endif
217

218
219
220
  /// Store tag and expression into a \ref PreGridFunctionOperator to create a \ref GridFunctionOperator
  template <class Tag, class Expr, class... QuadratureArgs>
  auto makeOperator(Tag tag, Expr&& expr, QuadratureArgs&&... args)
221
  {
222
223
    auto pqf = makePreQuadratureFactory(FWD(args)...);
    return PreGridFunctionOperator<Tag, TYPEOF(expr), TYPEOF(pqf)>{tag, FWD(expr), std::move(pqf)};
224
225
  }

226
227
228
229
  /** @} **/

#ifndef DOXYGEN

230
231
232
233
234
235
236
237
  /// The base-template for GridFunctionOperators
  /**
   * An operator can specialize this class, by deriving from \ref GridFunctionOperatorBase.
   * With the generic function \ref makeLocalOperator, an instance is created. To
   * distinguisch different GridFunction operators, a tag can be provided that has no
   * other effect.
   *
   * \tparam Tag  An Identifier for this GridFunctionOperator
238
   * \tparam LC  An Element or Intersection the operator is evaluated on
239
240
   * \tparam GridFct  A GridFunction evaluated in local coordinates on the bound element
   **/
241
  template <class Tag, class LC, class GridFct>
242
  class GridFunctionOperator
243
      : public GridFunctionOperatorBase<GridFunctionOperator<Tag, LC, GridFct>, LC, GridFct>
244
245
  {};

246
247
  template <class Context, class Tag, class GF, class QF>
  auto makeGridFunctionOperator(Tag tag, GF&& gf, QF&& qf)
248
  {
249
250
    GridFunctionOperator<Tag, Context, TYPEOF(gf)> gfo{tag, FWD(gf)};
    gfo.setQuadFactory(FWD(qf));
251
    return gfo;
252
253
  }

254

255
  /// Generate an \ref GridFunctionOperator from a PreOperator (tag, expr).
256
  /// @{
257
258
  template <class Context, class... Args, class GridView>
  auto makeLocalOperator(PreGridFunctionOperator<Args...> op, GridView const& gridView)
259
  {
260
261
    auto gf = makeGridFunction(std::move(op.expr), gridView);
    return makeGridFunctionOperator<Context>(op.tag, std::move(gf), std::move(op.quadFactory));
262
263
  }

264
265
  template <class Context, class... Args, class GridView>
  auto makeLocalOperator(std::reference_wrapper<PreGridFunctionOperator<Args...>> op_ref, GridView const& gridView)
266
  {
267
268
269
    PreGridFunctionOperator<Args...> const& op = op_ref;
    auto gf = makeGridFunction(std::ref(op.expr), gridView);
    return makeGridFunctionOperator<Context>(op.tag, std::move(gf), op.quadFactory);
270
  }
271
  /// @}
272

273
274
#endif // DOXYGEN

275
} // end namespace AMDiS