Constraints.hpp 6.05 KB
Newer Older
1
2
#pragma once

3
4
#include <algorithm>
#include <vector>
5
6
7
8
9
10

#include <boost/numeric/mtl/matrix/compressed2D.hpp>
#include <boost/numeric/mtl/matrix/inserter.hpp>
#include <boost/numeric/mtl/utility/property_map.hpp>
#include <boost/numeric/mtl/utility/range_wrapper.hpp>

11
#include <amdis/linearalgebra/Constraints.hpp>
12
13
14
#include <amdis/linearalgebra/SymmetryStructure.hpp>
#include <amdis/linearalgebra/mtl/MatrixBackend.hpp>
#include <amdis/linearalgebra/mtl/VectorBackend.hpp>
15
16
17

namespace AMDiS
{
Praetorius, Simon's avatar
Praetorius, Simon committed
18
19
  template <class T>
  struct Constraints<MTLSparseMatrix<T>>
20
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
21
22
    using Matrix = MTLSparseMatrix<T>;
    using Vector = MTLVector<T>;
23
24

    template <class BitVector>
25
    static void dirichletBC(Matrix& mat, Vector& sol, Vector& rhs, BitVector const& nodes, bool setDiagonal = true)
26
    {
27
28
29
30
31
32
      SymmetryStructure const symmetry = mat.symmetry();
      if (symmetry == SymmetryStructure::spd || symmetry == SymmetryStructure::symmetric || symmetry == SymmetryStructure::hermitian)
        symmetricDirichletBC(mat.matrix(), sol.vector(), rhs.vector(), nodes, setDiagonal);
      else
        unsymmetricDirichletBC(mat.matrix(), sol.vector(), rhs.vector(), nodes, setDiagonal);
    }
33

34
35
36
    template <class Mat, class Vec, class BitVector>
    static void symmetricDirichletBC(Mat& mat, Vec& sol, Vec& rhs, BitVector const& nodes, bool setDiagonal = true)
    {
37
38
39
40
41
42
43
44
45
46
47
48
49
      // Define the property maps
      auto row   = mtl::mat::row_map(mat);
      auto col   = mtl::mat::col_map(mat);
      auto value = mtl::mat::value_map(mat);

      // iterate over the matrix
      std::size_t slotSize = 0;
      for (auto r : mtl::rows_of(mat)) {   // rows or columns
        std::size_t rowSize = 0;
        for (auto i : mtl::nz_of(r)) {          // non-zeros within
          ++rowSize;
          if (nodes[row(i)]) {
            // set identity row
50
            value(i, 0);
51
52
          }
          else if (setDiagonal && nodes[col(i)]) {
53
54
            rhs[row(i)] -= value(i) * sol[col(i)];
            value(i, 0);
55
56
57
58
59
60
61
          }
        }
        slotSize = std::max(slotSize, rowSize);
      }

      // set diagonal entry
      if (setDiagonal) {
62
        mtl::mat::inserter<Mat, mtl::update_store<typename Mat::value_type> > ins(mat, slotSize);
63
64
        for (std::size_t i = 0; i < nodes.size(); ++i) {
          if (nodes[i])
65
            ins[i][i] = 1;
66
67
68
        }
      }

69
70
71
72
73
      // copy solution dirichlet data to rhs vector
      for (typename Vec::size_type i = 0; i < mtl::size(sol); ++i) {
        if (nodes[i])
          rhs[i] = sol[i];
      }
74
    }
75
76
77

