Assembler.inc.hpp 6.03 KB
Newer Older
1
2
3
4
#pragma once

namespace AMDiS {

5
6
7
template <class GlobalBasis>
  template <class MatrixOperators, class VectorOperators>
void Assembler<GlobalBasis>::assemble(
8
    GridView const& gv,
9
10
11
    SystemMatrixType& matrix,
    SystemVectorType& solution,
    SystemVectorType& rhs,
12
13
    MatrixOperators& matrixOperators,
    VectorOperators& rhsOperators)
14
{
15
  // 1. Update global bases
16
  globalBasis_.update(gv);
17

18
  // 2. init matrix and rhs vector and initialize dirichlet boundary conditions
19
  initMatrixVector(matrix, solution, rhs, matrixOperators, rhsOperators);
20
21

  // 3. traverse grid and assemble operators on the elements
22
  for (auto const& element : elements(globalBasis_.gridView()))
23
  {
24
25
26
    globalBasis_.bind(element);
    assembleElement(matrix, rhs, matrixOperators, rhsOperators);
    globalBasis_.unbind();
27
28
29
  }

  // 4. finish matrix insertion and apply dirichlet boundary conditions
30
  std::size_t nnz = finishMatrixVector(matrix, solution, rhs, matrixOperators, rhsOperators);
31
32
33
34
35

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


36
37
38
template <class GlobalBasis>
  template <class MatrixOperators, class VectorOperators>
void Assembler<GlobalBasis>::initMatrixVector(
39
40
41
    SystemMatrixType& matrix,
    SystemVectorType& solution,
    SystemVectorType& rhs,
42
43
    MatrixOperators& matrixOperators,
    VectorOperators& rhsOperators) const
44
{
45
46
47
  matrix.init(globalBasis_);
  solution.init(globalBasis_);
  rhs.init(globalBasis_);
48

49
50
51
52
53
  auto localView = globalBasis_.localView();
  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
  {
    if (rowNode.isLeaf)
      msg(globalBasis_.size(rowTreePath), " DOFs for Basis[", Dune::TypeTree::treePathIndex(rowTreePath), "]");
54

55
56
57
    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
    if (rhsOperators[rowNode].assemble(asmVector_))
      rhsOperators[rowNode].init(rowBasis);
58

59
    forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
60
    {
61
62
63
64
65
66
67
68
      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
      if (matrixOperators[rowNode][colNode].assemble(asmMatrix_))
        matrixOperators[rowNode][colNode].init(rowBasis, colBasis);

      // 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]);
69
70
71
72
73
    });
  });
}


74
75
76
77
78
79
80
template <class GlobalBasis>
  template <class MatrixOperators, class VectorOperators>
void Assembler<GlobalBasis>::assembleElement(
    SystemMatrixType& matrix,
    SystemVectorType& rhs,
    MatrixOperators& matrixOperators,
    VectorOperators& rhsOperators) const
81
{
82
83
  auto localView = globalBasis_.localView();
  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
84
  {
85
    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
86

87
88
89
90
    auto& rhsOp = rhsOperators[rowNode];
    if (rhsOp.assemble(asmVector_) && !rhsOp.empty()) {
      auto elementVector = makeElementVector(rowNode);
      set_to_zero(elementVector);
91

92
      assembleElementOperators(elementVector, rhs, rhsOp, rowBasis);
93
94
    }

95
96
97
98
99
    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);
100

101
102
        auto elementMatrix = makeElementMatrix(rowNode, colNode);
        set_to_zero(elementMatrix);
103

104
105
106
107
        assembleElementOperators(elementMatrix, matrix, matOp, rowBasis, colBasis);
      }
    });
  });
108
109
110
}


111
112
113
114
115
template <class GlobalBasis>
  template <class ElementContainer, class Container, class Operators, class... Bases>
void Assembler<GlobalBasis>::assembleElementOperators(
    ElementContainer& elementContainer,
    Container& container,
116
    Operators& operators,
117
    Bases const&... subBases) const
118
{
119
  if (operators.empty())
120
121
    return false;

122
123
  auto const& element = getElement(subBases...);
  auto const& gridView = subBasis.gridView();
124
125
126

  bool add = false;

127
  auto assemble_operators = [&](auto const& context, auto& operator_list) {
128
    for (auto scaled : operator_list) {
129
      bool add_op = scaled.op->assemble(gridView, context, subBases.localView()..., elementContainer, scaled.factor);
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
      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);
    }
  }

149
150
  if (!add)
    return;
151

152
  elementContainer.apply(subBases.localIndexSet()..., container);
153
154
155
}


156
157
158
template <class GlobalBasis>
  template <class MatrixOperators, class VectorOperators>
std::size_t Assembler<GlobalBasis>::finishMatrixVector(
159
160
161
    SystemMatrixType& matrix,
    SystemVectorType& solution,
    SystemVectorType& rhs,
162
163
    MatrixOperators& matrixOperators,
    VectorOperators& rhsOperators) const
164
{
165
166
167
168
  matrix.finish();

  auto localView = globalBasis_.localView();
  forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
169
  {
170
171
172
    auto& rhsOp = rhsOperators[rowNode];
    if (rhsOp.assemble(asmVector_))
      rhsOp.assembled = true;
173

174
    forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
175
    {
176
177
178
179
180
181
182
183
184
185
186
187
188
      auto& matOp = matrixOperators[rowNode][colNode];
      if (matOp.assemble(asmMatrix_))
        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]);
      //   }
      // }
189
190
191
    });
  });

192
  return matrix.getNnz();
193
194
195
}

} // end namespace AMDiS