LocalOperator.hpp 5.97 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
10
#include <amdis/common/TypeTraits.hpp>
#include <amdis/typetree/FiniteElementType.hpp>
11
12
13
14
15
16
17
18

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

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

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

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

38
    /// The geometry of the \ref Element
39
    using Geometry = typename Element::Geometry;
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

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

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


    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.
124
125
126
     *
     * \tparam RN  RowNode
     * \tparam CN  ColNode
127
     **/
128
129
    template <class RN, class CN>
    int getDegree(int order, int coeffDegree, RN const& rowNode, CN const& colNode) const
130
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
    {
      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
  };


179
  /// Generate a \ref LocalOperator from a PreOperator.
180
181
182
183
184
185
186
  /**
   * \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*/)
187
188
189
190
191
  {
    return localOp.derived();
  }

} // end namespace AMDiS