OGS
anonymous_namespace{EigenTools.cpp} Namespace Reference

Functions

void setRowsToZeroOffDiagonal (MathLib::EigenMatrix::RawMatrixType &A_eigen, const std::vector< MathLib::EigenMatrix::IndexType > &row_ids)
void applyKnownSolutionComplete (MathLib::EigenMatrix::RawMatrixType &A_eigen, MathLib::EigenVector::RawVectorType &b_eigen, const std::vector< MathLib::EigenMatrix::IndexType > &vec_knownX_id, const std::vector< double > &vec_knownX_x)
void applyKnownSolutionIncomplete (MathLib::EigenMatrix::RawMatrixType &A_eigen, MathLib::EigenVector::RawVectorType &b_eigen, const std::vector< MathLib::EigenMatrix::IndexType > &vec_knownX_id, const std::vector< double > &vec_knownX_x)

Function Documentation

◆ applyKnownSolutionComplete()

void anonymous_namespace{EigenTools.cpp}::applyKnownSolutionComplete ( MathLib::EigenMatrix::RawMatrixType & A_eigen,
MathLib::EigenVector::RawVectorType & b_eigen,
const std::vector< MathLib::EigenMatrix::IndexType > & vec_knownX_id,
const std::vector< double > & vec_knownX_x )

Definition at line 39 of file EigenTools.cpp.

44{
45 setRowsToZeroOffDiagonal(A_eigen, vec_knownX_id);
46
48 SpMat AT = A_eigen.transpose();
49
50 // Reserve space for at least one value (on the diagonal). For
51 // deactivated subdomains some rows and columns might end empty (in the
52 // A matrix) and so no space is reserved for the Dirichlet conditions in
53 // the transposed matrix. Then the coeffRef call will do costly
54 // reallocations.
55 AT.reserve(Eigen::VectorXi::Constant(A_eigen.rows(), 1));
56
57 for (std::size_t ix = 0; ix < vec_knownX_id.size(); ix++)
58 {
59 SpMat::Index const row_id = vec_knownX_id[ix];
60 auto const x = vec_knownX_x[ix];
61
62 // b_i -= A_eigen(i,k)*val, i!=k
63 // set column to zero, subtract from rhs
64 for (SpMat::InnerIterator it(AT, row_id); it; ++it)
65 {
66 if (it.col() == row_id)
67 {
68 continue;
69 }
70
71 b_eigen[it.col()] -= it.value() * x;
72 it.valueRef() = 0.0;
73 }
74
75 auto& c = AT.coeffRef(row_id, row_id);
76 if (c != 0.0)
77 {
78 b_eigen[row_id] = x * c;
79 }
80 else
81 {
82 b_eigen[row_id] = x;
83 c = 1.0;
84 }
85 }
86
87 A_eigen = AT.transpose();
88}
Eigen::SparseMatrix< double, Eigen::RowMajor > RawMatrixType
Definition EigenMatrix.h:31
void setRowsToZeroOffDiagonal(MathLib::EigenMatrix::RawMatrixType &A_eigen, const std::vector< MathLib::EigenMatrix::IndexType > &row_ids)

References setRowsToZeroOffDiagonal().

◆ applyKnownSolutionIncomplete()

void anonymous_namespace{EigenTools.cpp}::applyKnownSolutionIncomplete ( MathLib::EigenMatrix::RawMatrixType & A_eigen,
MathLib::EigenVector::RawVectorType & b_eigen,
const std::vector< MathLib::EigenMatrix::IndexType > & vec_knownX_id,
const std::vector< double > & vec_knownX_x )

Definition at line 90 of file EigenTools.cpp.

95{
96 INFO("partial Dirichlet BC application (rhs only).");
97
98 setRowsToZeroOffDiagonal(A_eigen, vec_knownX_id);
99
100 Eigen::VectorXd x_known{A_eigen.rows()};
101 x_known.setZero();
102
103 for (auto const& [row_id, x] :
104 ranges::views::zip(vec_knownX_id, vec_knownX_x))
105 {
106 x_known[row_id] = x;
107
108 // set rhs to Dirichlet value
109 auto const c = A_eigen.coeff(row_id, row_id);
110 if (c != 0)
111 {
112 b_eigen[row_id] = x * c;
113
114 // exclude diagonal (A_eigen(i,i)) from multiplication below
115 // (A_eigen * x_known)
116 A_eigen.coeffRef(row_id, row_id) = 0;
117 }
118 else
119 {
120 b_eigen[row_id] = x;
121 }
122 }
123
124 b_eigen -= A_eigen * x_known;
125}
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:36

References INFO(), and setRowsToZeroOffDiagonal().

◆ setRowsToZeroOffDiagonal()

void anonymous_namespace{EigenTools.cpp}::setRowsToZeroOffDiagonal ( MathLib::EigenMatrix::RawMatrixType & A_eigen,
const std::vector< MathLib::EigenMatrix::IndexType > & row_ids )

Definition at line 19 of file EigenTools.cpp.

22{
24
25 // A_eigen(k, j) = 0.
26 // set row to zero
27 for (auto row_id : row_ids)
28 {
29 for (SpMat::InnerIterator it(A_eigen, row_id); it; ++it)
30 {
31 if (it.col() != decltype(it.col())(row_id))
32 {
33 it.valueRef() = 0.0;
34 }
35 }
36 }
37}

Referenced by applyKnownSolutionComplete(), and applyKnownSolutionIncomplete().