LocalAssembler.hpp 7.47 KB
Newer Older
1
2
#pragma once

3
#include <functional>
4
#include <memory>
5
#include <type_traits>
6

7
8
9
10
11
12
#include <dune/geometry/quadraturerules.hh>

#include <dune/amdis/ContextGeometry.hpp>
#include <dune/amdis/LocalAssemblerBase.hpp>
#include <dune/amdis/common/Concepts.hpp>
#include <dune/amdis/utility/FiniteElementType.hpp>
13
14
15

namespace AMDiS
{
16
17
  /// Implementation of interface \ref LocalAssemblerBase
  template <class LocalContext, class Operator, class... Nodes>
18
  class LocalAssembler
19
      : public LocalAssemblerBase<LocalContext, Nodes...>
20
  {
21
22
23
24
25
26
    using Super = LocalAssemblerBase<LocalContext, Nodes...>;

  private:

    using Element = typename Super::Element;
    using ElementMatrixVector = typename Super::ElementMatrixVector;
27

28
29
    using Geometry = typename Super::Geometry;
    using LocalGeometry = typename Super::LocalGeometry;
30

31
    using Context = ContextGeometry<LocalContext, Geometry, LocalGeometry>;
32
33

  public:
34
35
36
37
38
39
40
41
42
43

    /// Constructor. Stores a copy of operator `op`.
    explicit LocalAssembler(Operator const& op)
      : opStorage_(std::make_unique<Operator>(op))
      , op_(*opStorage_)
    {}

    /// Constructor. Stores the reference to the operator.
    explicit LocalAssembler(std::reference_wrapper<Operator> const& op)
      : op_(op)
44
45
    {}

46
47
48
49
50
51
    /// \brief Implementation of \ref LocalAssemblerBase::bind.
    /**
     * Binds the operator `op_` to the `element` and `geometry` and
     * stores point to the `element` and `geometry`.
     **/
    virtual void bind(Element const& element, Geometry const& geometry) final
52
    {
53
      element_ = &element;
54
      geometry_ = &geometry;
55
      op_.bind(element, geometry);
56
    }
57

58
59
60
61
62
63
    /// \brief Implementation of \ref LocalAssemblerBase::unbind
    /**
     * Unbinds the operator `op_` and sets \ref element_ and \ref geometry_
     * to nullptr.
     **/
    virtual void unbind() final
64
    {
65
      op_.unbind();
66
      geometry_ = nullptr;
67
      element_ = nullptr;
68
69
    }

70
    /// Implementation of \ref LocalAssemblerBase::assemble
71
    /**
72
73
74
75
76
77
78
79
     * Provides a quadrature formula with necessary degree to integrate the operator,
     * stores geometry and localGeometry and a \ref ContextGeometry and calls
     * \ref calculateElementVector or \ref calculateElementMatrix on the
     * vector or matrix operator, respectively.
     *
     * The call to operator passes two optimization flags:
     * 1. RowNode and ColNode implement the same local finiteelement
     * 2. RowNode and ColNode are the same node in the globalBasis tree.
80
     **/
81
82
83
84
    virtual bool assemble(
        LocalContext const& localContext,
        ElementMatrixVector& elementMatrixVector,
        Nodes const&... nodes) final
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
      decltype(auto) localGeometry = getLocalGeometry(localContext);

      using QuadratureRules = Dune::QuadratureRules<typename Geometry::ctype, LocalContext::mydimension>;
      int degree = op_.getDegree(nodes...);
      auto const& quad = QuadratureRules::rule(localGeometry.type(), degree);

      Context data{localContext, getGeometry(), localGeometry};
      assembleImpl(data, quad, elementMatrixVector, nodes...);
      return true;
    }


#ifndef DOXYGEN

  protected: // implementation detail

