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

3
#include <dune/functions/functionspacebases/subspacebasis.hh>
4
5
#include <amdis/utility/TreePath.hpp>
#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
23
  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
  initMatrixVector(matrix, solution, rhs, asmMatrix, asmVector);

  auto localView = globalBasis_.localView();
  auto localIndexSet = globalBasis_.localIndexSet();
24

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

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

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

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

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

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

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

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

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

          matOp.bind(element, geometry);

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

74
          this->assembleElementOperators(element, matOp, matAssembler);
75
76
77
78
79
80
81
82
83
84
        }
      });
    });

    localIndexSet.bind(localView);

    // add element-matrix to system-matrix
    for (std::size_t i = 0; i < localView.size(); ++i) {
      auto const row = localIndexSet.index(i);
      for (std::size_t j = 0; j < localView.size(); ++j) {
85
        if (std::abs(elementMatrix(i,j)) > threshold<double>) {
86
87
88
          auto const col = localIndexSet.index(j);
          matrix(row,col) += elementMatrix(i,j);
        }
89
90
91
92
93
      }
    }

    // add element-vector to system-vector
    for (std::size_t i = 0; i < localView.size(); ++i) {
94
      if (std::abs(elementVector[i]) > threshold<double>) {
95
96
97
        auto const idx = localIndexSet.index(i);
        rhs[idx] += elementVector[i];
      }
98
99
    }

100
    // unbind all operators
101
    AMDiS::forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto&&) {
102
      rhsOperators_[rowNode].unbind();
103
      AMDiS::forEachNode_(localView.tree(), [&,this](auto const& colNode, auto&&) {
104
105
106
107
        matrixOperators_[rowNode][colNode].unbind();
      });
    });

108
109
    localIndexSet.unbind();
    localView.unbind();
110
111
112
  }

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

115
  msg("fill-in of assembled matrix: ", nnz);
116
117
118
}


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

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


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

  auto localView = globalBasis_.localView();
158
  AMDiS::forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
159
  {
160
#ifdef HAVE_EXTENDED_DUNE_FUNCTIONS
161
    if (rowNode.isLeaf)
162
      msg(globalBasis_.dimension(rowTreePath), " DOFs for Basis[", to_string(rowTreePath), "]");
163
#endif
164
165
166

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

167
    AMDiS::forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
168
169
170
    {
      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);

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


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

  auto localView = globalBasis_.localView();
190
  AMDiS::forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
191
  {
192
    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
193
    auto& rhsOp = rhsOperators_[rowNode];
194
    if (rhsOp.doAssemble(asmVector))
195
      rhsOp.assembled = true;
196

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

      // finish boundary condition
205
      for (auto bc : constraints_[rowNode][colNode])
206
        bc->finish(matrix, solution, rhs, rowBasis, colBasis);
207
208
209
    });
  });

210
  return matrix.getNnz();
211
212
213
}

} // end namespace AMDiS