Assembler.hpp 3.51 KB
Newer Older
1
2
3
#pragma once

#include <memory>
4
#include <tuple>
5

6
7
8
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>

9
#include <dune/amdis/DirichletBC.hpp>
10
#include <dune/amdis/LinearAlgebra.hpp>
11
#include <dune/amdis/Operator.hpp>
12
#include <dune/amdis/common/Mpl.hpp>
13
#include <dune/amdis/common/TypeDefs.hpp>
14
15
16

namespace AMDiS
{
17
18
19
  const Flag EQUAL_BASES      = 0x01L;
  const Flag EQUAL_LOCALBASES = 0x02L;

20
  template <class GlobalBasis>
21
22
  class Assembler
  {
23
    /// The grid view the global FE basis lives on
24
    using GridView = typename GlobalBasis::GridView;
25
26

  public:
27
    /// Constructor, stores a shared-pointer to the feSpaces
28
29
    Assembler(GlobalBasis& globalBasis,
              MatrixOperators<GlobalBasis>& matrixOperators,
30
31
              VectorOperators<GlobalBasis>& rhsOperators,
              Constraints<GlobalBasis>& constraints)
32
      : globalBasis_(globalBasis)
33
34
      , matrixOperators_(matrixOperators)
      , rhsOperators_(rhsOperators)
35
      , constraints_(constraints)
36
    {}
37

38
39
40
41
42
    void update(GridView const& gv)
    {
      globalBasis_.update(gv);
    }

43
    /// Assemble the linear system
44
    template <class SystemMatrixType, class SystemVectorType>
45
46
47
48
    void assemble(
        SystemMatrixType& matrix,
        SystemVectorType& solution,
        SystemVectorType& rhs,
49
        bool asmMatrix, bool asmVector);
50
51

  private:
52
    /// Sets the system to zero and initializes all operators and boundary conditions
53
    template <class SystemMatrixType, class SystemVectorType>
54
55
56
57
    void initMatrixVector(
        SystemMatrixType& matrix,
        SystemVectorType& solution,
        SystemVectorType& rhs,
58
        bool asmMatrix, bool asmVector) const;
59
60


61
62
63
64
    template <class ElementContainer, class Container, class Operators, class... Bases>
    void assembleElementOperators(
        ElementContainer& elementContainer,
        Container& container,
65
        Operators& operators,
66
        Bases const&... subBases) const;
67
68
69
70


    /// Finish insertion into the matrix and assembles boundary conditions
    /// Return the number of nonzeros assembled into the matrix
71
    template <class SystemMatrixType, class SystemVectorType>
72
73
74
75
    std::size_t finishMatrixVector(
        SystemMatrixType& matrix,
        SystemVectorType& solution,
        SystemVectorType& rhs,
76
77
        bool asmMatrix, bool asmVector) const;

78

79
80
81
82
83
84
    /// Return whether the matrix-block needs to be assembled
    template <class LocalView0, class... LovalViews>
    auto const& getElement(LocalView0 const& localView, LovalViews const&...) const
    {
      return localView.element();
    }
85
86

    /// Return whether the matrix-block needs to be assembled
87
88
    template <class LocalView0, class... LovalViews>
    auto const& getGridView(LocalView0 const& localView, LovalViews const&...) const
89
    {
90
      return globalBasis_.gridView();
91
    }
92

93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
  public:
    template <class RowNode, class ColNode>
    Flag optimizationFlags(RowNode const& rowNode, ColNode const& colNode) const
    {
      Flag flag;
      if (rowNode.treeIndex() == colNode.treeIndex())
        flag |= EQUAL_BASES;

      // NOTE: find a way to compare localBases directly
//       if (rowNode.finiteElement().localBasis().order() == colNode.finiteElement().localBasis().order())
//         flag |= EQUAL_LOCALBASES;

      return flag;
    }

108
  private:
109
110
111
    GlobalBasis& globalBasis_;
    MatrixOperators<GlobalBasis>& matrixOperators_;
    VectorOperators<GlobalBasis>& rhsOperators_;
112
    Constraints<GlobalBasis>& constraints_;
113
114

    //TODO: add caching of localBases
115
116
117
  };

} // end namespace AMDiS
118
119

#include "Assembler.inc.hpp"