Assembler.inc.hpp 7.81 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
41
42
43
44

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

      auto& rhsOp = rhsOperators_[rowNode];
      if (rhsOp.assemble(asmVector) && !rhsOp.empty())
49
        this->assembleElementOperators(elementVector, rhs, rhsOp, geometry, rowLocalView);
50
51
52
53
54
55
56

      forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
      {
        auto& matOp = matrixOperators_[rowNode][colNode];
        if (matOp.assemble(asmMatrix) && !matOp.empty()) {
          auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
          auto colLocalView = colBasis.localView();
57
          colLocalView.bind(element); // NOTE: Is this necessary?
58

59
          this->assembleElementOperators(elementMatrix, matrix, matOp, geometry, rowLocalView, colLocalView);
60
61
62
63
64
65
66
67
68
69
        }
      });
    });

    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) {
70
        if (std::abs(elementMatrix(i,j)) > threshold<double>) {
71
72
73
          auto const col = localIndexSet.index(j);
          matrix(row,col) += elementMatrix(i,j);
        }
74
75
76
77
78
      }
    }

    // add element-vector to system-vector
    for (std::size_t i = 0; i < localView.size(); ++i) {
79
      if (std::abs(elementVector[i]) > threshold<double>) {
80
81
82
        auto const idx = localIndexSet.index(i);
        rhs[idx] += elementVector[i];
      }
83
84
85
86
    }

    localIndexSet.unbind();
    localView.unbind();
87
88
89
  }

  // 4. finish matrix insertion and apply dirichlet boundary conditions
90
  std::size_t nnz = finishMatrixVector(matrix, solution, rhs, asmMatrix, asmVector);
91
92
93
94
95

  msg("fillin of assembled matrix: ", nnz);
}


96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
template <class Traits>
  template <class ElementContainer, class Container, class Operators, class Geometry, class LocalView>
void Assembler<Traits>::assembleElementOperators(
    ElementContainer& elementContainer,
    Container& container,
    Operators& operators,
    Geometry const& geometry,
    LocalView const& localView) const
{
  auto const& element = getElement(localView);
  auto const& gridView = getGridView(localView);

  bool add = false;

  auto assemble_operators = [&](auto const& context, auto& operator_list) {
    for (auto scaled : operator_list) {
      scaled.op->bind(element, geometry);
      bool add_op = scaled.op->assemble(context, elementContainer, localView.tree());
      scaled.op->unbind();
      add = add || add_op;
    }
  };

  // assemble element operators
  assemble_operators(element, operators.element);

  // assemble intersection operators
  if (!operators.intersection.empty()
      || (!operators.boundary.empty() && element.hasBoundaryIntersections()))
  {
    for (auto const& intersection : intersections(gridView, element)) {
      if (intersection.boundary())
        assemble_operators(intersection, operators.boundary);
      else
        assemble_operators(intersection, operators.intersection);
    }
  }
}


136
template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
137
  template <class ElementContainer, class Container, class Operators, class Geometry, class RowLocalView, class ColLocalView>
138
void Assembler<Traits>::assembleElementOperators(
139
140
    ElementContainer& elementContainer,
    Container& container,
141
    Operators& operators,
142
    Geometry const& geometry,
Praetorius, Simon's avatar
Praetorius, Simon committed
143
    RowLocalView const& rowLocalView, ColLocalView const& colLocalView) const
144
{
145
146
  auto const& element = getElement(rowLocalView, colLocalView);
  auto const& gridView = getGridView(rowLocalView, colLocalView);
147
148
149

  bool add = false;

150
  auto assemble_operators = [&](auto const& context, auto& operator_list) {
151
    for (auto scaled : operator_list) {
152
      scaled.op->bind(element, geometry);
153
      bool add_op = scaled.op->assemble(context, elementContainer, rowLocalView.tree(), colLocalView.tree());
154
      scaled.op->unbind();
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
      add = add || add_op;
    }
  };

  // assemble element operators
  assemble_operators(element, operators.element);

  // assemble intersection operators
  if (!operators.intersection.empty()
      || (!operators.boundary.empty() && element.hasBoundaryIntersections()))
  {
    for (auto const& intersection : intersections(gridView, element)) {
      if (intersection.boundary())
        assemble_operators(intersection, operators.boundary);
      else
        assemble_operators(intersection, operators.intersection);
    }
  }
}


176
template <class Traits>
177
  template <class SystemMatrixType, class SystemVectorType>
178
void Assembler<Traits>::initMatrixVector(
179
180
181
182
183
184
185
186
187
188
189
190
    SystemMatrixType& matrix,
    SystemVectorType& solution,
    SystemVectorType& rhs,
    bool asmMatrix, bool asmVector) const
{
  matrix.init(asmMatrix);
  solution.compress();
  rhs.compress();

  auto localView = globalBasis_.localView();
  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
  {
191
#ifdef HAVE_EXTENDED_DUNE_FUNCTIONS
192
    if (rowNode.isLeaf)
193
      msg(globalBasis_.dimension(rowTreePath), " DOFs for Basis[", to_string(rowTreePath), "]");
194
#endif
195
196
197
198
199
200
201

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

    forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
    {
      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);

202
      for (auto bc : constraints_[rowNode][colNode])
203
204
205
206
207
208
209
        bc->init(matrix, solution, rhs, rowBasis, colBasis);
    });
  });
  msg(globalBasis_.dimension(), " total DOFs");
}


210
template <class Traits>
211
  template <class SystemMatrixType, class SystemVectorType>
212
std::size_t Assembler<Traits>::finishMatrixVector(
213
214
215
    SystemMatrixType& matrix,
    SystemVectorType& solution,
    SystemVectorType& rhs,
216
    bool asmMatrix, bool asmVector) const
217
{
218
219
220
  matrix.finish();

  auto localView = globalBasis_.localView();
221
  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
222
  {
223
    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
224
225
    auto& rhsOp = rhsOperators_[rowNode];
    if (rhsOp.assemble(asmVector))
226
      rhsOp.assembled = true;
227

228
    forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
229
    {
230
      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
231
232
      auto& matOp = matrixOperators_[rowNode][colNode];
      if (matOp.assemble(asmMatrix))
233
234
235
        matOp.assembled = true;

      // finish boundary condition
236
      for (auto bc : constraints_[rowNode][colNode])
237
        bc->finish(matrix, solution, rhs, rowBasis, colBasis);
238
239
240
    });
  });

241
  return matrix.getNnz();
242
243
244
}

} // end namespace AMDiS