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

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

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

namespace AMDiS
{
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
  /**
   * \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.
   *
   * \tparam GridFunction The GridFunction, a LocalFunction is created from, and
   *                      that is evaluated at quadrature points.
   * \tparam QuadratureCreator A functor that provides a \ref Dune::QuadratureRule.
   *
   * **Requirements:**
   * - `GridFunction` models the \ref Concepts::GridFunction
   * - `QuadratureCreator` models \ref Concepts::Callable<Dune::GeometryType, LocalFunction, F>
   *   where `F` is a functor of the signature `int(int)` that calculates the
   *   degree of the (bi)linear-form. The argument passed to `F` is the polynomial
   *   order of the GridFunction.
   **/
  template <class GridFunction, class QuadratureCreator>
37
  class GridFunctionOperatorBase
38
  {
39
40
41
42
43
44
    using LocalFunction = decltype(localFunction(std::declval<GridFunction>()));
    using Element = typename GridFunction::EntitySet::Element;
    using Geometry = typename Element::Geometry;

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

45
46
47
48
49
50
51
  public:
    /// \brief Constructor. Stores a copy of `expr`.
    /**
     * 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.
     **/
52
    GridFunctionOperatorBase(GridFunction const& gridFct, QuadratureCreator creator, int order)
53
54
      : gridFct_(gridFct)
      , localFct_(localFunction(gridFct_))
55
      , quadCreator_(creator)
56
57
58
      , order_(order)
    {}

59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
    // Copy constructor
    GridFunctionOperatorBase(GridFunctionOperatorBase const& that)
      : gridFct_(that.gridFct_)
      , localFct_(localFunction(gridFct_)) // recreate localFct_ on new gridFct_
      , quadCreator_(that.quadCreator_)
      , order_(that.order_)
    {}

    // Move constructor
    GridFunctionOperatorBase(GridFunctionOperatorBase&& that)
      : gridFct_(std::move(that.gridFct_))
      , localFct_(localFunction(gridFct_)) // recreate localFct_ on new gridFct_
      , quadCreator_(std::move(that.quadCreator_))
      , order_(std::move(that.order_))
    {}

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
84
85
     **/
    void bind(Element const& element, Geometry const& geometry)
    {
      if (!bound_) {
86
        localFct_.bind(element);
87
88
89
90
91
92
93
94
95
96
        isSimplex_ = geometry.type().isSimplex();
        isAffine_ = geometry.affine();
        bound_ = true;
      }
    }

    /// Unbinds operator from element.
    void unbind()
    {
      bound_ = false;
97
      localFct_.unbind();
98
99
    }

100
    /// Return expression value at LocalCoordinates
101
    auto coefficient(LocalCoordinate const& local) const
102
    {
103
104
      assert( bound_ );
      return localFct_(local);
105
106
107
    }


108
109
110
111
112
113
114
115
116
117
    /// Create a quadrature rule using the \ref QuadratureCreator by calculating the
    /// quadrature order necessary to integrate the (bi)linear-form accurately.
    template <class... Nodes>
    decltype(auto) getQuadratureRule(Dune::GeometryType type, Nodes const&... nodes) const
    {
      auto degreeFunctor = [&,this](int order) -> int { return this->getDegree(order, nodes...); };
      return quadCreator_(type, localFct_, degreeFunctor);
    }

  protected:
118
119
120
121
122
123
124
    /// \brief Return the quadrature degree for a vector operator.
    /**
     * The quadrature degree that is necessary, to integrate the expression
     * multiplied with (possibly derivatives of) shape functions. Thus it needs
     * the order of derivatives, this operator implements.
     **/
    template <class Node>
125
    int getDegree(int coeffDegree, Node const& node) const
126
127
128
    {
      assert( bound_ );

129
      int psiDegree = polynomialDegree(node);
130

131
      int degree = psiDegree + coeffDegree;
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
      if (isSimplex_)
        degree -= order_;
      if (isAffine_)
        degree += 1; // oder += (order+1)

      return degree;
    }


    /// \brief Return the quadrature degree for a matrix operator.
    /**
     * The quadrature degree that is necessary, to integrate the expression
     * multiplied with (possibly derivatives of) shape functions. Thus it needs
     * the order of derivatives, this operator implements.
     **/
    template <class RowNode, class ColNode>
148
    int getDegree(int coeffDegree, RowNode const& rowNode, ColNode const& colNode) const
149
150
151
    {
      assert( bound_ );

152
153
      int psiDegree = polynomialDegree(rowNode);
      int phiDegree = polynomialDegree(colNode);
154

155
      int degree = psiDegree + phiDegree + coeffDegree;
156
157
158
159
160
161
162
163
164
      if (isSimplex_)
        degree -= order_;
      if (isAffine_)
        degree += 1; // oder += (order+1)

      return degree;
    }

  private:
165
166
    GridFunction gridFct_;
    LocalFunction localFct_;
167

168
169
    QuadratureCreator quadCreator_; //< a creator to provide a quadrature rule
    int order_;                     //< the derivative order of this operator
170
171
172
173
174
175
176

    bool isSimplex_ = false;    //< the bound element is a simplex
    bool isAffine_ = false;     //< the bound geometry is affine
    bool bound_ = false;        //< an element is bound to the operator
  };


177
  /// A default implementation of an GridFunctionOperator if no specialization is available.
178
179
180
181
182
  /**
   * An operator must implement either \ref calculateElementVector, or
   * \ref calculateElementMatrix, if it is a vector or matrix operator,
   * respectively.
   **/
183
  template <class Tag, class GridFct, class QuadratureCreator>
184
  class GridFunctionOperator
185
      : public GridFunctionOperatorBase<GridFct, QuadratureCreator>
186
187
  {
  public:
188
189
    GridFunctionOperator(Tag, GridFct const& gridFct, QuadratureCreator const& quadCreator)
      : GridFunctionOperatorBase<GridFct, QuadratureCreator>(gridFct, 0, quadCreator)
190
191
    {}

Praetorius, Simon's avatar
Praetorius, Simon committed
192
193
194
195
196
197
198
    /// \brief Assemble a local element vector on the element that is bound.
    /**
     * \param contextGeometry Container for context related data
     * \param quad A quadrature formula
     * \param elementVector The output element vector
     * \param node The node of the test-basis
     **/
199
200
    template <class Context, class QuadratureRule,
              class ElementVector, class Node>
Praetorius, Simon's avatar
Praetorius, Simon committed
201
202
203
204
    void calculateElementVector(Context const& contextGeometry,
                                QuadratureRule const& quad,
                                ElementVector& elementVector,
                                Node const& node)
205
206
207
208
    {
      error_exit("Needs to be implemented!");
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
209
210
211
212
213
214
215
216
217
218
    /// \brief Assemble a local element matrix on the element that is bound.
    /**
     * \param contextGeometry Container for context related data
     * \param quad A quadrature formula
     * \param elementMatrix The output element matrix
     * \param rowNode The node of the test-basis
     * \param colNode The node of the trial-basis
     * \param flag1 finiteelement is the same for test and trial
     * \param flag2 nodes are the same in the basis-tree
     **/
219
220
    template <class Context, class QuadratureRule,
              class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
Praetorius, Simon's avatar
Praetorius, Simon committed
221
222
223
224
225
226
227
    void calculateElementMatrix(Context const& contextGeometry,
                                QuadratureRule const& quad,
                                ElementMatrix& elementMatrix,
                                RowNode const& rowNode,
                                ColNode const& colNode,
                                std::integral_constant<bool, sameFE> flag1,
                                std::integral_constant<bool, sameNode> flag2)
228
229
230
231
232
    {
      error_exit("Needs to be implemented!");
    }
  };

233
234

#ifndef DOXYGEN
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
  namespace Concepts
  {
    namespace Definition
    {
      struct HasGridFunctionOrder
      {
        template <class F, class GridView>
        auto requires_(F&& f, GridView const& gridView)
          -> decltype( order(localFunction(makeGridFunction(f, gridView))) );
      };
    }

    template <class F, class GridView = Dune::YaspGrid<2>::LeafGridView>
    constexpr bool HasGridFunctionOrder = models<Definition::HasGridFunctionOrder(F, GridView)>;

  } // end namespace Concepts


253
254
255
256
257
258
259
260
261
262
263
264
  namespace tag
  {
    struct deduce {};
    struct order {};
    struct rule {};

  } // end namespace tag


  template <class Tag, class Expr, class QuadTag, class Rule = void>
  struct ExpressionPreOperator;

265
  template <class Tag, class Expr>
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
  struct ExpressionPreOperator<Tag, Expr, tag::deduce>
  {
    Tag tag;
    Expr expr;
  };

  template <class Tag, class Expr>
  struct ExpressionPreOperator<Tag, Expr, tag::order>
  {
    Tag tag;
    Expr expr;
    int order;
    Dune::QuadratureType::Enum qt;
  };

  template <class Tag, class Expr, class Rule>
  struct ExpressionPreOperator<Tag, Expr, tag::rule, Rule>
283
284
285
  {
    Tag tag;
    Expr expr;
286
    Rule const& rule;
287
  };
288
#endif
289

290
291

  /// Store tag and expression to create a \ref GridFunctionOperator
292
  template <class Tag, class Expr>
Praetorius, Simon's avatar
Praetorius, Simon committed
293
  auto makeOperator(Tag tag, Expr const& expr)
294
  {
295
296
    using RawExpr = Underlying_t<Expr>;
    static_assert(Concepts::HasGridFunctionOrder<RawExpr>,
297
      "Polynomial degree of expression can not be deduced. You need to provide an explicit value for polynomial order or a quadrature rule in 'makeOperator()'.");
298

Praetorius, Simon's avatar
Praetorius, Simon committed
299
    return ExpressionPreOperator<Tag, Expr, tag::deduce>{tag, expr};
300
301
  }

302
  /// Store tag and expression and polynomial order of expression to create a \ref GridFunctionOperator
303
  template <class Tag, class Expr>
Praetorius, Simon's avatar
Praetorius, Simon committed
304
  auto makeOperator(Tag tag, Expr const& expr, int order,
305
306
                    Dune::QuadratureType::Enum qt = Dune::QuadratureType::GaussLegendre)
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
307
    return ExpressionPreOperator<Tag, Expr, tag::order>{tag, expr, order, qt};
308
309
310
311
  }

