GridFunctionOperator.hpp 14.2 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
      , 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