DOFMatrixBase.inc.hpp 3.35 KB
Newer Older
1
2
3
4
#pragma once

#include <amdis/Assembler.hpp>
#include <amdis/LocalOperator.hpp>
5
6
#include <amdis/typetree/Visitor.hpp>
#include <amdis/utility/AssembleOperators.hpp>
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

namespace AMDiS {

template <class RB, class CB, class B>
void DOFMatrixBase<RB,CB,B>::
init(bool asmMatrix)
{
  backend_.init(*rowBasis_, *colBasis_, asmMatrix);

  auto const rowSize = rowBasis_->localView().maxSize();
  auto const colSize = colBasis_->localView().maxSize();
  elementMatrix_.resize(rowSize, colSize);
}


template <class RB, class CB, class B>
void DOFMatrixBase<RB,CB,B>::
24
finish(bool /*asmMatrix*/)
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
{
  backend_.finish();
}


template <class RB, class CB, class B>
void DOFMatrixBase<RB,CB,B>::
insert(RowLocalView const& rowLocalView, ColLocalView const& colLocalView,
       ElementMatrix const& elementMatrix)
{
  for (size_type i = 0; i < rowLocalView.size(); ++i) {
    size_type const row = flatMultiIndex(rowLocalView.index(i));
    for (size_type j = 0; j < colLocalView.size(); ++j) {
      if (std::abs(elementMatrix[i][j]) > threshold<double>) {
        size_type const col = flatMultiIndex(colLocalView.index(j));
        backend_.insert(row, col, elementMatrix[i][j]);
      }
    }
  }
}


template <class RB, class CB, class B>
48
  template <class ContextTag, class PreOperator, class RowTreePath, class ColTreePath>
49
void DOFMatrixBase<RB,CB,B>::
50
addOperator(ContextTag contextTag, PreOperator const& preOp,
51
52
53
54
55
56
57
58
59
60
61
62
            RowTreePath row, ColTreePath col)
{
  static_assert( Concepts::PreTreePath<RowTreePath>,
      "row must be a valid treepath, or an integer/index-constant");
  static_assert( Concepts::PreTreePath<ColTreePath>,
      "col must be a valid treepath, or an integer/index-constant");

  auto i = child(rowBasis_->localView().tree(), makeTreePath(row));
  auto j = child(colBasis_->localView().tree(), makeTreePath(col));

  using Context = typename ContextTag::type;
  auto op = makeLocalOperator<Context>(preOp, rowBasis_->gridView());
63
  auto localAssembler = makeUniquePtr(makeAssembler<Context>(std::move(op), i, j));
64

65
  operators_[i][j].push(contextTag, std::move(localAssembler));
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
}


template <class RB, class CB, class B>
void DOFMatrixBase<RB,CB,B>::
assemble(RowLocalView const& rowLocalView, ColLocalView const& colLocalView)
{
  elementMatrix_ = 0;
  auto const& gv = rowBasis_->gridView();
  auto const& element = rowLocalView.element();
  auto geometry = element.geometry();

  forEachNode_(rowLocalView.tree(), [&](auto const& rowNode, auto) {
    forEachNode_(colLocalView.tree(), [&](auto const& colNode, auto) {
      auto& matOp = operators_[rowNode][colNode];
Praetorius, Simon's avatar
Praetorius, Simon committed
81
      if (matOp) {
82
        matOp.bind(element, geometry);
83
        assembleOperators(gv, element, matOp, makeMatrixAssembler(rowNode, colNode, elementMatrix_));
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
        matOp.unbind();
      }
    });
  });

  insert(rowLocalView, colLocalView, elementMatrix_);
}


template <class RB, class CB, class B>
void DOFMatrixBase<RB,CB,B>::
assemble()
{
  auto rowLocalView = rowBasis_->localView();
  auto colLocalView = colBasis_->localView();

  init(true);
  for (auto const& element : elements(rowBasis_->gridView())) {
    rowLocalView.bind(element);
    if (rowBasis_ == colBasis_)
      assemble(rowLocalView, rowLocalView);
    else {
      colLocalView.bind(element);
      assemble(rowLocalView, colLocalView);
      colLocalView.unbind();
    }
    rowLocalView.unbind(element);
  }
  finish(true);
}


} // end namespace AMDiS