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

3
#include <dune/functions/functionspacebases/subspacebasis.hh>
4
#include <dune/amdis/utility/TreePath.hpp>
5
6
#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
  // 2. create a local matrix and vector
  std::size_t localSize = localView.maxSize();

  mtl::dense2D<double> elementMatrix(localSize, localSize);
  mtl::dense_vector<double> elementVector(localSize);
28
29

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

35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
    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

56
57
          assembleElementOperators(elementMatrix, matrix, matOp, rowLocalView, colLocalView,
                                   optimizationFlags(rowNode, colNode));
58
59
60
61
62
63
64
65
66
67
        }
      });
    });

    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) {
68
69
70
71
        if (elementMatrix(i,j) != 0.0) {
          auto const col = localIndexSet.index(j);
          matrix(row,col) += elementMatrix(i,j);
        }
72
73
74
75
76
      }
    }

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

    localIndexSet.unbind();
    localView.unbind();
85
86
87
  }

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

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


94
template <class GlobalBasis>
95
  template <class ElementContainer, class Container, class Operators, class... LocalViews>
96
97
98
void Assembler<GlobalBasis>::assembleElementOperators(
    ElementContainer& elementContainer,
    Container& container,
99
    Operators& operators,
100
    LocalViews const&... localViews) const
101
{
102
103
  auto const& element = getElement(localViews...);
  auto const& gridView = getGridView(localViews...);
104
105
106

  bool add = false;

107
  auto assemble_operators = [&](auto const& context, auto& operator_list) {
108
    for (auto scaled : operator_list) {
109
      bool add_op = scaled.op->assemble(context, elementContainer, localViews..., scaled.factor);
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
      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);
    }
  }
}


131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
template <class GlobalBasis>
  template <class SystemMatrixType, class SystemVectorType>
void Assembler<GlobalBasis>::initMatrixVector(
    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)
  {
    if (rowNode.isLeaf)
      msg(0, " DOFs for Basis[", to_string(rowTreePath), "]"); // TODO: add right values

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

    forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
    {
      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
      if (matrixOperators_[rowNode][colNode].assemble(asmMatrix))
        matrixOperators_[rowNode][colNode].init(rowBasis, colBasis);

      for (auto bc : constraints_[rowNode][colNode].scalar)
        bc->init(matrix, solution, rhs, rowBasis, colBasis);
      for (auto bc : constraints_[rowNode][colNode].vector)
        bc->init(matrix, solution, rhs, rowBasis, colBasis);
    });
  });
  msg(globalBasis_.dimension(), " total DOFs");
}


169
template <class GlobalBasis>
170
  template <class SystemMatrixType, class SystemVectorType>
171
std::size_t Assembler<GlobalBasis>::finishMatrixVector(
172
173
174
    SystemMatrixType& matrix,
    SystemVectorType& solution,
    SystemVectorType& rhs,
175
    bool asmMatrix, bool asmVector) const
176
{
177
178
179
  matrix.finish();

  auto localView = globalBasis_.localView();
180
  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
181
  {
182
    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
183
184
    auto& rhsOp = rhsOperators_[rowNode];
    if (rhsOp.assemble(asmVector))
185
      rhsOp.assembled = true;
186

187
    forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
188
    {
189
      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
190
191
      auto& matOp = matrixOperators_[rowNode][colNode];
      if (matOp.assemble(asmMatrix))
192
193
194
        matOp.assembled = true;

      // finish boundary condition
195
196
197
198
      for (auto bc : constraints_[rowNode][colNode].scalar)
        bc->finish(matrix, solution, rhs, rowBasis, colBasis);
      for (auto bc : constraints_[rowNode][colNode].vector)
        bc->finish(matrix, solution, rhs, rowBasis, colBasis);
199
200
201
    });
  });

202
  return matrix.getNnz();
203
204
205
}

} // end namespace AMDiS