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> 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> 24 template <
class BitVector>
25 static void dirichletBC(Matrix& mat,
Vector& sol,
Vector& rhs, BitVector
const& nodes,
bool setDiagonal =
true)
27 SymmetryStructure
const symmetry = mat.
symmetry();
28 if (symmetry == SymmetryStructure::spd || symmetry == SymmetryStructure::symmetric || symmetry == SymmetryStructure::hermitian)
34 template <
class Mat,
class Vec,
class BitVector>
35 static void symmetricDirichletBC(Mat& mat, Vec& sol, Vec& rhs, BitVector
const& nodes,
bool setDiagonal =
true)
38 auto row = mtl::mat::row_map(mat);
39 auto col = mtl::mat::col_map(mat);
40 auto value = mtl::mat::value_map(mat);
43 std::size_t slotSize = 0;
44 for (
auto r : mtl::rows_of(mat)) {
45 std::size_t rowSize = 0;
46 for (
auto i : mtl::nz_of(r)) {
52 else if (setDiagonal && nodes[col(i)]) {
53 rhs[row(i)] -= value(i) * sol[col(i)];
57 slotSize = std::max(slotSize, rowSize);
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) {
70 for (
typename Vec::size_type i = 0; i < mtl::size(sol); ++i) {
76 template <
class Mat,
class Vec,
class BitVector>
77 static void unsymmetricDirichletBC(Mat& mat, Vec& sol, Vec& rhs, BitVector
const& nodes,
bool setDiagonal =
true)
80 auto row = mtl::mat::row_map(mat);
81 auto col = mtl::mat::col_map(mat);
82 auto value = mtl::mat::value_map(mat);
85 for (
auto r : mtl::rows_of(mat)) {
86 if (nodes[r.value()]) {
87 for (
auto i : mtl::nz_of(r)) {
89 value(i, (setDiagonal && row(i) == col(i) ? 1 : 0) );
95 for (
typename Vec::size_type i = 0; i < mtl::size(sol); ++i) {
103 template <
class Associations>
104 static std::size_t at(Associations
const& m, std::size_t idx)
106 auto it = m.find(idx);
107 assert(it != m.end());
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)
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);
119 unsymmetricPeriodicBC(mat.
matrix(), sol.
vector(), rhs.
vector(), left, left2right, setDiagonal);
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)
131 template <
class Value>
134 std::size_t row, col;
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)
142 std::vector<Triplet<typename Mat::value_type>> rowValues;
143 rowValues.reserve(left2right.size()*std::size_t(mat.nnz()/(0.9*num_rows(mat))));
147 auto col = mtl::mat::col_map(mat);
148 auto value = mtl::mat::value_map(mat);
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());
157 for (
auto i : mtl::nz_of(r)) {
158 rowValues.push_back({right,col(i),value(i)});
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;
168 for (std::size_t i = 0; i < mtl::size(left); ++i) {
170 std::size_t j = at(left2right,i);
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