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

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

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

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

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


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

  bool add = false;

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


130
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
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");
}


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

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

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

      // finish boundary condition
194
195
196
197
      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);
198
199
200
    });
  });

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

} // end namespace AMDiS