Assembler.inc.hpp 6.17 KB
Newer Older
1
2
#pragma once

3
4
#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>
5
#include <dune/functions/functionspacebases/subspacebasis.hh>
6
#include <amdis/utility/TreePath.hpp>
7
#include <amdis/utility/Visitor.hpp>
8

9
#include <amdis/common/Math.hpp>
10

11
12
namespace AMDiS {

13
template <class Traits>
14
  template <class SystemMatrixType, class SystemVectorType>
15
void Assembler<Traits>::assemble(
16
17
18
    SystemMatrixType& matrix,
    SystemVectorType& solution,
    SystemVectorType& rhs,
19
    bool asmMatrix, bool asmVector)
20
{
21
22
23
24
  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
  initMatrixVector(matrix, solution, rhs, asmMatrix, asmVector);

  auto localView = globalBasis_.localView();
25

26
27
28
  // 2. create a local matrix and vector
  std::size_t localSize = localView.maxSize();

29
30
  Dune::DynamicMatrix<double> elementMatrix(localSize, localSize);
  Dune::DynamicVector<double> elementVector(localSize);
31
32

  // 3. traverse grid and assemble operators on the elements
33
  for (auto const& element : elements(globalBasis_.gridView()))
34
  {
35
36
    elementMatrix = 0;
    elementVector = 0;
37

38
    localView.bind(element);
39
    auto geometry = element.geometry();
40
41

    // traverse type-tree of global-basis
42
    forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
43
44
45
    {
      auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
      auto rowLocalView = rowBasis.localView();
46
      rowLocalView.bind(element);
47
48

      auto& rhsOp = rhsOperators_[rowNode];
49
      if (rhsOp.doAssemble(asmVector) && !rhsOp.empty()) {
50
51
52
53
        rhsOp.bind(element, geometry);

        auto vecAssembler = [&](auto const& context, auto& operator_list) {
          for (auto scaled : operator_list)
54
            scaled.op->assemble(context, rowLocalView.tree(), elementVector);
55
56
57
58
        };

        this->assembleElementOperators(element, rhsOp, vecAssembler);
      }
59

60
      forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
61
62
      {
        auto& matOp = matrixOperators_[rowNode][colNode];
63
        if (matOp.doAssemble(asmMatrix) && !matOp.empty()) {
64
65
          auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
          auto colLocalView = colBasis.localView();
66
67
68
69
70
71
          colLocalView.bind(element);

          matOp.bind(element, geometry);

          auto matAssembler = [&](auto const& context, auto& operator_list) {
            for (auto scaled : operator_list)
72
              scaled.op->assemble(context, rowLocalView.tree(), colLocalView.tree(), elementMatrix);
73
          };
74

75
          this->assembleElementOperators(element, matOp, matAssembler);
76
77
78
79
        }
      });
    });

80
81
82
    // add element-matrix to system-matrix and element-vector to rhs
    matrix.insert(localView, localView, elementMatrix);
    rhs.insert(localView, elementVector);
83

84
    // unbind all operators
85
    forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto&&) {
86
      rhsOperators_[rowNode].unbind();
87
      forEachNode_(localView.tree(), [&,this](auto const& colNode, auto&&) {
88
89
90
91
        matrixOperators_[rowNode][colNode].unbind();
      });
    });

92
    localView.unbind();
93
94
95
  }

  // 4. finish matrix insertion and apply dirichlet boundary conditions
96
  std::size_t nnz = finishMatrixVector(matrix, solution, rhs, asmMatrix, asmVector);
97

98
  msg("fill-in of assembled matrix: ", nnz);
99
100
101
}


102
template <class Traits>
103
  template <class Element, class Operators, class ElementAssembler>
104
void Assembler<Traits>::assembleElementOperators(
105
    Element const& element,
106
    Operators& operators,
107
    ElementAssembler const& localAssembler) const
108
{
109
  // assemble element operators
110
  localAssembler(element, operators.element);
111
112
113
114
115

  // assemble intersection operators
  if (!operators.intersection.empty()
      || (!operators.boundary.empty() && element.hasBoundaryIntersections()))
  {
116
    for (auto const& intersection : intersections(globalBasis_.gridView(), element)) {
117
      if (intersection.boundary())
118
        localAssembler(intersection, operators.boundary);
119
      else
120
        localAssembler(intersection, operators.intersection);
121
122
123
124
125
    }
  }
}


126
template <class Traits>
127
  template <class SystemMatrixType, class SystemVectorType>
128
void Assembler<Traits>::initMatrixVector(
129
130
131
132
133
134
135
136
    SystemMatrixType& matrix,
    SystemVectorType& solution,
    SystemVectorType& rhs,
    bool asmMatrix, bool asmVector) const
{
  matrix.init(asmMatrix);
  solution.compress();
  rhs.compress();
137
138
  if (asmVector)
    rhs = 0;
139
140

  auto localView = globalBasis_.localView();
141
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
142
  {
143
#ifdef HAVE_EXTENDED_DUNE_FUNCTIONS
144
    if (rowNode.isLeaf)
145
      msg(globalBasis_.dimension(rowTreePath), " DOFs for Basis[", to_string(rowTreePath), "]");
146
#endif
147
148
149

    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);

150
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
151
152
153
    {
      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);

154
      for (auto bc : constraints_[rowNode][colNode])
155
156
157
158
159
160
161
        bc->init(matrix, solution, rhs, rowBasis, colBasis);
    });
  });
  msg(globalBasis_.dimension(), " total DOFs");
}


162
template <class Traits>
163
  template <class SystemMatrixType, class SystemVectorType>
164
std::size_t Assembler<Traits>::finishMatrixVector(
165
166
167
    SystemMatrixType& matrix,
    SystemVectorType& solution,
    SystemVectorType& rhs,
168
    bool asmMatrix, bool asmVector) const
169
{
170
171
172
  matrix.finish();

  auto localView = globalBasis_.localView();
173
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
174
  {
175
    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
176
    auto& rhsOp = rhsOperators_[rowNode];
177
    if (rhsOp.doAssemble(asmVector))
178
      rhsOp.assembled = true;
179

180
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
181
    {
182
      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
183
      auto& matOp = matrixOperators_[rowNode][colNode];
184
      if (matOp.doAssemble(asmMatrix))
185
186
187
        matOp.assembled = true;

      // finish boundary condition
188
      for (auto bc : constraints_[rowNode][colNode])
189
        bc->finish(matrix, solution, rhs, rowBasis, colBasis);
190
191
192
    });
  });

193
  return matrix.nnz();
194
195
196
}

} // end namespace AMDiS