  /// Store tag and expression and a quadrature rule to create a \ref GridFunctionOperator
  template <class Tag, class Expr, class ctype, int dim>
Praetorius, Simon's avatar
Praetorius, Simon committed
312
  auto makeOperator(Tag tag, Expr const& expr, Dune::QuadratureRule<ctype,dim> const& rule)
313
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
314
    return ExpressionPreOperator<Tag, Expr, tag::rule, Dune::QuadratureRule<ctype,dim>>{tag, expr, rule};
315
316
  }

317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
  /** @} **/

#ifndef DOXYGEN

  namespace Impl
  {
    // Deduce polynomial order of expression automatically. Create standard quadrature rule with degree
    // build up of operator derivative order, and polynomial degrees of trial/test functions.
    template <class Tag, class Expr, class GridView>
    auto makeQuadCreator(ExpressionPreOperator<Tag, Expr, tag::deduce> const& /*op*/, GridView const& /*gridView*/)
    {
      using QuadratureRules = Dune::QuadratureRules<typename GridView::ctype, GridView::dimension>;
      return [](auto type, auto&& localFct, auto getDegree) -> auto const&
      {
        return QuadratureRules::rule(type, getDegree(order(localFct)));
      };
    }

    // Provide polynomial order of expression explicitly
    template <class Tag, class Expr, class GridView>
    auto makeQuadCreator(ExpressionPreOperator<Tag, Expr, tag::order> const& op, GridView const& /*gridView*/)
    {
      using QuadratureRules = Dune::QuadratureRules<typename GridView::ctype, GridView::dimension>;
      return [order=op.order,qt=op.qt](auto type, auto&& /*localFct*/, auto getDegree) -> auto const&
      {
        return QuadratureRules::rule(type, getDegree(order), qt);
      };
    }

    // Provide quadrature rule explicitly
    template <class Tag, class Expr, class Rule, class GridView>
    auto makeQuadCreator(ExpressionPreOperator<Tag, Expr, tag::rule, Rule> const& op, GridView const& /*gridView*/)
    {
      return [&rule=op.rule](auto&&...) -> auto const&
      {
        return rule;
      };
    }

  } // end namespace Impl


