AMDiS  0.3
The Adaptive Multi-Dimensional Simulation Toolbox
Constraints.hpp
1 #pragma once
2 
3 #include <algorithm>
4 #include <vector>
5 
6 #include <boost/numeric/mtl/matrix/compressed2D.hpp>
7 #include <boost/numeric/mtl/matrix/inserter.hpp>
8 #include <boost/numeric/mtl/utility/property_map.hpp>
9 #include <boost/numeric/mtl/utility/range_wrapper.hpp>
10 
11 #include <amdis/linearalgebra/Constraints.hpp>
12 #include <amdis/linearalgebra/SymmetryStructure.hpp>
13 #include <amdis/linearalgebra/mtl/MatrixBackend.hpp>
14 #include <amdis/linearalgebra/mtl/VectorBackend.hpp>
15 
16 namespace AMDiS
17 {
18  template <class T>
20  {
21  using Matrix = MTLSparseMatrix<T>;
22  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  SymmetryStructure const symmetry = mat.symmetry();
28  if (symmetry == SymmetryStructure::spd || symmetry == SymmetryStructure::symmetric || symmetry == SymmetryStructure::hermitian)
29  symmetricDirichletBC(mat.matrix(), sol.vector(), rhs.vector(), nodes, setDiagonal);
30  else
31  unsymmetricDirichletBC(mat.matrix(), sol.vector(), rhs.vector(), nodes, setDiagonal);
32  }
33 
34  template <class Mat, class Vec, class BitVector>
35  static void symmetricDirichletBC(Mat& mat, Vec& sol, Vec& rhs, BitVector const& nodes, bool setDiagonal = true)
36  {
37  // Define the property maps
38  auto row = mtl::mat::row_map(mat);
39  auto col = mtl::mat::col_map(mat);
40  auto value = mtl::mat::value_map(mat);
41 
42  // iterate over the matrix
43  std::size_t slotSize = 0;
44  for (auto r : mtl::rows_of(mat)) { // rows or columns
45  std::size_t rowSize = 0;
46  for (auto i : mtl::nz_of(r)) { // non-zeros within
47  ++rowSize;
48  if (nodes[row(i)]) {
49  // set identity row
50  value(i, 0);
51  }
52  else if (setDiagonal && nodes[col(i)]) {
53  rhs[row(i)] -= value(i) * sol[col(i)];
54  value(i, 0);
55  }
56  }
57  slotSize = std::max(slotSize, rowSize);
58  }
59 
60  // set diagonal entry
61  if (setDiagonal) {
62  mtl::mat::inserter<Mat, mtl::update_store<typename Mat::value_type> > ins(mat, slotSize);
63  for (std::size_t i = 0; i < nodes.size(); ++i) {
64  if (nodes[i])
65  ins[i][i] = 1;
66  }
67  }
68 
69  // copy solution dirichlet data to rhs vector
70  for (typename Vec::size_type i = 0; i < mtl::size(sol); ++i) {
71  if (nodes[i])
72  rhs[i] = sol[i];
73  }
74  }
75 
76  template <class Mat, class Vec, class BitVector>
77  static void unsymmetricDirichletBC(Mat& mat, Vec& sol, Vec& rhs, BitVector const& nodes, bool setDiagonal = true)
78  {
79  // Define the property maps
80  auto row = mtl::mat::row_map(mat);
81  auto col = mtl::mat::col_map(mat);
82  auto value = mtl::mat::value_map(mat);
83 
84  // iterate over the matrix
85  for (auto r : mtl::rows_of(mat)) { // rows of the matrix
86  if (nodes[r.value()]) {
87  for (auto i : mtl::nz_of(r)) { // non-zeros within
88  // set identity row
89  value(i, (setDiagonal && row(i) == col(i) ? 1 : 0) );
90  }
91  }
92  }
93 
94  // copy solution dirichlet data to rhs vector
95  for (typename Vec::size_type i = 0; i < mtl::size(sol); ++i) {
96  if (nodes[i])
97  rhs[i] = sol[i];
98  }
99  }
100 
101 
102 
103  template <class Associations>
104  static std::size_t at(Associations const& m, std::size_t idx)
105  {
106  auto it = m.find(idx);
107  assert(it != m.end());
108  return it->second;
109  }
110 
111  template <class BitVector, class Associations>
112  static void periodicBC(Matrix& mat, Vector& sol, Vector& rhs, BitVector const& left, Associations const& left2right,
113  bool setDiagonal = true)
114  {
115  SymmetryStructure const symmetry = mat.symmetry();
116  if (symmetry == SymmetryStructure::spd || symmetry == SymmetryStructure::symmetric || symmetry == SymmetryStructure::hermitian)
117  symmetricPeriodicBC(mat.matrix(), sol.vector(), rhs.vector(), left, left2right, setDiagonal);
118  else
119  unsymmetricPeriodicBC(mat.matrix(), sol.vector(), rhs.vector(), left, left2right, setDiagonal);
120  }
121 
122 
123  template <class Mat, class Vec, class BitVector, class Associations>
124  static void symmetricPeriodicBC(Mat& mat, Vec& sol, Vec& rhs, BitVector const& left, Associations const& left2right,
125  bool setDiagonal = true)
126  {
127  error_exit("Not implemented");
128  }
129 
130 
131  template <class Value>
132  struct Triplet
133  {
134  std::size_t row, col;
135  Value value;
136  };
137 
138  template <class Mat, class Vec, class BitVector, class Associations>
139  static void unsymmetricPeriodicBC(Mat& mat, Vec& sol, Vec& rhs, BitVector const& left, Associations const& left2right,
140  bool setDiagonal = true)
141  {
142  std::vector<Triplet<typename Mat::value_type>> rowValues;
143  rowValues.reserve(left2right.size()*std::size_t(mat.nnz()/(0.9*num_rows(mat))));
144 
145  // Define the property maps
146  // auto row = mtl::mat::row_map(mat);
147  auto col = mtl::mat::col_map(mat);
148  auto value = mtl::mat::value_map(mat);
149 
150  // iterate over the matrix
151  std::size_t slotSize = 0;
152  for (auto r : mtl::rows_of(mat)) {
153  if (left[r.value()]) {
154  slotSize = std::max(slotSize, std::size_t(mat.nnz_local(r.value())));
155  std::size_t right = at(left2right,r.value());
156 
157  for (auto i : mtl::nz_of(r)) {
158  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  for (auto const& t : rowValues)
166  ins[t.row][t.col] += t.value;
167 
168  for (std::size_t i = 0; i < mtl::size(left); ++i) {
169  if (left[i]) {
170  std::size_t j = at(left2right,i);
171  if (setDiagonal) {
172  ins[i][i] = 1;
173  ins[i][j] = -1;
174  }
175 
176  rhs[j] += rhs[i];
177  rhs[i] = 0;
178 
179  sol[j] = sol[i];
180  }
181  }
182  }
183 
184  };
185 
186 } // end namespace AMDiS
BaseMatrix & matrix()
Return a reference to the data-matrix matrix.
Definition: MatrixBackend.hpp:46
SymmetryStructure symmetry() const
Symmetry of the matrix entries.
Definition: MatrixBackend.hpp:131
void error_exit(std::string const &str, Args &&... args)
print a message and exit
Definition: Output.hpp:142
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
The basic container that stores a base matrix.
Definition: MatrixBackend.hpp:21
The basic container that stores a base vector data.
Definition: VectorBackend.hpp:16
Definition: Constraints.hpp:13
BaseVector const & vector() const
Return the data-vector vector_.
Definition: VectorBackend.hpp:42