16 #include <deal.II/base/function.h> 17 #include <deal.II/base/quadrature.h> 18 #include <deal.II/base/work_stream.h> 19 #include <deal.II/base/geometry_info.h> 20 #include <deal.II/base/quadrature.h> 21 #include <deal.II/dofs/dof_handler.h> 22 #include <deal.II/dofs/dof_accessor.h> 23 #include <deal.II/dofs/dof_tools.h> 24 #include <deal.II/fe/fe.h> 25 #include <deal.II/fe/fe_values.h> 26 #include <deal.II/fe/mapping_q1.h> 27 #include <deal.II/grid/tria_iterator.h> 28 #include <deal.II/hp/fe_values.h> 29 #include <deal.II/hp/mapping_collection.h> 30 #include <deal.II/numerics/matrix_tools.h> 31 #include <deal.II/lac/vector.h> 32 #include <deal.II/lac/block_vector.h> 33 #include <deal.II/lac/full_matrix.h> 34 #include <deal.II/lac/sparse_matrix.h> 35 #include <deal.II/lac/block_sparse_matrix.h> 37 #ifdef DEAL_II_WITH_PETSC 38 # include <deal.II/lac/petsc_parallel_sparse_matrix.h> 39 # include <deal.II/lac/petsc_sparse_matrix.h> 40 # include <deal.II/lac/petsc_parallel_vector.h> 41 # include <deal.II/lac/petsc_vector.h> 42 # include <deal.II/lac/petsc_parallel_block_sparse_matrix.h> 45 #ifdef DEAL_II_WITH_TRILINOS 46 # include <deal.II/lac/trilinos_sparse_matrix.h> 47 # include <deal.II/lac/trilinos_vector.h> 48 # include <deal.II/lac/trilinos_block_sparse_matrix.h> 49 # include <deal.II/lac/trilinos_block_vector.h> 60 DEAL_II_NAMESPACE_OPEN
69 template <
typename Iterator>
70 bool column_less_than(
const typename Iterator::value_type p,
71 const unsigned int column)
73 return (p.column() < column);
78 template <
typename number>
84 const bool eliminate_columns)
95 if (boundary_values.size() == 0)
107 number first_nonzero_diagonal_entry = 1;
108 for (
unsigned int i=0; i<n_dofs; ++i)
116 typename std::map<types::global_dof_index,number>::const_iterator dof = boundary_values.begin(),
117 endd = boundary_values.end();
118 for (; dof != endd; ++dof)
128 p = matrix.
begin(dof_number);
129 p != matrix.
end(dof_number); ++p)
130 if (p->column() != dof_number)
149 new_rhs = dof->second * matrix.
diag_element(dof_number);
150 right_hand_side(dof_number) = new_rhs;
154 matrix.
set (dof_number, dof_number,
155 first_nonzero_diagonal_entry);
156 new_rhs = dof->second * first_nonzero_diagonal_entry;
157 right_hand_side(dof_number) = new_rhs;
168 if (eliminate_columns)
173 const number diagonal_entry = matrix.
diag_element(dof_number);
184 q = matrix.
begin(dof_number)+1;
185 q != matrix.
end(dof_number); ++q)
193 const unsigned int column)
213 (p->column() == dof_number),
214 ExcMessage(
"This function is trying to access an element of the " 215 "matrix that doesn't seem to exist. Are you using a " 216 "nonsymmetric sparsity pattern? If so, you are not " 217 "allowed to set the eliminate_column argument of this " 218 "function, see the documentation."));
221 right_hand_side(row) -=
static_cast<number
>(p->value()) /
222 diagonal_entry * new_rhs;
230 solution(dof_number) = dof->second;
236 template <
typename number>
242 const bool eliminate_columns)
264 if (boundary_values.size() == 0)
276 number first_nonzero_diagonal_entry = 0;
277 for (
unsigned int diag_block=0; diag_block<blocks; ++diag_block)
279 for (
unsigned int i=0; i<matrix.
block(diag_block,diag_block).
n(); ++i)
282 first_nonzero_diagonal_entry
289 if (first_nonzero_diagonal_entry != 0)
294 if (first_nonzero_diagonal_entry == 0)
295 first_nonzero_diagonal_entry = 1;
298 typename std::map<types::global_dof_index,number>::const_iterator dof = boundary_values.begin(),
299 endd = boundary_values.end();
312 for (; dof != endd; ++dof)
321 const std::pair<unsigned int,types::global_dof_index>
331 for (
unsigned int block_col=0; block_col<blocks; ++block_col)
333 p = (block_col == block_index.first ?
334 matrix.
block(block_index.first,block_col).
begin(block_index.second) + 1 :
335 matrix.
block(block_index.first,block_col).
begin(block_index.second));
336 p != matrix.
block(block_index.first,block_col).
end(block_index.second);
354 if (matrix.
block(block_index.first, block_index.first)
356 new_rhs = dof->second *
357 matrix.
block(block_index.first, block_index.first)
361 matrix.
block(block_index.first, block_index.first)
363 = first_nonzero_diagonal_entry;
364 new_rhs = dof->second * first_nonzero_diagonal_entry;
366 right_hand_side.
block(block_index.first)(block_index.second)
379 if (eliminate_columns)
384 const number diagonal_entry
385 = matrix.
block(block_index.first,block_index.first)
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 q = (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);
440 const unsigned int column)
447 if (this_matrix.
begin(row)->column()
450 p = this_matrix.
begin(row);
453 this_matrix.
end(row),
459 this_matrix.
end(row),
473 Assert ((p->column() == block_index.second) &&
474 (p != this_matrix.
end(row)),
478 right_hand_side.
block(block_row)(row)
480 diagonal_entry * new_rhs;
489 solution.
block(block_index.first)(block_index.second) = dof->second;
497 template <
typename number>
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(),
508 Assert (local_dof_indices.size() == local_matrix.
n(),
511 Assert (local_dof_indices.size() == local_rhs.
size(),
517 if (boundary_values.size() == 0)
543 number average_diagonal = 0;
544 const unsigned int n_local_dofs = local_dof_indices.size();
545 for (
unsigned int i=0; i<n_local_dofs; ++i)
547 const typename std::map<types::global_dof_index, number>::const_iterator
548 boundary_value = boundary_values.find (local_dof_indices[i]);
549 if (boundary_value != boundary_values.end())
553 for (
unsigned int j=0; j<n_local_dofs; ++j)
555 local_matrix(i,j) = 0;
562 if (local_matrix(i,i) == 0.)
566 if (average_diagonal == 0.)
568 unsigned int nonzero_diagonals = 0;
569 for (
unsigned int k=0; k<n_local_dofs; ++k)
570 if (local_matrix(k,k) != 0.)
572 average_diagonal += std::fabs(local_matrix(k,k));
575 if (nonzero_diagonals != 0)
576 average_diagonal /= nonzero_diagonals;
578 average_diagonal = 0;
584 if (average_diagonal == 0.)
585 average_diagonal = 1.;
587 local_matrix(i,i) = average_diagonal;
590 local_matrix(i,i) = std::fabs(local_matrix(i,i));
594 local_rhs(i) = local_matrix(i,i) * boundary_value->second;
598 if (eliminate_columns ==
true)
600 for (
unsigned int row=0; row<n_local_dofs; ++row)
603 local_rhs(row) -= local_matrix(row,i) * boundary_value->second;
604 local_matrix(row,i) = 0;
615 #include "matrix_tools.inst" 618 DEAL_II_NAMESPACE_CLOSE
Iterator lower_bound(Iterator first, Iterator last, const T &val)
const BlockIndices & get_row_indices() const
void set(const size_type i, const size_type j, const number value)
unsigned int n_block_cols() const
SparsityPatternType & block(const size_type row, const size_type column)
const_iterator begin() const
static::ExceptionBase & ExcMessage(std::string arg1)
number diag_element(const size_type i) const
unsigned int global_dof_index
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const BlockIndices & get_block_indices() const
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
const Accessor< number, Constness > & value_type
const_iterator end() const
static::ExceptionBase & ExcNotQuadratic()
static::ExceptionBase & ExcBlocksDontMatch()
const BlockSparsityPattern & get_sparsity_pattern() const
const BlockIndices & get_column_indices() const
BlockType & block(const unsigned int row, const unsigned int column)
BlockType & block(const unsigned int i)
unsigned int n_block_rows() const
static::ExceptionBase & ExcInternalError()