LocalOperator.hpp 6.01 KB
Newer Older
1
2
3
4
5
#pragma once

#include <cassert>
#include <type_traits>

6
7
// #include <amdis/GridFunctions.hpp>
#include <amdis/ContextGeometry.hpp>
8
#include <amdis/Output.hpp>
9
#include <amdis/common/TypeTraits.hpp>
10
#include <amdis/functions/Order.hpp>
11
#include <amdis/typetree/FiniteElementType.hpp>
12
13
14
15
16
17
18
19

namespace AMDiS
{
  /**
   * \addtogroup operators
   * @{
   **/

20
  /// \brief The main implementation of an operator to be used in a \ref Assembler.
21
  /**
22
23
   * The CRTP Base class for local operators.
   *
24
25
   * \tparam Derived  The class that derives from this base class
   * \tparam LC       The type of the element or intersection the operator is evaluated on
26
   **/
27
  template <class Derived, class LC>
28
29
  class LocalOperator
  {
30
31
    using ContextType = Impl::ContextTypes<LC>;

32
  public:
33
    /// The element or intersection the operator is assembled on
34
35
    using LocalContext = LC;

36
    /// The codim=0 grid entity
37
38
    using Element = typename ContextType::Entity;

39
    /// The geometry of the \ref Element
40
    using Geometry = typename Element::Geometry;
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76

    /// Initialize the local operator on the current gridView
    template <class GridView>
    void init(GridView const& gridView)
    {
      derived().init_impl(gridView);
    }

    /// \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.
     *
     * By default, it binds the \ref localFct_ to the `element`.
     **/
    void bind(Element const& element, Geometry const& geometry)
    {
      if (!bound_) {
        isSimplex_ = geometry.type().isSimplex();
        isAffine_ = geometry.affine();
        bound_ = true;
      }

      derived().bind_impl(element, geometry);
    }

    /// Unbinds operator from element.
    void unbind()
    {
      derived().unbind_impl();
      bound_ = false;
    }

    /// \brief Assemble a local element matrix on the element that is bound.
    /**
77
     * \param contextGeo Container for geometry related data, \ref ContextGeometry
78
79
80
81
     * \param rowNode The node of the test-basis
     * \param colNode The node of the trial-basis
     * \param elementMatrix The output element matrix
     **/
82
83
    template <class CG, class RN, class CN, class Mat>
    void calculateElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
84
    {
85
      assert( bound_ );
86
      derived().getElementMatrix(contextGeo, rowNode, colNode, elementMatrix);
87
88
89
90
    }

    /// \brief Assemble a local element vector on the element that is bound.
    /**
91
     * \param contextGeo Container for geometry related data \ref ContextGeometry
92
93
94
     * \param node The node of the test-basis
     * \param elementVector The output element vector
     **/
95
96
    template <class CG, class Node, class Vec>
    void calculateElementVector(CG const& contextGeo, Node const& node, Vec& elementVector)
97
    {
98
      assert( bound_ );
99
      derived().getElementVector(contextGeo, node, elementVector);
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
    }


    Derived      & derived()       { return static_cast<Derived&>(*this); }
    Derived const& derived() const { return static_cast<Derived const&>(*this); }


  protected:

    // default implementation. Can be overridden in the derived classes
    template <class GridView>
    void init_impl(GridView const& /*gridView*/) {}

    // default implementation. Can be overridden in the derived classes
    template <class Element, class Geometry>
    void bind_impl(Element const& /*element*/, Geometry const& /*geometry*/) {}

    // default implementation. Can be overridden in the derived classes
    void unbind_impl() {}

    /// \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.
125
126
127
     *
     * \tparam RN  RowNode
     * \tparam CN  ColNode
128
     **/
129
    template <class RN, class CN>
130
    int getDegree(int derivOrder, int coeffDegree, RN const& rowNode, CN const& colNode) const
131
132
133
134
135
    {
      assert( bound_ );
      test_warning(coeffDegree >= 0,
        "polynomial order of coefficient function not determined. Use order 4 by default.");

136
137
      int psiDegree = order(rowNode);
      int phiDegree = order(colNode);
138
139
140

      int degree = psiDegree + phiDegree + (coeffDegree >= 0 ? coeffDegree : 4);
      if (isSimplex_)
141
        degree -= derivOrder;
142
      if (isAffine_)
143
        degree += 1; // oder += (derivOrder+1)
144
145
146
147
148
149
150
151
152
153
154

      return degree;
    }

    /// \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>
155
    int getDegree(int derivOrder, int coeffDegree, Node const& node) const
156
157
158
159
160
    {
      assert( bound_ );
      test_warning(coeffDegree >= 0,
        "polynomial order of coefficient function not determined. Use order 4 by default.");

161
      int psiDegree = order(node);
162
163
164

      int degree = psiDegree + (coeffDegree >= 0 ? coeffDegree : 4);
      if (isSimplex_)
165
        degree -= derivOrder;
166
      if (isAffine_)
167
        degree += 1; // oder += (derivOrder+1)
168
169
170
171
172
173
174
175
176
177
178
179

      return degree;
    }

  protected:

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


180
  /// Generate a \ref LocalOperator from a PreOperator.
181
182
183
184
185
186
187
  /**
   * \tparam Derived  Implementation of LocalOperator
   * \tparam LC       LocalContext
   * \tparam GV       GridView
   **/
  template <class Derived, class LC, class GV>
  auto makeLocalOperator(LocalOperator<Derived, LC> const& localOp, GV const& /*gridView*/)
188
189
190
191
192
  {
    return localOp.derived();
  }

} // end namespace AMDiS