Assembler.inc.hpp 6.93 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
#include <dune/amdis/Math.hpp>

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
45
46
47
48

    // 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())
49
        assembleElementOperators(elementVector, rhs, rhsOp, geometry, rowLocalView);
50
51
52
53
54
55
56
57
58

      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

59
          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
template <class Traits>
97
  template <class ElementContainer, class Container, class Operators, class Geometry, class... LocalViews>
98
void Assembler<Traits>::assembleElementOperators(
99
100
    ElementContainer& elementContainer,
    Container& container,
101
    Operators& operators,
102
    Geometry const& geometry,
103
    LocalViews const&... localViews) const
104
{
105
106
  auto const& element = getElement(localViews...);
  auto const& gridView = getGridView(localViews...);
107
108
109

  bool add = false;

110
  auto assemble_operators = [&](auto const& context, auto& operator_list) {
111
    for (auto scaled : operator_list) {
112
113
114
      scaled.op->bind(element, geometry);
      bool add_op = scaled.op->assemble(context, elementContainer, localViews.tree()..., scaled.factor);
      scaled.op->unbind();
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
      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>
137
  template <class SystemMatrixType, class SystemVectorType>
138
void Assembler<Traits>::initMatrixVector(
139
140
141
142
143
144
145
146
147
148
149
150
151
    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)
152
      msg(globalBasis_.dimension(rowTreePath), " DOFs for Basis[", to_string(rowTreePath), "]");
153
154

    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
155
156
    //if (rhsOperators_[rowNode].assemble(asmVector))
    //  rhsOperators_[rowNode].init(rowBasis);
157
158
159
160

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

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


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

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

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

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

207
  return matrix.getNnz();
208
209
210
}

} // end namespace AMDiS