Flow123d  master-f44eb46
local_system.hh
Go to the documentation of this file.
1 
2 #ifndef LOCAL_SYSTEM_HH_
3 #define LOCAL_SYSTEM_HH_
4 
5 #include <armadillo>
6 #include "system/index_types.hh"
7 
8 class LinSys;
9 class LocalConstraint;
10 
11 /** Local system class is meant to be used for local assembly and then pass to global linear system.
12  * The key idea is to take care of known solution values (Dirichlet boundary conditions) in a common way.
13  *
14  * A] Usage of the class consists of 5 steps:
15  * 1) create local system, set local/global DoFs.
16  * 2) possibly setup sparsity pattern for the local system
17  * 3) set matrix and RHS entries
18  * 4) set all known values (Dirichlet BC)
19  * (the order of items 3 and 4 does not matter)
20  * 5) possibly eliminate the known solution
21  * (possibly fix the diagonal entries of the local system, where Dirichlet BC is set)
22  * 6) pass the local system to the global system, possibly with the map local_to_global
23  *
24  * B] Scenario with computation of Schur complement:
25  * 1) do A1, do not have to set dof indices (if not doing B3)
26  * 2) do A2, A3
27  * 3) possibly do A4, A5, if there is known solution for the dofs not included in the Schur complement
28  * 3) create another LocalSystem for Schur complement (local_schur), possibly set dof indices
29  * 4) possibly do A2 on local_schur
30  * 4) compute Schur complement, passing the prepared local_schur
31  * 5) if not done before, set dof indices
32  * 6) possibly do A4, A5 on local_schur
33  * 7) do A6 with local_schur
34  *
35  * C] Reconstruction from the Schur complement
36  * 1) do A1, A3
37  * 2) reconstruct the full solution
38  *
39  * TODO:
40  * - set dofs as references
41  * - done: use arma::vec<int> to keep dofs (efficient, can be constructed from raw arrays)
42  * - done: use local dof indeces to set solution
43  * - rename set_solution to set_constraint
44  * -
45  */
47 {
48 public:
49  /**
50  * Global row and col indices. Are public and can be freely set.
51  * Nevertheless one can also provide reference to already existing arrays through
52  * specific constructor or reset function.
53  */
55 
56  /**
57  * @brief Default constructor.
58  *
59  * Object must be initialized by subsequent call of reset(nrows, ncols).
60  */
61  LocalSystem();
62 
63  /** @brief Constructor.
64  *
65  * @p nrows is number of rows of local system
66  * @p ncols is number of columns of local system
67  */
68  LocalSystem(uint nrows, uint ncols);
69 
70  /// Resets the matrix, RHS, dofs to zero and clears solution settings
71  void reset();
72 
73  /// Resize and reset.
74  void reset(uint nrows, uint ncols);
75 
76  /**
77  * Resize and reset. Set dofs vectors to reuse arrays provided by given vectors.
78  * Given vectors can not be changed until next call to of any reset function.
79  */
80  void reset(const LocDofVec &row_dofs, const LocDofVec &col_dofs);
81 
82  const arma::mat& get_matrix() {return matrix;}
83  const arma::vec& get_rhs() {return rhs;}
84 
85  /** @brief Set the position and value of known solution. E.g. Dirichlet boundary condition.
86  *
87  * @p loc_dofs is local row index in solution vector
88  * @p solution is the values of the solution
89  * @p diag_val is preferred diagonal value on the solution row
90  */
91  void set_solution(uint loc_dof, double solution, double diag=0.0);
92 
93  /** @brief Set the position and value of known solution. E.g. Dirichlet boundary condition.
94  *
95  * @p LocalConstraint data object holds all needed data used in previous method
96  */
97  void set_solution(LocalConstraint &loc_constraint);
98 
99  void set_solution_row(uint loc_row, double solution, double diag=0.0);
100 
101  void set_solution_col(uint loc_col, double solution);
102 
103 
104  /**
105  * When finished with assembly of the local system,
106  * this function eliminates all the known dofs.
107  *
108  * It is skipped if there is not any solution dof set.
109  *
110  * During elimination, the (global) diagonal entries on the rows, where the solution is set, might be zero.
111  * Therefore it is necessary to set a proper value to the diagonal entry
112  * and respective RHS entry, such that the given solution holds.
113  * If preferred diagonal value has been set by @p set_solution then it is used.
114  *
115  * Calling this function after the assembly of local system
116  * and before passing the local system to the global one
117  * is finished is users's responsibility.
118  */
119  void eliminate_solution();
120 
121  /** @brief Adds a single entry into the local system.
122  *
123  * @p row is local row index of local system
124  * @p col is local column index of local system
125  * @p mat_val is matrix entry value
126  * @p rhs_val is RHS entry value
127  */
128  void add_value(uint row, uint col,
129  double mat_val, double rhs_val);
130 
131  /** @brief Matrix entry.
132  * Adds a single entry into the local system matrix.
133  *
134  * @p row is local row index of local system
135  * @p col is local column index of local system
136  * @p mat_val is matrix entry value
137  */
138  void add_value(uint row, uint col,
139  double mat_val);
140 
141  /** @brief RHS entry.
142  * Adds a single entry into the local system RHS.
143  *
144  * @p row is local row index of local system
145  * @p rhs_val is RHS entry value
146  */
147  void add_value(uint row, double rhs_val);
148 
150  void set_rhs(arma::vec rhs);
151 
152  /// Sets the sparsity pattern for the local system.
153  /** Due to petsc options: MatSetOption(matrix_, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE)
154  * all zeros will be thrown away from the system.
155  * If we do not want some zero entries in the system matrix to be thrown away,
156  * we can set these entries with this almost zero value.
157  *
158  * Almost_zero values will be set in all entries: sp(i,j) != 0
159  */
160  void set_sparsity(const arma::umat & sp);
161 
162  /** @brief Computes Schur complement of the local system: S = C - B * invA * Bt
163  * Applicable for square matrices.
164  * It can be called either after eliminating Dirichlet dofs,
165  * or the Dirichlet dofs can be set on the Schur complement
166  * and the elimination done on the Schur complement.
167  *
168  * @p offset index of the first row/column of submatrix C (size of A)
169  * @p schur (output) LocalSystem with Schur complement
170  * @p negative if true, the schur complement (including its rhs) is multiplied by -1.0
171  */
172  void compute_schur_complement(uint offset, LocalSystem& schur, bool negative=false) const;
173 
174  /** @brief Reconstructs the solution from the Schur complement solution: x = invA*b - invA * Bt * schur_solution
175  * Applicable for square matrices.
176  *
177  * @p offset index of the first row/column of submatrix C (size of A)
178  * @p schur_solution solution of the Schur complement
179  * @p reconstructed_solution (output) reconstructed solution of the complementary variable
180  */
181  void reconstruct_solution_schur(uint offset, const arma::vec &schur_solution, arma::vec& reconstructed_solution) const;
182 
183  /// Due to petsc options: MatSetOption(matrix_, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE)
184  /// all zeros will be thrown away from the system.
185  /// If we do not want some zero entries in the system matrix to be thrown away,
186  /// we can set these entries with this almost zero value.
187  ///
188  /// This is done for example when BC values are eliminated and later the BC changes to different type (e.g. seepage).
189  /// Another case is keeping the structure of matrix unchanged for the schur complements -
190  /// for that we fill the whole diagonal (escpecially block C in darcy flow) with artificial zeros.
191  static constexpr double almost_zero = std::numeric_limits<double>::min();
192 
193 
194 protected:
195  void set_size(uint nrows, uint ncols);
196 
197  arma::mat matrix; ///< local system matrix
198  arma::vec rhs; ///< local system RHS
199 
200  arma::mat sparsity; ///< sparsity pattern
201 
202  /// Number of rows/cols to be eliminated due to known solution.
204  LocDofVec elim_rows; /// Rows indices of rows to be eliminated.
205  LocDofVec elim_cols; /// Cols indices of cols to be eliminated.
206  arma::vec solution_rows; /// Values of the known solution (for row dofs).
207  arma::vec solution_cols; /// Values of the known solution (for col dofs).
208  arma::vec diag_rows; /// Prefered values on the diagonal after elimination.
209 
210 
211 
212  friend class LinSys;
213 };
214 
215 #endif // LOCAL_SYSTEM_HH_
void set_sparsity(const arma::umat &sp)
Sets the sparsity pattern for the local system.
LocDofVec elim_cols
Rows indices of rows to be eliminated.
void set_solution_col(uint loc_col, double solution)
LocalSystem()
Default constructor.
Definition: local_system.cc:10
void set_solution(uint loc_dof, double solution, double diag=0.0)
Set the position and value of known solution. E.g. Dirichlet boundary condition.
Definition: local_system.cc:82
LocDofVec col_dofs
Definition: local_system.hh:54
void set_size(uint nrows, uint ncols)
Definition: local_system.cc:29
uint n_elim_rows
Number of rows/cols to be eliminated due to known solution.
void eliminate_solution()
const arma::mat & get_matrix()
Definition: local_system.hh:82
void reconstruct_solution_schur(uint offset, const arma::vec &schur_solution, arma::vec &reconstructed_solution) const
Reconstructs the solution from the Schur complement solution: x = invA*b - invA * Bt * schur_solution...
arma::mat sparsity
sparsity pattern
void set_rhs(arma::vec rhs)
const arma::vec & get_rhs()
Definition: local_system.hh:83
arma::vec rhs
local system RHS
void add_value(uint row, uint col, double mat_val, double rhs_val)
Adds a single entry into the local system.
void set_matrix(arma::mat matrix)
arma::vec solution_rows
Cols indices of cols to be eliminated.
arma::vec solution_cols
Values of the known solution (for row dofs).
arma::vec diag_rows
Values of the known solution (for col dofs).
static constexpr double almost_zero
LocDofVec row_dofs
Definition: local_system.hh:54
void compute_schur_complement(uint offset, LocalSystem &schur, bool negative=false) const
Computes Schur complement of the local system: S = C - B * invA * Bt Applicable for square matrices....
arma::mat matrix
local system matrix
LocDofVec elim_rows
void reset()
Resets the matrix, RHS, dofs to zero and clears solution settings.
Definition: local_system.cc:46
void set_solution_row(uint loc_row, double solution, double diag=0.0)
Definition: local_system.cc:98
arma::Col< IntIdx > LocDofVec
Definition: index_types.hh:28
unsigned int uint
ArmaMat< double, N, M > mat
Definition: armor.hh:888
ArmaVec< double, N > vec
Definition: armor.hh:885