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

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

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

namespace AMDiS
{
14
15
16
17
18
19
20
21
22
23
24
25
26
  /**
   * \addtogroup operators
   * @{
   **/

  /// \brief The main implementation of an operator to be used in a \ref LocalAssembler.
  /**
   * 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.
   *
27
   * \tparam LocalContext The Element or Intersection type
28
29
30
31
   * \tparam GridFunction The GridFunction, a LocalFunction is created from, and
   *                      that is evaluated at quadrature points.
   *
   * **Requirements:**
32
   * - `LocalContext` models either Entity (of codim 0) or Intersection
33
34
   * - `GridFunction` models the \ref Concepts::GridFunction
   **/
35
  template <class Derived, class LocalContext, class GridFunction>
36
  class GridFunctionOperatorBase
37
      : public LocalOperator<Derived, LocalContext>
38
  {
39
    using LocalFunction = decltype(localFunction(std::declval<GridFunction>()));
40
41

    using Element = typename Impl::ContextTypes<LocalContext>::Entity;
42
43
44
45
    using Geometry = typename Element::Geometry;

    using LocalCoordinate = typename GridFunction::EntitySet::LocalCoordinate;

46
    using QuadFactory = QuadratureFactory<typename Geometry::ctype, LocalContext::mydimension, LocalFunction>;
47

48
  public:
49
    /// \brief Constructor. Stores a copy of `gridFct`.
50
51
52
53
54
    /**
     * An expression operator takes an expression, following the interface of
     * \ref ExpressionBase, and stores a copy. Additionally, it gets the
     * differentiation order, to calculate the quadrature degree in \ref getDegree.
     **/
55
    GridFunctionOperatorBase(GridFunction const& gridFct, int termOrder)
56
57
      : gridFct_(gridFct)
      , localFct_(localFunction(gridFct_))
58
      , termOrder_(termOrder)
59
60
    {}

61
62
63
64
    // Copy constructor
    GridFunctionOperatorBase(GridFunctionOperatorBase const& that)
      : gridFct_(that.gridFct_)
      , localFct_(localFunction(gridFct_)) // recreate localFct_ on new gridFct_
65
      , termOrder_(that.termOrder_)
66
67
68
69
70
71
    {}

    // Move constructor
    GridFunctionOperatorBase(GridFunctionOperatorBase&& that)
      : gridFct_(std::move(that.gridFct_))
      , localFct_(localFunction(gridFct_)) // recreate localFct_ on new gridFct_
72
      , termOrder_(std::move(that.termOrder_))
73
74
    {}

75
76
77
78
79
80
    /// \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.
     *
81
     * By default, it binds the \ref localFct_ to the `element`.
82
     **/
83
    void bind_impl(Element const& element, Geometry const& geometry)
84
    {
85
      assert( bool(quadFactory_) );
86
87
      localFct_.bind(element);
      quadFactory_->bind(localFct_);
88
89
90
    }

    /// Unbinds operator from element.
91
    void unbind_impl()
92
    {
93
      localFct_.unbind();
94
95
    }

96
97
98
    template <class PreQuadFactory>
    void setQuadFactory(PreQuadFactory const& pre)
    {
99
      quadFactory_ = makeQuadratureFactoryPtr<typename Geometry::ctype, LocalContext::mydimension, LocalFunction>(pre);
100
101
102
    }

  protected:
103
    /// Return expression value at LocalCoordinates
104
    auto coefficient(LocalCoordinate const& local) const
105
    {
106
      assert( this->bound_ );
107
      return localFct_(local);
108
109
    }

110
111
112
    /// Create a quadrature rule using the \ref QuadratureCreator by calculating the
    /// quadrature order necessary to integrate the (bi)linear-form accurately.
    template <class... Nodes>
113
    auto const& getQuadratureRule(Dune::GeometryType type, Nodes const&... nodes) const
114
    {
115
116
117
      assert( bool(quadFactory_) );
      int quadDegree = this->getDegree(termOrder_, quadFactory_->order(), nodes...);
      return quadFactory_->rule(type, quadDegree);
118
119
120
    }

  private:
121
    /// The gridFunction to be used within the operator
122
    GridFunction gridFct_;
123
124

    /// localFunction associated with gridFunction. Mus be updated whenever gridFunction changes.
125
    LocalFunction localFct_;
126

127
128
    /// Assign each element type a quadrature rule
    std::shared_ptr<QuadFactory> quadFactory_ = std::make_shared<QuadFactory>();
129

130
131
    /// the derivative order of this operator (in {0, 1, 2})
    int termOrder_ = 0;
132
133
134
  };


135
136
  /// \brief The transposed operator, implemented in term of its transposed by
  /// calling \ref getElementMatrix with inverted arguments.
137
138
  template <class Derived, class Transposed>
  class GridFunctionOperatorTransposed
139
      : public LocalOperator<Derived, typename Transposed::LocalContext>
140
  {
141
142
143
    template <class T, class... Args>
    using Constructable = decltype( new T(std::declval<Args>()...) );

144
  public:
145
146
147
148
    template <class... Args,
      std::enable_if_t<Dune::Std::is_detected<Constructable, Transposed, Args...>::value, int> = 0>
    GridFunctionOperatorTransposed(Args&&... args)
      : transposedOp_(std::forward<Args>(args)...)
149
150
    {}

151
152
    template <class Element, class Geometry>
    void bind_impl(Element const& element, Geometry const& geometry)
153
    {
154
      transposedOp_.bind_impl(element, geometry);
155
156
    }

157
    void unbind_impl()
158
    {
159
      transposedOp_.unbind_impl();
160
    }
161

162
163
    template <class PreQuadFactory>
    void setQuadFactory(PreQuadFactory const& pre)
164
    {
165
      transposedOp_.setQuadFactory(pre);
166
167
    }

168
169
170
171
172
173
174
175
176
177
    template <class Context, class RowNode, class ColNode, class ElementMatrix>
    void getElementMatrix(Context const& context,
                          RowNode const& rowNode,
                          ColNode const& colNode,
                          ElementMatrix& elementMatrix)
    {
      auto elementMatrixTransposed = trans(elementMatrix);
      transposedOp_.getElementMatrix(
        context, colNode, rowNode, elementMatrixTransposed);
    }
178

179
180
181
  private:
    Transposed transposedOp_;
  };
182

183
184
  template <class Tag, class PreGridFct, class PreQuadFactory>
  struct PreGridFunctionOperator
185
186
  {
    Tag tag;
187
188
    PreGridFct expr;
    PreQuadFactory quadFactory;
189
190
  };

191
  /// Store tag and expression to create a \ref GridFunctionOperator
192
193
  template <class Tag, class PreGridFct, class... QuadratureArgs>
  auto makeOperator(Tag tag, PreGridFct const& expr, QuadratureArgs&&... args)
194
  {
195
196
    auto preQuadFactory = makePreQuadratureFactory(std::forward<QuadratureArgs>(args)...);
    return PreGridFunctionOperator<Tag, PreGridFct, decltype(preQuadFactory)>{tag, expr, preQuadFactory};
197
198
  }

199
200
201
202
  /** @} **/

#ifndef DOXYGEN

203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
  /// 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
   * \tparam LocalContext  An Element or Intersection the operator is evaluated on
   * \tparam GridFct  A GridFunction evaluated in local coordinates on the bound element
   **/
  template <class Tag, class LocalContext, class GridFct>
  class GridFunctionOperator
      : public GridFunctionOperatorBase<GridFunctionOperator<Tag, LocalContext, GridFct>, LocalContext, GridFct>
  {};


220
  /// Generate an \ref GridFunctionOperator from a PreOperator (tag, expr).
221
  /// @{
222
  template <class LocalContext, class Tag, class... Args, class GridView>
223
  auto makeLocalOperator(PreGridFunctionOperator<Tag, Args...> const& op, GridView const& gridView)
224
  {
225
    auto gridFct = makeGridFunction(op.expr, gridView);
226
    using GridFctOp = GridFunctionOperator<Tag, LocalContext, decltype(gridFct)>;
227
228
229
    GridFctOp localOperator{op.tag, gridFct};
    localOperator.setQuadFactory(op.quadFactory);
    return localOperator;
230
231
  }

232
  template <class LocalContext, class Tag, class... Args, class GridView>
233
  auto makeLocalOperator(std::reference_wrapper<PreGridFunctionOperator<Tag, Args...>> op, GridView const& gridView)
234
  {
235
    PreGridFunctionOperator<Tag, Args...> const& op_ref = op;
236
    auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView);
237
    using GridFctOp = GridFunctionOperator<Tag, LocalContext, decltype(gridFct)>;
238
239
240
    GridFctOp localOperator{op_ref.tag, gridFct};
    localOperator.setQuadFactory(op_ref.quadFactory);
    return localOperator;
241
  }