    template <class Mat, class Vec, class BitVector>
    static void unsymmetricDirichletBC(Mat& mat, Vec& sol, Vec& rhs, BitVector const& nodes, bool setDiagonal = true)
78
79
80
81
82
83
84
85
86
87
88
    {
      // Define the property maps
      auto row   = mtl::mat::row_map(mat);
      auto col   = mtl::mat::col_map(mat);
      auto value = mtl::mat::value_map(mat);

      // iterate over the matrix
      for (auto r : mtl::rows_of(mat)) {   // rows of the matrix
        if (nodes[r.value()]) {
          for (auto i : mtl::nz_of(r)) {          // non-zeros within
            // set identity row
89
            value(i, (setDiagonal && row(i) == col(i) ? 1 : 0) );
90
91
92
93
          }
        }
      }

94
95
96
97
98
      // copy solution dirichlet data to rhs vector
      for (typename Vec::size_type i = 0; i < mtl::size(sol); ++i) {
        if (nodes[i])
          rhs[i] = sol[i];
      }
99
    }
100

101
102
103
104
105
106
107
108
109
110
111


    template <class Associations>
    static std::size_t at(Associations const& m, std::size_t idx)
    {
      auto it = m.find(idx);
      assert(it != m.end());
      return it->second;
    }

    template <class BitVector, class Associations>
112
113
114
115
116
117
118
119
120
121
122
123
124
125
    static void periodicBC(Matrix& mat, Vector& sol, Vector& rhs, BitVector const& left, Associations const& left2right,
                           bool setDiagonal = true)
    {
      SymmetryStructure const symmetry = mat.symmetry();
      if (symmetry == SymmetryStructure::spd || symmetry == SymmetryStructure::symmetric || symmetry == SymmetryStructure::hermitian)
        symmetricPeriodicBC(mat.matrix(), sol.vector(), rhs.vector(), left, left2right, setDiagonal);
      else
        unsymmetricPeriodicBC(mat.matrix(), sol.vector(), rhs.vector(), left, left2right, setDiagonal);
    }


    template <class Mat, class Vec, class BitVector, class Associations>
    static void symmetricPeriodicBC(Mat& mat, Vec& sol, Vec& rhs, BitVector const& left, Associations const& left2right,
                                    bool setDiagonal = true)
126
127
128
    {
      error_exit("Not implemented");
    }
129
130


Praetorius, Simon's avatar
Praetorius, Simon committed
131
    template <class Value>
132
133
134
    struct Triplet
    {
      std::size_t row, col;
Praetorius, Simon's avatar
Praetorius, Simon committed
135
      Value value;
136
137
138
139
140
    };

    template <class Mat, class Vec, class BitVector, class Associations>
    static void unsymmetricPeriodicBC(Mat& mat, Vec& sol, Vec& rhs, BitVector const& left, Associations const& left2right,
                                      bool setDiagonal = true)
141
    {
142
      std::vector<Triplet<typename Mat::value_type>> rowValues;
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
      rowValues.reserve(left2right.size()*std::size_t(mat.nnz()/(0.9*num_rows(mat))));

      // Define the property maps
      // auto row   = mtl::mat::row_map(mat);
      auto col   = mtl::mat::col_map(mat);
      auto value = mtl::mat::value_map(mat);

      // iterate over the matrix
      std::size_t slotSize = 0;
      for (auto r : mtl::rows_of(mat)) {
        if (left[r.value()]) {
          slotSize = std::max(slotSize, std::size_t(mat.nnz_local(r.value())));
          std::size_t right = at(left2right,r.value());

          for (auto i : mtl::nz_of(r)) {
            rowValues.push_back({right,col(i),value(i)});
159
            value(i, 0);
160
161
162
163
          }
        }
      }

164
      mtl::mat::inserter<Mat, mtl::update_plus<typename Mat::value_type> > ins(mat, 2*slotSize);
165
166
167
      for (auto const& t : rowValues)
        ins[t.row][t.col] += t.value;

168
169
170
171
172
173
      for (std::size_t i = 0; i < mtl::size(left); ++i) {
        if (left[i]) {
          std::size_t j = at(left2right,i);
          if (setDiagonal) {
            ins[i][i] = 1;
            ins[i][j] = -1;
174
          }
175
176
177
178
179

          rhs[j] += rhs[i];
          rhs[i] = 0;

          sol[j] = sol[i];
180
181
182
        }
      }
    }
183

184
185
186
  };

} // end namespace AMDiS