Flow123d  JS_before_hm-62-ge56f9d5
local_system.cc
Go to the documentation of this file.
1 
2 #include "local_system.hh"
3 
4 #include <armadillo>
5 #include "system/sys_vector.hh"
6 
7 
9 {}
10 
11 
12 LocalSystem::LocalSystem(unsigned int nrows, unsigned int ncols)
13 : row_dofs(nrows),
14  col_dofs(ncols),
15  matrix(nrows, ncols),
16  rhs(nrows),
17  sparsity(nrows,ncols),
18  elim_rows(nrows),
19  elim_cols(ncols),
20  solution_rows(nrows),
21  solution_cols(ncols),
22  diag_rows(nrows)
23 {
24  reset();
25 }
26 
27 void LocalSystem::set_size(unsigned int nrows, unsigned int ncols)
28 {
29  row_dofs.set_size(nrows);
30  col_dofs.set_size(ncols);
31  matrix.set_size(nrows, ncols);
32  rhs.set_size(nrows);
33  elim_rows.set_size(nrows);
34  elim_cols.set_size(ncols);
35  solution_rows.set_size(nrows);
36  solution_cols.set_size(ncols);
37  diag_rows.set_size(nrows);
38  // destroy previous sparsity pattern
39  sparsity.zeros(nrows,ncols);
40 }
41 
42 
44 {
45  // zeros in local system
46  matrix.zeros();
47  rhs.zeros();
48  // drop all dirichlet values
50 }
51 
52 
53 void LocalSystem::reset(arma::uword nrows, arma::uword ncols)
54 {
55  matrix.set_size(nrows, ncols);
56  rhs.set_size(nrows);
57  row_dofs.resize(matrix.n_rows);
58  col_dofs.resize(matrix.n_cols);
59  // destroy previous sparsity pattern
60  sparsity.zeros(nrows,ncols);
61  reset();
62 }
63 
64 
65 
66 void LocalSystem::reset(const DofVec &rdofs, const DofVec &cdofs)
67 {
68  set_size(rdofs.n_rows, cdofs.n_rows);
69  reset();
70  row_dofs = rdofs;
71  col_dofs = cdofs;
72 }
73 
74 
75 
76 void LocalSystem::set_solution(unsigned int loc_dof, double solution, double diag)
77 {
78  // check that dofs are same
79  //ASSERT_DBG( arma::all(row_dofs == col_dofs) );
80  set_solution_row(loc_dof, solution, diag);
81  set_solution_col(loc_dof, solution);
82 }
83 
84 void LocalSystem::set_solution_row(uint loc_row, double solution, double diag) {
85  elim_rows[n_elim_rows]=loc_row;
86  solution_rows[n_elim_rows] = solution;
87  diag_rows[n_elim_rows] = diag;
88  n_elim_rows++;
89 }
90 
91 void LocalSystem::set_solution_col(uint loc_col, double solution) {
92  elim_cols[n_elim_cols]=loc_col;
93  solution_cols[n_elim_cols] = solution;
94  n_elim_cols++;
95 }
96 
97 /*
98 
99 void LocalSystem::set_solution(const DofVec & loc_rows, const arma::vec &solution, const arma::vec &diag) {
100  ASSERT_EQ_DBG(loc_rows.n_rows(), solution)
101  set_solution_rows(loc_rows, solution, diag);
102  set_solution_cols()
103 }
104 void LocalSystem::set_solution_rows(DofVec & loc_rows, const arma::vec &solution, const arma::vec &diag);
105 void LocalSystem::set_solution_cols(DofVec & loc_cols, const arma::vec &solution);
106 
107 */
108 
110 {
111  //if there is solution set, eliminate:
112  if (n_elim_rows || n_elim_cols )
113  {
114  //DebugOut().fmt("elim rows: {} elim_cols: {}", n_elim_rows, n_elim_cols);
115 
116  arma::mat tmp_mat = matrix;
117  arma::vec tmp_rhs = rhs;
118 
119  unsigned int ic, ir, row, col;
120 
121  // eliminate columns
122  for(ic=0; ic < n_elim_cols; ic++) {
123  col = elim_cols[ic];
124  tmp_rhs -= solution_cols[ic] * tmp_mat.col( col );
125  tmp_mat.col( col ).zeros();
126  }
127 
128  // eliminate rows
129  for(ir=0; ir < n_elim_rows; ir++) {
130  row = elim_rows[ir];
131  tmp_rhs( row ) = 0.0;
132  tmp_mat.row( row ).zeros();
133 
134  // fix global diagonal
135  for(ic=0; ic < n_elim_cols; ic++) {
136  col = elim_cols[ic];
137  if (row_dofs[row] == col_dofs[col]) {
138  ASSERT_DBG(fabs(solution_rows[ir] - solution_cols[ic]) <1e-12 );
139  // if preferred value is not set, then try using matrix value
140  double new_diagonal = matrix(row, col);
141 
142  if (diag_rows[ir] != 0.0) // if preferred value is set
143  new_diagonal = diag_rows[ir];
144  else if(new_diagonal == 0.0) // if an assembled value is not available
145  new_diagonal = 1.0;
146  // double new_diagonal = fabs(matrix(sol_row,col));
147  // if (new_diagonal == 0.0) {
148  // if (matrix.is_square()) {
149  // new_diagonal = arma::sum( abs(matrix.diag())) / matrix.n_rows;
150  // } else {
151  // new_diagonal = arma::accu( abs(matrix) ) / matrix.n_elem;
152  // }
153  // }
154  tmp_mat(row,col) = new_diagonal;
155  tmp_rhs(row) = new_diagonal * solution_rows[ir];
156 
157  }
158  }
159  }
160 
161  matrix = tmp_mat;
162  rhs = tmp_rhs;
163  n_elim_cols=n_elim_rows=0;
164  }
165 
166  // filling almost_zero according to sparsity pattern
167  ASSERT_EQ_DBG(matrix.n_rows, sparsity.n_rows);
168  ASSERT_EQ_DBG(matrix.n_cols, sparsity.n_cols);
169  matrix = matrix + sparsity;
170 
171  //DebugOut() << matrix;
172  //DebugOut() << rhs;
173 }
174 
175 
176 void LocalSystem::add_value(unsigned int row, unsigned int col, double mat_val, double rhs_val)
177 {
178  ASSERT_DBG(row < matrix.n_rows);
179  ASSERT_DBG(col < matrix.n_cols);
180  ASSERT_DBG(sparsity(row,col))(row)(col).error("Violation of sparsity pattern.");
181 
182  matrix(row, col) += mat_val;
183  rhs(row) += rhs_val;
184 }
185 
186 void LocalSystem::add_value(unsigned int row, unsigned int col, double mat_val)
187 {
188  ASSERT_DBG(row < matrix.n_rows);
189  ASSERT_DBG(col < matrix.n_cols);
190  ASSERT_DBG(sparsity(row,col))(row)(col).error("Violation of sparsity pattern.");
191 
192  matrix(row, col) += mat_val;
193 }
194 
195 void LocalSystem::add_value(unsigned int row, double rhs_val)
196 {
197  ASSERT_DBG(row < matrix.n_rows);
198 
199  rhs(row) += rhs_val;
200 }
201 
202 
204  ASSERT_EQ_DBG(matrix.n_rows, m.n_rows);
205  ASSERT_EQ_DBG(matrix.n_cols, m.n_cols);
206  matrix = m;
207 }
208 
210  ASSERT_EQ_DBG(matrix.n_rows, r.n_rows);
211  rhs = r;
212 }
213 
214 void LocalSystem::set_sparsity(const arma::umat & sp)
215 {
216  ASSERT_EQ_DBG(sparsity.n_rows, sp.n_rows);
217  ASSERT_EQ_DBG(sparsity.n_cols, sp.n_cols);
218 
219  sparsity.zeros();
220  for(unsigned int i=0; i < sp.n_rows; i++)
221  for(unsigned int j=0; j < sp.n_cols; j++)
222  if( sp(i,j) != 0 )
223  sparsity(i,j) = almost_zero;
224 // sparsity.print("sparsity");
225 }
void set_solution_col(uint loc_col, double solution)
Definition: local_system.cc:91
#define ASSERT_EQ_DBG(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:332
Mat< double, N, M > mat
Definition: armor.hh:214
DofVec elim_cols
unsigned int uint
void set_matrix(arma::mat matrix)
void set_sparsity(const arma::umat &sp)
Sets the sparsity pattern for the local system.
static constexpr double almost_zero
unsigned int n_elim_cols
void set_size(unsigned int nrows, unsigned int ncols)
Definition: local_system.cc:27
Mat< double, N, 1 > vec
Definition: armor.hh:211
DofVec col_dofs
Definition: local_system.hh:37
void set_solution(unsigned int 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:76
arma::Col< LongIdx > DofVec
Definition: local_system.hh:30
arma::vec rhs
local system RHS
arma::vec solution_rows
unsigned int n_elim_rows
arma::mat sparsity
sparsity pattern
void reset()
Resets the matrix, RHS, dofs to zero and clears solution settings.
Definition: local_system.cc:43
arma::vec diag_rows
LocalSystem()
Default constructor.
Definition: local_system.cc:8
arma::mat matrix
local system matrix
void set_solution_row(uint loc_row, double solution, double diag=0.0)
Definition: local_system.cc:84
Vector classes to support both Iterator, index and Id access and creating co-located vectors...
void set_rhs(arma::vec rhs)
#define ASSERT_DBG(expr)
arma::vec solution_cols
DofVec elim_rows
void add_value(unsigned int row, unsigned int col, double mat_val, double rhs_val)
Adds a single entry into the local system.
void eliminate_solution()
DofVec row_dofs
Definition: local_system.hh:37