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

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

7
#include <amdis/common/Math.hpp>
8

9
10
namespace AMDiS {

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

  auto localView = globalBasis_.localView();
23

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

27
28
  Impl::ElementMatrix elementMatrix(localSize, localSize);
  Impl::ElementVector elementVector(localSize);
29
30

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

36
    localView.bind(element);
37
    auto geometry = element.geometry();
38
39

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

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

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

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

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

          matOp.bind(element, geometry);

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

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

    // add element-matrix to system-matrix
    for (std::size_t i = 0; i < localView.size(); ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
80
      auto const row = localView.index(i);
81
      for (std::size_t j = 0; j < localView.size(); ++j) {
82
        if (std::abs(elementMatrix(i,j)) > threshold<double>) {
Praetorius, Simon's avatar
Praetorius, Simon committed
83
          auto const col = localView.index(j);
84
85
          matrix(row,col) += elementMatrix(i,j);
        }
86
87
88
89
90
      }
    }

    // add element-vector to system-vector
    for (std::size_t i = 0; i < localView.size(); ++i) {
91
      if (std::abs(elementVector[i]) > threshold<double>) {
Praetorius, Simon's avatar
Praetorius, Simon committed
92
        auto const idx = localView.index(i);
93
94
        rhs[idx] += elementVector[i];
      }
95
96
    }

97
    // unbind all operators
98
    forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto&&) {
99
      rhsOperators_[rowNode].unbind();
100
      forEachNode_(localView.tree(), [&,this](auto const& colNode, auto&&) {
101
102
103
104
        matrixOperators_[rowNode][colNode].unbind();
      });
    });

105
    localView.unbind();
106
107
108
  }

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

111
  msg("fill-in of assembled matrix: ", nnz);
112
113
114
}


115
template <class Traits>
116
  template <class Element, class Operators, class ElementAssembler>
117
void Assembler<Traits>::assembleElementOperators(
118
    Element const& element,
119
    Operators& operators,
120
    ElementAssembler const& localAssembler) const
121
{
122
  // assemble element operators
123
  localAssembler(element, operators.element);
124
125
126
127
128

  // assemble intersection operators
  if (!operators.intersection.empty()
      || (!operators.boundary.empty() && element.hasBoundaryIntersections()))
  {
129
    for (auto const& intersection : intersections(globalBasis_.gridView(), element)) {
130
      if (intersection.boundary())
131
        localAssembler(intersection, operators.boundary);
132
      else
133
        localAssembler(intersection, operators.intersection);
134
135
136
137
138
    }
  }
}


139
template <class Traits>
140
  template <class SystemMatrixType, class SystemVectorType>
141
void Assembler<Traits>::initMatrixVector(
142
143
144
145
146
147
148
149
    SystemMatrixType& matrix,
    SystemVectorType& solution,
    SystemVectorType& rhs,
    bool asmMatrix, bool asmVector) const
{
  matrix.init(asmMatrix);
  solution.compress();
  rhs.compress();
150
151
  if (asmVector)
    rhs = 0;
152
153

  auto localView = globalBasis_.localView();
154
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
155
  {
156
#ifdef HAVE_EXTENDED_DUNE_FUNCTIONS
157
    if (rowNode.isLeaf)
158
      msg(globalBasis_.dimension(rowTreePath), " DOFs for Basis[", to_string(rowTreePath), "]");
159
#endif
160
161
162

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

163
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
164
165
166
    {
      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);

167
      for (auto bc : constraints_[rowNode][colNode])
168
169
170
171
172
173
174
        bc->init(matrix, solution, rhs, rowBasis, colBasis);
    });
  });
  msg(globalBasis_.dimension(), " total DOFs");
}


175
template <class Traits>
176
  template <class SystemMatrixType, class SystemVectorType>
177
std::size_t Assembler<Traits>::finishMatrixVector(
178
179
180
    SystemMatrixType& matrix,
    SystemVectorType& solution,
    SystemVectorType& rhs,
181
    bool asmMatrix, bool asmVector) const
182
{
183
184
185
  matrix.finish();

  auto localView = globalBasis_.localView();
186
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
187
  {
188
    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
189
    auto& rhsOp = rhsOperators_[rowNode];
190
    if (rhsOp.doAssemble(asmVector))
191
      rhsOp.assembled = true;
192

193
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
194
    {
195
      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
196
      auto& matOp = matrixOperators_[rowNode][colNode];
197
      if (matOp.doAssemble(asmMatrix))
198
199
200
        matOp.assembled = true;

      // finish boundary condition
201
      for (auto bc : constraints_[rowNode][colNode])
202
        bc->finish(matrix, solution, rhs, rowBasis, colBasis);
203
204
205
    });
  });

206
  return matrix.nnz();
207
208
209
}

} // end namespace AMDiS