242
243
  /// @}

244

245
  /// Generate a shared_ptr to \ref GridFunctionOperator from a PreOperator (tag, expr).
246
  /// @{
247
  template <class LocalContext, class Tag, class... Args, class GridView>
248
  auto makeLocalOperatorPtr(PreGridFunctionOperator<Tag, Args...> const& op, GridView const& gridView)
249
  {
250
    auto gridFct = makeGridFunction(op.expr, gridView);
251
    using GridFctOp = GridFunctionOperator<Tag, LocalContext, decltype(gridFct)>;
252
253
254
    auto localOperator = std::make_shared<GridFctOp>(op.tag, gridFct);
    localOperator->setQuadFactory(op.quadFactory);
    return localOperator;
255
256
  }

257
  template <class LocalContext, class Tag, class... Args, class GridView>
258
  auto makeLocalOperatorPtr(std::reference_wrapper<PreGridFunctionOperator<Tag, Args...>> op, GridView const& gridView)
259
  {
260
    PreGridFunctionOperator<Tag, Args...> const& op_ref = op;
261
    auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView);
262
    using GridFctOp = GridFunctionOperator<Tag, LocalContext, decltype(gridFct)>;
263
264
265
    auto localOperator = std::make_shared<GridFctOp>(op_ref.tag, gridFct);
    localOperator->setQuadFactory(op_ref.quadFactory);
    return localOperator;
266
  }
267
  /// @}
268

269
270
#endif // DOXYGEN

271
} // end namespace AMDiS