Flow123d  master-f44eb46
local_system.cc
Go to the documentation of this file.
1 
2 #include "local_system.hh"
3 #include "la/local_constraint.hh"
4 
5 #include <armadillo>
6 #include "system/asserts.hh"
7 #include "system/index_types.hh"
8 
9 
11 {}
12 
13 
15 : row_dofs(nrows),
16  col_dofs(ncols),
17  matrix(nrows, ncols),
18  rhs(nrows),
19  sparsity(nrows,ncols),
20  elim_rows(nrows),
21  elim_cols(ncols),
22  solution_rows(nrows),
23  solution_cols(ncols),
24  diag_rows(nrows)
25 {
26  reset();
27 }
28 
29 void LocalSystem::set_size(uint nrows, uint ncols)
30 {
31  row_dofs.set_size(nrows);
32  col_dofs.set_size(ncols);
33  matrix.set_size(nrows, ncols);
34  rhs.set_size(nrows);
35  elim_rows.set_size(nrows);
36  elim_cols.set_size(ncols);
37  solution_rows.set_size(nrows);
38  solution_cols.set_size(ncols);
39  diag_rows.set_size(nrows);
40  // destroy previous sparsity pattern
41  sparsity.set_size(nrows,ncols);
42  sparsity.fill(almost_zero);
43 }
44 
45 
47 {
48  // zeros in local system
49  matrix.zeros();
50  rhs.zeros();
51  // drop all dirichlet values
53 }
54 
55 
56 void LocalSystem::reset(uint nrows, uint ncols)
57 {
58  if(matrix.n_rows != nrows || matrix.n_cols != ncols)
59  {
60  set_size(nrows, ncols);
61  }
62 
63  reset();
64 }
65 
66 
67 
68 void LocalSystem::reset(const LocDofVec &rdofs, const LocDofVec &cdofs)
69 {
70  if(matrix.n_rows != rdofs.n_rows || matrix.n_cols != cdofs.n_rows)
71  {
72  set_size(rdofs.n_rows, cdofs.n_rows);
73  }
74 
75  reset();
76  row_dofs = rdofs;
77  col_dofs = cdofs;
78 }
79 
80 
81 
82 void LocalSystem::set_solution(uint loc_dof, double solution, double diag)
83 {
84  // check that dofs are same
85  //ASSERT( arma::all(row_dofs == col_dofs) );
86  set_solution_row(loc_dof, solution, diag);
87  set_solution_col(loc_dof, solution);
88 }
89 
91 {
92  if (loc_constraint.type_ == ConstraintType::rows || loc_constraint.type_ == ConstraintType::both)
93  set_solution_row(loc_constraint.loc_dof_, loc_constraint.solution_, loc_constraint.diag_);
94  if (loc_constraint.type_ == ConstraintType::cols || loc_constraint.type_ == ConstraintType::both)
95  set_solution_col(loc_constraint.loc_dof_, loc_constraint.solution_);
96 }
97 
98 void LocalSystem::set_solution_row(uint loc_row, double solution, double diag) {
99  elim_rows[n_elim_rows]=loc_row;
100  solution_rows[n_elim_rows] = solution;
101  diag_rows[n_elim_rows] = diag;
102  n_elim_rows++;
103 }
104 
105 void LocalSystem::set_solution_col(uint loc_col, double solution) {
106  elim_cols[n_elim_cols]=loc_col;
107  solution_cols[n_elim_cols] = solution;
108  n_elim_cols++;
109 }
110 
112 {
113  //if there is solution set, eliminate:
114  if (n_elim_rows || n_elim_cols )
115  {
116  //DebugOut().fmt("elim rows: {} elim_cols: {}", n_elim_rows, n_elim_cols);
117 
118  arma::mat tmp_mat = matrix;
119  arma::vec tmp_rhs = rhs;
120 
121  uint ic, ir, row, col;
122 
123  // eliminate columns
124  for(ic=0; ic < n_elim_cols; ic++) {
125  col = elim_cols[ic];
126  tmp_rhs -= solution_cols[ic] * tmp_mat.col( col );
127  tmp_mat.col( col ).zeros();
128  }
129 
130  // eliminate rows
131  for(ir=0; ir < n_elim_rows; ir++) {
132  row = elim_rows[ir];
133  tmp_rhs( row ) = 0.0;
134  tmp_mat.row( row ).zeros();
135 
136  // fix global diagonal
137  for(ic=0; ic < n_elim_cols; ic++) {
138  col = elim_cols[ic];
139  if (row_dofs[row] == col_dofs[col]) {
140  ASSERT(fabs(solution_rows[ir] - solution_cols[ic]) <1e-12 );
141  // if preferred value is not set, then try using matrix value
142  double new_diagonal = matrix(row, col);
143 
144  if (diag_rows[ir] != 0.0) // if preferred value is set
145  new_diagonal = diag_rows[ir];
146  else if(new_diagonal == 0.0) // if an assembled value is not available
147  new_diagonal = 1.0;
148  // double new_diagonal = fabs(matrix(sol_row,col));
149  // if (new_diagonal == 0.0) {
150  // if (matrix.is_square()) {
151  // new_diagonal = arma::sum( abs(matrix.diag())) / matrix.n_rows;
152  // } else {
153  // new_diagonal = arma::accu( abs(matrix) ) / matrix.n_elem;
154  // }
155  // }
156  tmp_mat(row,col) = new_diagonal;
157  tmp_rhs(row) = new_diagonal * solution_rows[ir];
158 
159  }
160  }
161  }
162 
163  matrix = tmp_mat;
164  rhs = tmp_rhs;
166  }
167 
168  // filling almost_zero according to sparsity pattern
169  ASSERT_EQ(matrix.n_rows, sparsity.n_rows);
170  ASSERT_EQ(matrix.n_cols, sparsity.n_cols);
171  matrix = matrix + sparsity;
172 
173  //DebugOut() << matrix;
174  //DebugOut() << rhs;
175 }
176 
177 
178 void LocalSystem::add_value(uint row, uint col, double mat_val, double rhs_val)
179 {
180  ASSERT(row < matrix.n_rows);
181  ASSERT(col < matrix.n_cols);
182  ASSERT(sparsity(row,col))(row)(col).error("Violation of sparsity pattern.");
183 
184  matrix(row, col) += mat_val;
185  rhs(row) += rhs_val;
186 }
187 
188 void LocalSystem::add_value(uint row, uint col, double mat_val)
189 {
190  ASSERT(row < matrix.n_rows);
191  ASSERT(col < matrix.n_cols);
192  ASSERT(sparsity(row,col))(row)(col).error("Violation of sparsity pattern.");
193 
194  matrix(row, col) += mat_val;
195 }
196 
197 void LocalSystem::add_value(uint row, double rhs_val)
198 {
199  ASSERT(row < matrix.n_rows);
200 
201  rhs(row) += rhs_val;
202 }
203 
204 
206  ASSERT_EQ(matrix.n_rows, m.n_rows);
207  ASSERT_EQ(matrix.n_cols, m.n_cols);
208  matrix = m;
209 }
210 
212  ASSERT_EQ(matrix.n_rows, r.n_rows);
213  rhs = r;
214 }
215 
216 void LocalSystem::set_sparsity(const arma::umat & sp)
217 {
218  ASSERT_EQ(sparsity.n_rows, sp.n_rows);
219  ASSERT_EQ(sparsity.n_cols, sp.n_cols);
220 
221  sparsity.zeros();
222  for(uint i=0; i < sp.n_rows; i++)
223  for(uint j=0; j < sp.n_cols; j++)
224  if( sp(i,j) != 0 )
225  sparsity(i,j) = almost_zero;
226 // sparsity.print("sparsity");
227 }
228 
229 void LocalSystem::compute_schur_complement(uint offset, LocalSystem& schur, bool negative) const
230 {
231  // only for square matrix
232  ASSERT_EQ(matrix.n_rows, matrix.n_cols)("Cannot compute Schur complement for non-square matrix.");
233  arma::uword n = matrix.n_rows - 1;
234  ASSERT_LT(offset, n)("Schur complement (offset) dimension mismatch.");
235 
236  // B * invA
237  arma::mat BinvA = matrix.submat(offset, 0, n, offset-1) * matrix.submat(0, 0, offset-1, offset-1).i();
238 
239  // Schur complement S = C - B * invA * Bt
240  schur.matrix = matrix.submat(offset, offset, n, n) - BinvA * matrix.submat(0, offset, offset-1, n);
241  schur.rhs = rhs.subvec(offset, n) - BinvA * rhs.subvec(0, offset-1);
242 
243  if(negative)
244  {
245  schur.matrix = -1.0 * schur.matrix;
246  schur.rhs = -1.0 * schur.rhs;
247  }
248 }
249 
250 void LocalSystem::reconstruct_solution_schur(uint offset, const arma::vec &schur_solution, arma::vec& reconstructed_solution) const
251 {
252  // only for square matrix
253  ASSERT_EQ(matrix.n_rows, matrix.n_cols)("Cannot compute Schur complement for non-square matrix.");
254  arma::uword n = matrix.n_rows - 1;
255  ASSERT_LT(offset, n)("Schur complement (offset) dimension mismatch.");
256 
257  reconstructed_solution.set_size(offset);
258  // invA
259  arma::mat invA = matrix.submat(0, 0, offset-1, offset-1).i();
260 
261  // x = invA*b - invA * Bt * schur_solution
262  reconstructed_solution = invA * rhs.subvec(0,offset-1) - invA * matrix.submat(0, offset, offset-1, n) * schur_solution;
263 }
Definitions of ASSERTS.
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:301
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
ConstraintType type_
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()
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)
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