Constraints.hpp 2.11 KB
Newer Older
1
2
3
4
5
6
#pragma once

#include <list>

#include <Eigen/SparseCore>

7
#include <amdis/common/Index.hpp>
8
#include <amdis/linearalgebra/Constraints.hpp>
9
10
#include <amdis/linearalgebra/eigen/MatrixBackend.hpp>
#include <amdis/linearalgebra/eigen/VectorBackend.hpp>
11
12
13

namespace AMDiS
{
Praetorius, Simon's avatar
Praetorius, Simon committed
14
15
  template <class T, int O>
  struct Constraints<EigenSparseMatrix<T,O>>
16
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
17
18
    using Matrix = EigenSparseMatrix<T,O>;
    using Vector = EigenVector<T>;
19
20

    template <class BitVector>
21
    static void dirichletBC(Matrix& mat, Vector& sol, Vector& rhs, BitVector const& nodes, bool setDiagonal = true)
22
    {
23
24
25
26
27
28
29
      clearDirichletRow(mat.matrix(), nodes, setDiagonal);

      // copy solution dirichlet data to rhs vector
      for (typename Vector::size_type i = 0; i < sol.vector().size(); ++i) {
        if (nodes[i])
          rhs.vector()[i] = sol.vector()[i];
      }
30
31
32
    }

    template <class BitVector, class Associations>
33
    static void periodicBC(Matrix& mat, Vector& sol, Vector& rhs, BitVector const& left, Associations const& left2right,
34
35
36
37
38
39
      bool setDiagonal = true)
    {
      error_exit("Not implemented");
    }

  protected:
Praetorius, Simon's avatar
Praetorius, Simon committed
40
    template <class BitVector>
41
    static void clearDirichletRow(Eigen::SparseMatrix<T, Eigen::ColMajor>& mat, BitVector const& nodes, bool setDiagonal)
42
    {
43
44
45
      using Mat = Eigen::SparseMatrix<T, Eigen::ColMajor>;
      for (typename Mat::Index c = 0; c < mat.outerSize(); ++c) {
        for (typename Mat::InnerIterator it(mat, c); it; ++it) {
46
47
48
49
50
51
52
53
          if (nodes[it.row()]) {
            it.valueRef() = (setDiagonal && it.row() == it.col() ? T(1) : T(0));
            break;
          }
        }
      }
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
54
    template <class BitVector>
55
    static void clearDirichletRow(Eigen::SparseMatrix<T, Eigen::RowMajor>& mat, BitVector const& nodes, bool setDiagonal)
56
    {
57
58
      using Mat = Eigen::SparseMatrix<T, Eigen::RowMajor>;
      for (typename Mat::Index r = 0; r < mat.outerSize(); ++r) {
59
        if (nodes[r]) {
60
          for (typename Mat::InnerIterator it(mat, r); it; ++it) {
61
62
63
64
65
66
67
68
            it.valueRef() = (setDiagonal && it.row() == it.col() ? T(1) : T(0));
          }
        }
      }
    }
  };

} // end namespace AMDiS