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