LocalOperator.hpp 7.05 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#pragma once

#include <cassert>
#include <type_traits>

#include <amdis/GridFunctions.hpp>
#include <amdis/common/Utility.hpp>
#include <amdis/utility/FiniteElementType.hpp>

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

  /// \brief The main implementation of an operator to be used in a \ref LocalAssembler.
  /**
19
20
   * The CRTP Base class for local operators.
   *
21
   * \tparam Derived           The class that derives from this base class
22
   * \tparam LocalContextType  The type of the element or intersection the operator is evaluated on
23
   **/
24
  template <class Derived, class LocalContextType>
25
26
27
  class LocalOperator
  {
  public:
28
    /// The element or intersection the operator is assembled on
29
    using LocalContext = LocalContextType;
30
    /// The codim=0 grid entity
31
    using Element = typename Impl::ContextTypes<LocalContext>::Entity;
32
    /// The geometry of the \ref Element
33
    using Geometry = typename Element::Geometry;
34
35
36
37
38
39
40
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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115

    /// 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.
    /**
     * \param context Container for geometry related data, \ref ContextGeometry
     * \param rowNode The node of the test-basis
     * \param colNode The node of the trial-basis
     * \param elementMatrix The output element matrix
     **/
    template <class Context, class RowNode, class ColNode, class ElementMatrix>
    void calculateElementMatrix(Context const& context,
                                RowNode const& rowNode,
                                ColNode const& colNode,
                                ElementMatrix& elementMatrix)
    {
      derived().getElementMatrix(context, rowNode, colNode, elementMatrix);
    }

    /// \brief Assemble a local element vector on the element that is bound.
    /**
     * \param context Container for geometry related data \ref ContextGeometry
     * \param node The node of the test-basis
     * \param elementVector The output element vector
     **/
    template <class Context, class Node, class ElementVector>
    void calculateElementVector(Context const& context,
                                Node const& node,
                                ElementVector& elementVector)
    {
      derived().getElementVector(context, node, elementVector);
    }


    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() {}

116
    // default implementation called by \ref calculateElementMatrix
117
    template <class Context, class RowNode, class ColNode, class ElementMatrix>
118
119
120
121
    void getElementMatrix(Context const& /*context*/,
                          RowNode const& /*rowNode*/,
                          ColNode const& /*colNode*/,
                          ElementMatrix& /*elementMatrix*/)
122
123
124
125
    {
      error_exit("Needs to be implemented by derived class!");
    }

126
    // default implementation called by \ref calculateElementVector
127
    template <class Context, class Node, class ElementVector>
128
129
130
    void getElementVector(Context const& /*context*/,
                          Node const& /*node*/,
                          ElementVector& /*elementVector*/)
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
    {
      error_exit("Needs to be implemented by derived class!");
    }

    /// \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>
    int getDegree(int order, int coeffDegree, RowNode const& rowNode, ColNode const& colNode) const
    {
      assert( bound_ );
      test_warning(coeffDegree >= 0,
        "polynomial order of coefficient function not determined. Use order 4 by default.");

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

      int degree = psiDegree + phiDegree + (coeffDegree >= 0 ? coeffDegree : 4);
      if (isSimplex_)
        degree -= order;
      if (isAffine_)
        degree += 1; // oder += (order+1)

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

      int psiDegree = polynomialDegree(node);

      int degree = psiDegree + (coeffDegree >= 0 ? coeffDegree : 4);
      if (isSimplex_)
        degree -= order;
      if (isAffine_)
        degree += 1; // oder += (order+1)

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


192
  /// Generate a \ref LocalOperator from a PreOperator.
193
194
  template <class Derived, class LocalContext, class GridView>
  auto makeLocalOperator(LocalOperator<Derived, LocalContext> const& localOp, GridView const& /*gridView*/)
195
196
197
198
  {
    return localOp.derived();
  }

199
  /// Generate a shared_ptr to a \ref LocalOperator from a PreOperator.
200
201
  template <class Derived, class LocalContext, class GridView>
  auto makeLocalOperatorPtr(LocalOperator<Derived, LocalContext> const& localOp, GridView const& /*gridView*/)
202
  {
203
    return std::make_unique<Derived>(localOp.derived());
204
205
206
  }

} // end namespace AMDiS