359
  /// Generate an \ref GridFunctionOperator from a PreOperator (tag, expr).
360
  /// @{
361
362
  template <class Tag, class... Args, class GridView>
  auto makeGridOperator(ExpressionPreOperator<Tag, Args...> const& op, GridView const& gridView)
363
  {
364
    auto gridFct = makeGridFunction(op.expr, gridView);
365
366
367
    auto quadCreator = Impl::makeQuadCreator(op, gridView);
    using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
    return GridFctOp{op.tag, gridFct, quadCreator};
368
369
  }

370
371
  template <class Tag, class... Args, class GridView>
  auto makeGridOperator(std::reference_wrapper<ExpressionPreOperator<Tag, Args...>> op, GridView const& gridView)
372
  {
373
    ExpressionPreOperator<Tag, Args...> const& op_ref = op;
374
    auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView);
375
376
377
    auto quadCreator = Impl::makeQuadCreator(op_ref, gridView);
    using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
    return GridFctOp{op_ref.tag, gridFct, quadCreator};
378
  }
379
380
  /// @}

381

382
  /// Generate a shared_ptr to \ref GridFunctionOperator from a PreOperator (tag, expr).
383
  /// @{
384
385
  template <class Tag, class... Args, class GridView>
  auto makeGridOperatorPtr(ExpressionPreOperator<Tag, Args...> const& op, GridView const& gridView)
386
  {
387
    auto gridFct = makeGridFunction(op.expr, gridView);
388
389
390
    auto quadCreator = Impl::makeQuadCreator(op, gridView);
    using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
    return std::make_shared<GridFctOp>(op.tag, gridFct, quadCreator);
391
392
  }

393
394
  template <class Tag, class... Args, class GridView>
  auto makeGridOperatorPtr(std::reference_wrapper<ExpressionPreOperator<Tag, Args...>> op, GridView const& gridView)
395
  {
396
    ExpressionPreOperator<Tag, Args...> const& op_ref = op;
397
    auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView);
398
399
400
    auto quadCreator = Impl::makeQuadCreator(op_ref, gridView);
    using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
    return std::make_shared<GridFctOp>(op_ref.tag, gridFct, quadCreator);
401
  }
402
  /// @}
403

404
405
#endif // DOXYGEN

406
} // end namespace AMDiS