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

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

6
7
#include <dune/amdis/GridFunctions.hpp>
#include <dune/amdis/common/Utility.hpp>
8
#include <dune/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
59
60
61
62
63
64
      , order_(order)
    {}

    /// \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.
     *
65
     * By default, it binds the \ref localFct_ to the `element`.
66
67
68
69
     **/
    void bind(Element const& element, Geometry const& geometry)
    {
      if (!bound_) {
70
        localFct_.bind(element);
71
72
73
74
75
76
77
78
79
80
        isSimplex_ = geometry.type().isSimplex();
        isAffine_ = geometry.affine();
        bound_ = true;
      }
    }

    /// Unbinds operator from element.
    void unbind()
    {
      bound_ = false;
81
      localFct_.unbind();
82
83
    }

84
    /// Return expression value at LocalCoordinates
85
    auto coefficient(LocalCoordinate const& local) const
86
    {
87
88
      assert( bound_ );
      return localFct_(local);
89
90
91
    }


92
93
94
95
96
97
98
99
100
101
    /// 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:
102
103
104
105
106
107
108
    /// \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>
109
    int getDegree(int coeffDegree, Node const& node) const
110
111
112
    {
      assert( bound_ );

113
      int psiDegree = polynomialDegree(node);
114

115
      int degree = psiDegree + coeffDegree;
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
      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>
132
    int getDegree(int coeffDegree, RowNode const& rowNode, ColNode const& colNode) const
133
134
135
    {
      assert( bound_ );

136
137
      int psiDegree = polynomialDegree(rowNode);
      int phiDegree = polynomialDegree(colNode);
138

139
      int degree = psiDegree + phiDegree + coeffDegree;
140
141
142
143
144
145
146
147
148
      if (isSimplex_)
        degree -= order_;
      if (isAffine_)
        degree += 1; // oder += (order+1)

      return degree;
    }

  private:
149
150
    GridFunction gridFct_;
    LocalFunction localFct_;
151

152
153
    QuadratureCreator quadCreator_; //< a creator to provide a quadrature rule
    int order_;                     //< the derivative order of this operator
154
155
156
157
158
159
160

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


161
  /// A default implementation of an GridFunctionOperator if no specialization is available.
162
163
164
165
166
  /**
   * An operator must implement either \ref calculateElementVector, or
   * \ref calculateElementMatrix, if it is a vector or matrix operator,
   * respectively.
   **/
167
  template <class Tag, class GridFct, class QuadratureCreator>
168
  class GridFunctionOperator
169
      : public GridFunctionOperatorBase<GridFct, QuadratureCreator>
170
171
  {
  public:
172
173
    GridFunctionOperator(Tag, GridFct const& gridFct, QuadratureCreator const& quadCreator)
      : GridFunctionOperatorBase<GridFct, QuadratureCreator>(gridFct, 0, quadCreator)
174
175
    {}

Praetorius, Simon's avatar
Praetorius, Simon committed
176
177
178
179
180
181
182
    /// \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
     **/
183
184
    template <class Context, class QuadratureRule,
              class ElementVector, class Node>
Praetorius, Simon's avatar
Praetorius, Simon committed
185
186
187
188
    void calculateElementVector(Context const& contextGeometry,
                                QuadratureRule const& quad,
                                ElementVector& elementVector,
                                Node const& node)
189
190
191
192
    {
      error_exit("Needs to be implemented!");
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
193
194
195
196
197
198
199
200
201
202
    /// \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
     **/
203
204
    template <class Context, class QuadratureRule,
              class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
Praetorius, Simon's avatar
Praetorius, Simon committed
205
206
207
208
209
210
211
    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)
212
213
214
215
216
    {
      error_exit("Needs to be implemented!");
    }
  };

217
218

#ifndef DOXYGEN
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
  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


237
238
239
240
241
242
243
244
245
246
247
248
  namespace tag
  {
    struct deduce {};
    struct order {};
    struct rule {};

  } // end namespace tag


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

249
  template <class Tag, class Expr>
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
  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>
267
268
269
  {
    Tag tag;
    Expr expr;
270
    Rule const& rule;
271
  };
272
#endif
273

274
275

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

Praetorius, Simon's avatar
Praetorius, Simon committed
283
    return ExpressionPreOperator<Tag, Expr, tag::deduce>{tag, expr};
284
285
  }

286
  /// Store tag and expression and polynomial order of expression to create a \ref GridFunctionOperator
287
  template <class Tag, class Expr>
Praetorius, Simon's avatar
Praetorius, Simon committed
288
  auto makeOperator(Tag tag, Expr const& expr, int order,
289
290
                    Dune::QuadratureType::Enum qt = Dune::QuadratureType::GaussLegendre)
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
291
    return ExpressionPreOperator<Tag, Expr, tag::order>{tag, expr, order, qt};
292
293
294
295
  }

  /// 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
296
  auto makeOperator(Tag tag, Expr const& expr, Dune::QuadratureRule<ctype,dim> const& rule)
297
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
298
    return ExpressionPreOperator<Tag, Expr, tag::rule, Dune::QuadratureRule<ctype,dim>>{tag, expr, rule};
299
300
  }

301
302
303
304
305
306
307
308
309
310
311
312
313
314
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
  /** @} **/

#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


343
  /// Generate an \ref GridFunctionOperator from a PreOperator (tag, expr).
344
  /// @{
345
346
  template <class Tag, class... Args, class GridView>
  auto makeGridOperator(ExpressionPreOperator<Tag, Args...> const& op, GridView const& gridView)
347
  {
348
    auto gridFct = makeGridFunction(op.expr, gridView);
349
350
351
    auto quadCreator = Impl::makeQuadCreator(op, gridView);
    using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
    return GridFctOp{op.tag, gridFct, quadCreator};
352
353
  }

354
355
  template <class Tag, class... Args, class GridView>
  auto makeGridOperator(std::reference_wrapper<ExpressionPreOperator<Tag, Args...>> op, GridView const& gridView)
356
  {
357
    ExpressionPreOperator<Tag, Args...> const& op_ref = op;
358
    auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView);
359
360
361
    auto quadCreator = Impl::makeQuadCreator(op_ref, gridView);
    using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
    return GridFctOp{op_ref.tag, gridFct, quadCreator};
362
  }
363
364
  /// @}

365

366
  /// Generate a shared_ptr to \ref GridFunctionOperator from a PreOperator (tag, expr).
367
  /// @{
368
369
  template <class Tag, class... Args, class GridView>
  auto makeGridOperatorPtr(ExpressionPreOperator<Tag, Args...> const& op, GridView const& gridView)
370
  {
371
    auto gridFct = makeGridFunction(op.expr, gridView);
372
373
374
    auto quadCreator = Impl::makeQuadCreator(op, gridView);
    using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
    return std::make_shared<GridFctOp>(op.tag, gridFct, quadCreator);
375
376
  }

377
378
  template <class Tag, class... Args, class GridView>
  auto makeGridOperatorPtr(std::reference_wrapper<ExpressionPreOperator<Tag, Args...>> op, GridView const& gridView)
379
  {
380
    ExpressionPreOperator<Tag, Args...> const& op_ref = op;
381
    auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView);
382
383
384
    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);
385
  }
386
  /// @}
387

388
389
#endif // DOXYGEN

390
} // end namespace AMDiS