GridFunctionOperator.hpp 13.4 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
9
10
11
#include <dune/amdis/utility/GetDegree.hpp>

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
113
114
    {
      assert( bound_ );

      int psiDegree = getPolynomialDegree(node);

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
136
137
138
    {
      assert( bound_ );

      int psiDegree = getPolynomialDegree(rowNode);
      int phiDegree = getPolynomialDegree(colNode);

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
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
    {}

    /// Assemble a local element vector on the element that is bound.
    template <class Context, class QuadratureRule,
              class ElementVector, class Node>
    void calculateElementVector(Context const&,         //< container for context related data
                                QuadratureRule const&,  //< a quadrature formula
                                ElementVector&,         //< the output element vector
                                Node const&             //< the node of the test-basis
                                )
    {
      error_exit("Needs to be implemented!");
    }

    /// Assemble a local element matrix on the element that is bound.
    template <class Context, class QuadratureRule,
              class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
    void calculateElementMatrix(Context const&,         //< container for context related data
                                QuadratureRule const&,  //< a quadrature formula
                                ElementMatrix&,         //< the output element matrix
                                RowNode const&,         //< the node of the test-basis
                                ColNode const&,         //< the node of the trial-basis
                                std::integral_constant<bool, sameFE>,   //< finiteelement is the same for test and trial
                                std::integral_constant<bool, sameNode>  //< nodes are the same in the basis-tree
                                )
    {
      error_exit("Needs to be implemented!");
    }
  };

204
205

#ifndef DOXYGEN
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
  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


224
225
226
227
228
229
230
231
232
233
234
235
  namespace tag
  {
    struct deduce {};
    struct order {};
    struct rule {};

  } // end namespace tag


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

236
  template <class Tag, class Expr>
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
  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>
254
255
256
  {
    Tag tag;
    Expr expr;
257
    Rule const& rule;
258
  };
259
#endif
260

261
262

  /// Store tag and expression to create a \ref GridFunctionOperator
263
264
265
  template <class Tag, class Expr>
  auto makeOperator(Tag t, Expr const& expr)
  {
266
267
    using RawExpr = Underlying_t<Expr>;
    static_assert(Concepts::HasGridFunctionOrder<RawExpr>,
268
      "Polynomial degree of expression can not be deduced. You need to provide an explicit value for polynomial order or a quadrature rule in 'makeOperator()'.");
269

270
    return ExpressionPreOperator<Tag, Expr, tag::deduce>{t, expr};
271
272
  }

273
  /// Store tag and expression and polynomial order of expression to create a \ref GridFunctionOperator
274
  template <class Tag, class Expr>
275
276
277
278
279
280
281
282
283
  auto makeOperator(Tag t, Expr const& expr, int order,
                    Dune::QuadratureType::Enum qt = Dune::QuadratureType::GaussLegendre)
  {
    return ExpressionPreOperator<Tag, Expr, tag::order>{t, expr, order, qt};
  }

  /// Store tag and expression and a quadrature rule to create a \ref GridFunctionOperator
  template <class Tag, class Expr, class ctype, int dim>
  auto makeOperator(Tag t, Expr const& expr, Dune::QuadratureRule<ctype,dim> const& rule)
284
  {
285
    return ExpressionPreOperator<Tag, Expr, tag::rule, Dune::QuadratureRule<ctype,dim>>{t, expr, rule};
286
287
  }

288
289
290
291
292
293
294
295
296
297
298
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
  /** @} **/

#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


330
  /// Generate an \ref GridFunctionOperator from a PreOperator (tag, expr).
331
  /// @{
332
333
  template <class Tag, class... Args, class GridView>
  auto makeGridOperator(ExpressionPreOperator<Tag, Args...> const& op, GridView const& gridView)
334
  {
335
    auto gridFct = makeGridFunction(op.expr, gridView);
336
337
338
    auto quadCreator = Impl::makeQuadCreator(op, gridView);
    using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
    return GridFctOp{op.tag, gridFct, quadCreator};
339
340
  }

341
342
  template <class Tag, class... Args, class GridView>
  auto makeGridOperator(std::reference_wrapper<ExpressionPreOperator<Tag, Args...>> op, GridView const& gridView)
343
  {
344
    ExpressionPreOperator<Tag, Args...> const& op_ref = op;
345
    auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView);
346
347
348
    auto quadCreator = Impl::makeQuadCreator(op_ref, gridView);
    using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
    return GridFctOp{op_ref.tag, gridFct, quadCreator};
349
  }
350
351
  /// @}

352

353
  /// Generate a shared_ptr to \ref GridFunctionOperator from a PreOperator (tag, expr).
354
  /// @{
355
356
  template <class Tag, class... Args, class GridView>
  auto makeGridOperatorPtr(ExpressionPreOperator<Tag, Args...> const& op, GridView const& gridView)
357
  {
358
    auto gridFct = makeGridFunction(op.expr, gridView);
359
360
361
    auto quadCreator = Impl::makeQuadCreator(op, gridView);
    using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
    return std::make_shared<GridFctOp>(op.tag, gridFct, quadCreator);
362
363
  }

364
365
  template <class Tag, class... Args, class GridView>
  auto makeGridOperatorPtr(std::reference_wrapper<ExpressionPreOperator<Tag, Args...>> op, GridView const& gridView)
366
  {
367
    ExpressionPreOperator<Tag, Args...> const& op_ref = op;
368
    auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView);
369
370
371
    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);
372
  }
373
  /// @}
374

375
376
#endif // DOXYGEN

377
} // end namespace AMDiS