76 const std::map<types::global_dof_index, number> &boundary_values,
80 const bool eliminate_columns)
86 Assert(matrix.n() == matrix.m(),
91 if (boundary_values.empty())
103 number first_nonzero_diagonal_entry = 1;
104 for (
unsigned int i = 0; i < n_dofs; ++i)
105 if (matrix.diag_element(i) != number())
107 first_nonzero_diagonal_entry = matrix.diag_element(i);
112 typename std::map<types::global_dof_index, number>::const_iterator
113 dof = boundary_values.begin(),
114 endd = boundary_values.end();
115 for (; dof != endd; ++dof)
125 matrix.begin(dof_number);
126 p != matrix.end(dof_number);
128 if (p->column() != dof_number)
145 if (matrix.diag_element(dof_number) != number())
147 new_rhs = dof->second * matrix.diag_element(dof_number);
148 right_hand_side(dof_number) = new_rhs;
152 matrix.set(dof_number, dof_number, first_nonzero_diagonal_entry);
153 new_rhs = dof->second * first_nonzero_diagonal_entry;
154 right_hand_side(dof_number) = new_rhs;
165 if (eliminate_columns)
170 const number diagonal_entry = matrix.diag_element(dof_number);
181 matrix.begin(dof_number) + 1;
182 q != matrix.end(dof_number);
192 const unsigned int column) =
210 Assert((p != matrix.end(row)) && (p->column() == dof_number),
212 "This function is trying to access an element of the "
213 "matrix that doesn't seem to exist. Are you using a "
214 "nonsymmetric sparsity pattern? If so, you are not "
215 "allowed to set the eliminate_column argument of this "
216 "function, see the documentation."));
219 right_hand_side(row) -=
220 static_cast<number
>(p->value()) / diagonal_entry * new_rhs;
228 solution(dof_number) = dof->second;
237 const std::map<types::global_dof_index, number> &boundary_values,
241 const bool eliminate_columns)
243 const unsigned int blocks = matrix.n_block_rows();
250 Assert(matrix.get_sparsity_pattern().get_row_indices() ==
251 matrix.get_sparsity_pattern().get_column_indices(),
253 Assert(matrix.get_sparsity_pattern().get_column_indices() ==
256 Assert(matrix.get_sparsity_pattern().get_row_indices() ==
262 if (boundary_values.empty())
274 number first_nonzero_diagonal_entry = 0;
275 for (
unsigned int diag_block = 0; diag_block < blocks; ++diag_block)
277 for (
unsigned int i = 0; i < matrix.block(diag_block, diag_block).n();
279 if (matrix.block(diag_block, diag_block).diag_element(i) != number{})
281 first_nonzero_diagonal_entry =
282 matrix.block(diag_block, diag_block).diag_element(i);
288 if (first_nonzero_diagonal_entry != number{})
293 if (first_nonzero_diagonal_entry == number{})
294 first_nonzero_diagonal_entry = 1;
297 typename std::map<types::global_dof_index, number>::const_iterator
298 dof = boundary_values.begin(),
299 endd = boundary_values.end();
301 matrix.get_sparsity_pattern();
311 for (; dof != endd; ++dof)
319 const std::pair<unsigned int, types::global_dof_index> block_index =
329 for (
unsigned int block_col = 0; block_col < blocks; ++block_col)
331 (block_col == block_index.first ?
332 matrix.block(block_index.first, block_col)
333 .begin(block_index.second) +
335 matrix.block(block_index.first, block_col)
336 .begin(block_index.second));
337 p != matrix.block(block_index.first, block_col)
338 .end(block_index.second);
356 if (matrix.block(block_index.first, block_index.first)
357 .diag_element(block_index.second) != number{})
359 dof->second * matrix.block(block_index.first, block_index.first)
360 .diag_element(block_index.second);
363 matrix.block(block_index.first, block_index.first)
364 .diag_element(block_index.second) = first_nonzero_diagonal_entry;
365 new_rhs = dof->second * first_nonzero_diagonal_entry;
367 right_hand_side.
block(block_index.first)(block_index.second) = new_rhs;
379 if (eliminate_columns)
384 const number diagonal_entry =
385 matrix.block(block_index.first, block_index.first)
386 .diag_element(block_index.second);
409 for (
unsigned int block_row = 0; block_row < blocks; ++block_row)
414 sparsity_pattern.
block(block_row, block_index.first);
417 matrix.block(block_row, block_index.first);
419 matrix.block(block_index.first, block_row);
425 (block_index.first == block_row ?
426 transpose_matrix.
begin(block_index.second) + 1 :
427 transpose_matrix.
begin(block_index.second));
428 q != transpose_matrix.
end(block_index.second);
441 const unsigned int column) =
450 if (this_matrix.
begin(row)->column() ==
452 p = this_matrix.
begin(row);
455 this_matrix.
end(row),
461 this_matrix.
end(row),
475 Assert((p->column() == block_index.second) &&
476 (p != this_matrix.
end(row)),
480 right_hand_side.
block(block_row)(row) -=
481 number(p->value()) / diagonal_entry * new_rhs;
490 solution.
block(block_index.first)(block_index.second) = dof->second;
499 const std::map<types::global_dof_index, number> &boundary_values,
500 const std::vector<types::global_dof_index> &local_dof_indices,
503 const bool eliminate_columns)
505 Assert(local_dof_indices.size() == local_matrix.
m(),
507 Assert(local_dof_indices.size() == local_matrix.
n(),
509 Assert(local_dof_indices.size() == local_rhs.
size(),
514 if (boundary_values.empty())
540 number average_diagonal = 0;
541 const unsigned int n_local_dofs = local_dof_indices.size();
542 for (
unsigned int i = 0; i < n_local_dofs; ++i)
544 const typename std::map<types::global_dof_index, number>::const_iterator
545 boundary_value = boundary_values.find(local_dof_indices[i]);
546 if (boundary_value != boundary_values.end())
550 for (
unsigned int j = 0; j < n_local_dofs; ++j)
552 local_matrix(i, j) = 0;
559 if (local_matrix(i, i) == number{})
563 if (average_diagonal == number{})
565 unsigned int nonzero_diagonals = 0;
566 for (
unsigned int k = 0; k < n_local_dofs; ++k)
567 if (local_matrix(k, k) != number{})
569 average_diagonal +=
std::abs(local_matrix(k, k));
572 if (nonzero_diagonals != 0)
573 average_diagonal /= nonzero_diagonals;
575 average_diagonal = 0;
581 if (average_diagonal == number{})
582 average_diagonal = 1.;
584 local_matrix(i, i) = average_diagonal;
587 local_matrix(i, i) =
std::abs(local_matrix(i, i));
591 local_rhs(i) = local_matrix(i, i) * boundary_value->second;
595 if (eliminate_columns ==
true)
597 for (
unsigned int row = 0; row < n_local_dofs; ++row)
601 local_matrix(row, i) * boundary_value->second;
602 local_matrix(row, i) = 0;