DirichletBC.inc.hpp 1.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#pragma once

#include <type_traits>

#include <dune/common/hybridutilities.hh>
#include <dune/functions/functionspacebases/boundarydofs.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>

#include <amdis/Output.hpp>
#include <amdis/linear_algebra/Constraints.hpp>

namespace AMDiS {

template <class D, class R>
  template <class RowBasis, class ColBasis>
void DirichletBC<D,R>::
init(RowBasis const& rowBasis, ColBasis const& colBasis)
{
  using Dune::Functions::forEachBoundaryDOF;
  using LV = typename ColBasis::LocalView;
  using I = typename ColBasis::GridView::Intersection;
  dirichletNodes_.resize(colBasis.dimension());
  forEachBoundaryDOF(colBasis, [&](int localIndex, LV const& localView, I const& intersection) {
    dirichletNodes_[localView.index(localIndex)] = onBoundary(intersection);
  });
}


template <class D, class R>
  template <class Matrix, class VectorX, class VectorB, class RowNode, class ColNode>
void DirichletBC<D,R>::
fillBoundaryCondition(Matrix& matrix, VectorX& solution, VectorB& rhs,
                      RowNode const& rowNode, ColNode const& colNode)
{
  auto columns = Constraints<Matrix>::dirichletBC(matrix, dirichletNodes_);

  Dune::Hybrid::ifElse(std::is_same<RangeType_t<RowNode>, R>{},
    [&](auto id) {
      auto subBasis = Dune::Functions::subspaceBasis(rhs.basis(), rowNode.treePath());
      Dune::Functions::interpolate(subBasis, rhs, values_, dirichletNodes_);
    });

  // subtract columns of dirichlet nodes from rhs
  for (auto const& triplet : columns)
    rhs[triplet.row] -= triplet.value * solution[triplet.col];

  Dune::Hybrid::ifElse(std::is_same<RangeType_t<ColNode>, R>{},
    [&](auto id) {
      auto subBasis = Dune::Functions::subspaceBasis(solution.basis(), colNode.treePath());
      Dune::Functions::interpolate(subBasis, solution, values_, dirichletNodes_);
    });
}

} // end namespace AMDiS