DirichletBC.inc.hpp 2.13 KB
Newer Older
1
2
3
#pragma once

#include <dune/functions/functionspacebases/interpolate.hh>
4
#include <dune/functions/functionspacebases/subspacebasis.hh>
5
#include <dune/amdis/LinearAlgebra.hpp>
6
#include <dune/amdis/linear_algebra/HierarchicWrapper.hpp>
7
8
9

namespace AMDiS
{
10
11
12
13
  template <class WorldVector, class Range>
    template <class RowBasis>
  void DirichletBC<WorldVector, Range>::initImpl(
      RowBasis const& rowBasis, Dune::TypeTree::LeafNodeTag)
14
15
  {
    using Dune::Functions::interpolate;
16

17
18
19
    if (!initialized_) {
      interpolate(rowBasis, hierarchicVectorWrapper(dirichletNodes_), predicate_);
      initialized_ = true;
20
    }
21
  }
22

23
24
25
26
27
28
  template <class WorldVector, class Range>
    template <class RowBasis>
  void DirichletBC<WorldVector, Range>::initImpl(
      RowBasis const& rowBasis, Dune::TypeTree::PowerNodeTag)
  {
    using Dune::Functions::interpolate;
29

30
    if (!initialized_) {
31
32
33
34
      auto tp = rowBasis.prefixPath();
      auto const& basis = rowBasis.rootBasis();
      for (std::size_t i = 0; i < degree(rowBasis.localView().tree()); ++i)
        interpolate(Dune::Functions::subspaceBasis(basis, push_back(tp,i)),
35
36
          hierarchicVectorWrapper(dirichletNodes_), predicate_);
      initialized_ = true;
37
38
39
40
41
    }
  }


  template <class WorldVector, class Range>
42
43
    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis, class ValueCat>
  void DirichletBC<WorldVector, Range>::finishImpl(
44
      Matrix& matrix, VectorX& solution, VectorB& rhs,
45
      RowBasis const& rowBasis, ColBasis const& colBasis, ValueCat)
46
47
  {
    using Dune::Functions::interpolate;
48

49
50
    test_exit_dbg(initialized_, "Boundary condition not initialized!");
    auto columns = matrix.applyDirichletBC(dirichletNodes_);
51

52
53
54
55
    interpolate(rowBasis, hierarchicVectorWrapper(rhs.getVector()), values_,
      hierarchicVectorWrapper(dirichletNodes_));
    interpolate(colBasis, hierarchicVectorWrapper(solution.getVector()), values_,
      hierarchicVectorWrapper(dirichletNodes_));
56

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

62
63

} // end namespace AMDiS