Assembler.inc.hpp 6.89 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
  // 2. create a local matrix and vector
  std::size_t localSize = localView.maxSize();

26
27
  Impl::ElementMatrix elementMatrix(localSize, localSize);
  Impl::ElementVector 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
    localView.bind(element);
36
    auto geometry = element.geometry();
37
38
39
40
41
42
43
44
45
46

    // 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())
47
        assembleElementOperators(elementVector, rhs, rhsOp, geometry, rowLocalView);
48
49
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();
          colLocalView.bind(element); // NOTE: Is this necessary

57
          assembleElementOperators(elementMatrix, matrix, matOp, geometry, rowLocalView, colLocalView);
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 Geometry, class... LocalViews>
96
97
98
void Assembler<GlobalBasis>::assembleElementOperators(
    ElementContainer& elementContainer,
    Container& container,
99
    Operators& operators,
100
    Geometry const& geometry,
101
    LocalViews const&... localViews) const
102
{
103
104
  auto const& element = getElement(localViews...);
  auto const& gridView = getGridView(localViews...);
105
106
107

  bool add = false;

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


134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
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);
153
154
    //if (rhsOperators_[rowNode].assemble(asmVector))
    //  rhsOperators_[rowNode].init(rowBasis);
155
156
157
158

    forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
    {
      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
159
160
      //if (matrixOperators_[rowNode][colNode].assemble(asmMatrix))
      //  matrixOperators_[rowNode][colNode].init(rowBasis, colBasis);
161
162
163
164
165
166
167
168
169
170
171

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


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

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

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

      // finish boundary condition
198
199
200
201
      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);
202
203
204
    });
  });

205
  return matrix.getNnz();
206
207
208
}

} // end namespace AMDiS