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

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

7
8
namespace AMDiS {

9
template <class GlobalBasis>
10
  template <class SystemMatrixType, class SystemVectorType>
11
void Assembler<GlobalBasis>::assemble(
12
13
14
    SystemMatrixType& matrix,
    SystemVectorType& solution,
    SystemVectorType& rhs,
15
    bool asmMatrix, bool asmVector)
16
{
17
18
19
20
21
  // 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();
22

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

  mtl::dense2D<double> elementMatrix(localSize, localSize);
  set_to_zero(elementMatrix);

  mtl::dense_vector<double> elementVector(localSize);
  set_to_zero(elementVector);
31
32

  // 3. traverse grid and assemble operators on the elements
33
  for (auto const& element : elements(globalBasis_.gridView()))
34
  {
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
    localView.bind(element);

    // 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();
      rowLocalView.bind(element); // NOTE: Is this necessary

      auto& rhsOp = rhsOperators_[rowNode];
      if (rhsOp.assemble(asmVector) && !rhsOp.empty())
        assembleElementOperators(elementVector, rhs, rhsOp, rowLocalView);

      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();
          colLocalView.bind(element); // NOTE: Is this necessary

          assembleElementOperators(elementMatrix, matrix, matOp, rowLocalView, colLocalView);
        }
      });
    });

    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) {
        auto const col = localIndexSet.index(j);
        matrix(row,col) += elementMatrix(i,j);
      }
    }

    // add element-vector to system-vector
    for (std::size_t i = 0; i < localView.size(); ++i) {
      auto const idx = localIndexSet.index(i);
      rhs[idx] += elementVector[i];
    }

    localIndexSet.unbind();
    localView.unbind();
80
81
82
  }

  // 4. finish matrix insertion and apply dirichlet boundary conditions
83
  std::size_t nnz = finishMatrixVector(matrix, solution, rhs, asmMatrix, asmVector);
84
85
86
87
88

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


89
template <class GlobalBasis>
90
  template <class SystemMatrixType, class SystemVectorType>
91
void Assembler<GlobalBasis>::initMatrixVector(
92
93
94
    SystemMatrixType& matrix,
    SystemVectorType& solution,
    SystemVectorType& rhs,
95
    bool asmMatrix, bool asmVector) const
96
{
97
98
99
  matrix.init(asmMatrix);
  solution.compress();
  rhs.compress();
100

101
102
103
104
  auto localView = globalBasis_.localView();
  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
  {
    if (rowNode.isLeaf)
105
      msg(0, " DOFs for Basis[", 1, "]"); // TODO: add right values
106

107
    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
108
109
    if (rhsOperators_[rowNode].assemble(asmVector))
      rhsOperators_[rowNode].init(rowBasis);
110

111
    forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
112
    {
113
      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
114
115
      if (matrixOperators_[rowNode][colNode].assemble(asmMatrix))
        matrixOperators_[rowNode][colNode].init(rowBasis, colBasis);
116
117
118
119
120

      // init boundary condition
      // for (int c = 0; c < nComponents; ++c)
      //   for (auto bc : matrixOperators[R][c].dirichlet)
      //     bc->init(c == C, matrix(_r, _c), solution[_c], rhs[_r]);
121
122
123
124
125
    });
  });
}


126
template <class GlobalBasis>
127
  template <class ElementContainer, class Container, class Operators, class... LocalViews>
128
129
130
void Assembler<GlobalBasis>::assembleElementOperators(
    ElementContainer& elementContainer,
    Container& container,
131
    Operators& operators,
132
    LocalViews const&... localViews) const
133
{
134
135
  auto const& element = getElement(localViews...);
  auto const& gridView = getGridView(localViews...);
136
137
138

  bool add = false;

139
  auto assemble_operators = [&](auto const& context, auto& operator_list) {
140
    for (auto scaled : operator_list) {
141
      bool add_op = scaled.op->assemble(context, localViews..., elementContainer, scaled.factor);
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
      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);
    }
  }
}


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

  auto localView = globalBasis_.localView();
174
  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto /*rowTreePath*/)
175
  {
176
177
    auto& rhsOp = rhsOperators_[rowNode];
    if (rhsOp.assemble(asmVector))
178
      rhsOp.assembled = true;
179

180
    forEachNode(localView.tree(), [&,this](auto const& colNode, auto /*colTreePath*/)
181
    {
182
183
      auto& matOp = matrixOperators_[rowNode][colNode];
      if (matOp.assemble(asmMatrix))
184
185
186
187
188
189
190
191
192
193
194
        matOp.assembled = true;

      // finish boundary condition
      // for (int c = 0; c < nComponents; ++c) {
      //   for (int r = 0; r < nComponents; ++r) {
      //     if (r != R && c != C)
      //       continue;
      //     for (auto bc : matrixOperators[r][c].dirichlet)
      //       bc->finish(r == R, c == C, matrix(_r, _c), solution[_c], rhs[_r]);
      //   }
      // }
195
196
197
    });
  });

198
  return matrix.getNnz();
199
200
201
}

} // end namespace AMDiS