OGS
EigenTools.cpp
Go to the documentation of this file.
1
11#include "EigenTools.h"
12
13#include "EigenVector.h"
14
15namespace MathLib
16{
18 EigenMatrix& A, EigenVector& b, EigenVector& /*x*/,
19 const std::vector<EigenMatrix::IndexType>& vec_knownX_id,
20 const std::vector<double>& vec_knownX_x)
21{
22 using SpMat = EigenMatrix::RawMatrixType;
23 static_assert(SpMat::IsRowMajor, "matrix is assumed to be row major!");
24
25 auto& A_eigen = A.getRawMatrix();
26 auto& b_eigen = b.getRawVector();
27
28 // A_eigen(k, j) = 0.
29 // set row to zero
30 for (auto row_id : vec_knownX_id)
31 {
32 for (SpMat::InnerIterator it(A_eigen, row_id); it; ++it)
33 {
34 if (it.col() != decltype(it.col())(row_id))
35 {
36 it.valueRef() = 0.0;
37 }
38 }
39 }
40
41 SpMat AT = A_eigen.transpose();
42
43 // Reserve space for at least one value (on the diagonal). For deactivated
44 // subdomains some rows and columns might end empty (in the A matrix) and so
45 // no space is reserved for the Dirichlet conditions in the transposed
46 // matrix. Then the coeffRef call will do costly reallocations.
47 AT.reserve(Eigen::VectorXi::Constant(A_eigen.rows(), 1));
48
49 for (std::size_t ix = 0; ix < vec_knownX_id.size(); ix++)
50 {
51 SpMat::Index const row_id = vec_knownX_id[ix];
52 auto const x = vec_knownX_x[ix];
53
54 // b_i -= A_eigen(i,k)*val, i!=k
55 // set column to zero, subtract from rhs
56 for (SpMat::InnerIterator it(AT, row_id); it; ++it)
57 {
58 if (it.col() == row_id)
59 {
60 continue;
61 }
62
63 b_eigen[it.col()] -= it.value() * x;
64 it.valueRef() = 0.0;
65 }
66
67 auto& c = AT.coeffRef(row_id, row_id);
68 if (c != 0.0)
69 {
70 b_eigen[row_id] = x * c;
71 }
72 else
73 {
74 b_eigen[row_id] = x;
75 c = 1.0;
76 }
77 }
78
79 A_eigen = AT.transpose();
80}
81
82} // namespace MathLib
Eigen::SparseMatrix< double, Eigen::RowMajor > RawMatrixType
Definition EigenMatrix.h:31
RawMatrixType & getRawMatrix()
Global vector based on Eigen vector.
Definition EigenVector.h:25
RawVectorType & getRawVector()
return a raw Eigen vector object
void applyKnownSolution(EigenMatrix &A, EigenVector &b, EigenVector &, const std::vector< EigenMatrix::IndexType > &vec_knownX_id, const std::vector< double > &vec_knownX_x)