    // matrix assembling
    template <class QuadratureRule, class ElementMatrix, class RowNode, class ColNode>
    void assembleImpl(Context const& context,
                      QuadratureRule const& quad,
                      ElementMatrix& elementMatrix,
                      RowNode const& rowNode, ColNode const& colNode)
    {
      assembleImpl(context, quad, elementMatrix, rowNode, colNode,
        std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>{});
111
112
    }

113
114
115
116
117
118
119
    // local finite elements are the same for rowNode and colNode
    template <class QuadratureRule, class ElementMatrix, class RowNode, class ColNode>
    void assembleImpl(Context const& context,
                      QuadratureRule const& quad,
                      ElementMatrix& elementMatrix,
                      RowNode const& rowNode, ColNode const& colNode,
                      std::true_type /*sameFE*/)
120
    {
121
122
123
124
125
126
127
128
129
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
      if (rowNode.treeIndex() == colNode.treeIndex())
        op_.calculateElementMatrix(
          context, quad, elementMatrix, rowNode, colNode, std::true_type{}, std::true_type{});
      else
        op_.calculateElementMatrix(
          context, quad, elementMatrix, rowNode, colNode, std::true_type{}, std::false_type{});
    }

    // local finite elements are different for rowNode and colNode
    template <class QuadratureRule, class ElementMatrix, class RowNode, class ColNode>
    void assembleImpl(Context const& context,
                      QuadratureRule const& quad,
                      ElementMatrix& elementMatrix,
                      RowNode const& rowNode, ColNode const& colNode,
                      std::false_type /*sameFE*/)
    {
      op_.calculateElementMatrix(
        context, quad, elementMatrix, rowNode, colNode, std::false_type{}, std::false_type{});
    }


    // vector assembling
    template <class QuadratureRule, class ElementVector, class Node>
    void assembleImpl(Context const& context,
                      QuadratureRule const& quad,
                      ElementVector& elementVector,
                      Node const& node)
    {
      op_.calculateElementVector(context, quad, elementVector, node);
    }

#endif // DOXYGEN


  public:

    Element const& getElement() const
    {
      assert( element_ );
      return *element_;
161
162
163
164
    }

    Geometry const& getGeometry() const
    {
165
166
      assert( geometry_ );
      return *geometry_;
167
168
    }

169
170
171
172
173
174
175
176
177
178
    decltype(auto) getLocalGeometry(LocalContext const& context) const
    {
      return getLocalGeometry(context, std::is_same<LocalContext, Element>{});
    }


  private: // implementation detail

    // localGeometry for Entities
    Geometry const& getLocalGeometry(Element const& element, std::true_type) const
179
    {
180
      return *geometry_;
181
    }
182

183
184
    // localGeometry for intersections
    LocalGeometry getLocalGeometry(LocalContext const& intersection, std::false_type) const
185
    {
186
      return intersection.geometryInInside();
187
188
    }

189

190
191
  private:

192
193
    std::unique_ptr<Operator> opStorage_;     //< the stored operator, implementing \ref ExpressionOperatorBase
    Operator& op_;
194

195
196
    Element const* element_ = nullptr;
    Geometry const* geometry_ = nullptr;
197
198
  };

199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227

  /// Generate a \ref LocalAssembler on a given `LocalContext` (element or intersection)
  template <class LocalContext, class Operator, class... Nodes,
    std::enable_if_t<not traits::is_reference_wrapper<Operator>::value, int> = 0>
  auto makeLocalAssembler(Operator const& op, Nodes const&...)
  {
    return LocalAssembler<LocalContext, Operator, Nodes...>{op};
  }

  template <class LocalContext, class Operator, class... Nodes>
  auto makeLocalAssembler(std::reference_wrapper<Operator> const& op, Nodes const&...)
  {
    return LocalAssembler<LocalContext, Operator, Nodes...>{op};
  }

  /// Generate a shared_ptr to \ref LocalAssembler on a given `LocalContext` (element or intersection)
  template <class LocalContext, class Operator, class... Nodes,
    std::enable_if_t<not traits::is_reference_wrapper<Operator>::value, int> = 0>
  auto makeLocalAssemblerPtr(Operator const& op, Nodes const&...)
  {
    return std::make_shared<LocalAssembler<LocalContext, Operator, Nodes...>>(op);
  }

  template <class LocalContext, class Operator, class... Nodes>
  auto makeLocalAssemblerPtr(std::reference_wrapper<Operator> const& op, Nodes const&...)
  {
    return std::make_shared<LocalAssembler<LocalContext, Operator, Nodes...>>(op);
  }

228
} // end namespace AMDiS