FirstOrderTestGradTrial.hpp 4.05 KB
Newer Older
1 2 3 4
#pragma once

#include <type_traits>

5
#include <amdis/GridFunctionOperator.hpp>
6
#include <amdis/utility/LocalBasisCache.hpp>
7
#include <amdis/Output.hpp>
8
#include <amdis/common/StaticSize.hpp>
9 10 11

namespace AMDiS
{
12 13 14 15 16
  /**
   * \addtogroup operators
   * @{
   **/

17 18 19
  namespace tag
  {
    struct test_gradtrial {};
20
    struct grad_trial {};
21 22 23
  }


24
  /// first-order operator \f$ \langle\psi, \mathbf{b}\cdot\nabla\phi\rangle \f$
25 26 27 28
  template <class LocalContext, class GridFct>
  class GridFunctionOperator<tag::test_gradtrial, LocalContext, GridFct>
      : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_gradtrial, LocalContext, GridFct>,
                                        LocalContext, GridFct>
29
  {
30
    static const int dow = ContextGeometry<LocalContext>::dow;
31
    using Self = GridFunctionOperator;
32
    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
33

34
    static_assert( Size_v<typename GridFct::Range> == dow, "Expression must be of vector type." );
35 36

  public:
37 38
    GridFunctionOperator(tag::test_gradtrial, GridFct const& expr)
      : Super(expr, 1)
39 40
    {}

41 42 43 44 45
    template <class Context, class RowNode, class ColNode, class ElementMatrix>
    void getElementMatrix(Context const& context,
                          RowNode const& rowNode,
                          ColNode const& colNode,
                          ElementMatrix& elementMatrix)
46
    {
47 48
      static_assert(RowNode::isLeaf && ColNode::isLeaf,
        "Operator can be applied to Leaf-Nodes only.");
49

50
      auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
51 52
      auto const& rowLocalFE = rowNode.finiteElement();
      auto const& colLocalFE = colNode.finiteElement();
53 54
      std::size_t rowSize = rowLocalFE.size();
      std::size_t colSize = colLocalFE.size();
55

56
      using RangeFieldType = typename NodeQuadCache<ColNode>::Traits::RangeFieldType;
57

58 59
      NodeQuadCache<RowNode> rowCache(rowLocalFE.localBasis());
      NodeQuadCache<ColNode> colCache(colLocalFE.localBasis());
60

61 62
      auto const& shapeValuesCache = rowCache.evaluateFunctionAtQP(context, quad);
      auto const& shapeGradientsCache = colCache.evaluateJacobianAtQP(context, quad);
63
      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
64
        // Position of the current quadrature point in the reference element
65
        decltype(auto) local = context.local(quad[iq].position());
66 67

        // The transposed inverse Jacobian of the map from the reference element to the element
68
        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
69 70

        // The multiplicative factor in the integral transformation formula
71
        const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
72
        const auto b = Super::coefficient(local);
73

74 75
        // the values of the shape functions on the reference element at the quadrature point
        auto const& shapeValues = shapeValuesCache[iq];
76 77

        // The gradients of the shape functions on the reference element
78
        auto const& shapeGradients = shapeGradientsCache[iq];
79 80

        // Compute the shape function gradients on the real element
81 82
        using WorldVector = FieldVector<RangeFieldType,Context::dow>;
        thread_local std::vector<WorldVector> colGradients;
83
        colGradients.resize(shapeGradients.size());
84

85
        for (std::size_t i = 0; i < colGradients.size(); ++i)
86
          jacobian.mv(shapeGradients[i][0], colGradients[i]);
87

88
        for (std::size_t j = 0; j < colSize; ++j) {
89 90
          const auto local_j = colNode.localIndex(j);
          const auto value = factor * (b * colGradients[j]);
91
          for (std::size_t i = 0; i < rowSize; ++i) {
92 93
            const auto local_i = rowNode.localIndex(i);
            elementMatrix[local_i][local_j] += value * shapeValues[i];
94 95 96 97 98 99 100
          }
        }
      }

    }
  };

101 102 103 104 105 106 107 108

  /// Create a first-order term with derivative on test-function
  template <class Expr, class... QuadratureArgs>
  auto fot(Expr&& expr, tag::grad_trial, QuadratureArgs&&... args)
  {
    return makeOperator(tag::test_gradtrial{}, FWD(expr), FWD(args)...);
  }

109 110
  /** @} **/

111
} // end namespace AMDiS