#include <base/exceptions.h>
#include <base/thread_management.h>
+#include <dofs/function_map.h>
#include <map>
template <int dim> class DoFHandler;
template <int dim> class MGDoFHandler;
template <int dim> class FEValues;
-template <int dim> class FunctionMap;
*
* @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2001
*/
-//TODO: [WB] remove template arg and make functions static once FunctionMap is globalized
-template <int dim>
class MatrixCreator
{
public:
* See the general doc of this class
* for more information.
*/
+ template <int dim>
static void create_mass_matrix (const Mapping<dim> &mapping,
const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
* function, see above, with
* @p{mapping=MappingQ1<dim>()}.
*/
+ template <int dim>
static void create_mass_matrix (const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
* See the general doc of this
* class for more information.
*/
+ template <int dim>
static void create_mass_matrix (const Mapping<dim> &mapping,
const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
* function, see above, with
* @p{mapping=MappingQ1<dim>()}.
*/
+ template <int dim>
static void create_mass_matrix (const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
* See the general doc of this
* class for more information.
*/
+ template <int dim>
static
void create_boundary_mass_matrix (const Mapping<dim> &mapping,
const DoFHandler<dim> &dof,
const Quadrature<dim-1> &q,
SparseMatrix<double> &matrix,
- const typename FunctionMap<dim>::type &boundary_functions,
+ const typename FunctionMap<dim>::type &boundary_functions,
Vector<double> &rhs_vector,
std::vector<unsigned int>&dof_to_boundary_mapping,
const Function<dim> * const a = 0);
+ /**
+ * Same function, but for 1d.
+ */
+ static
+ void create_boundary_mass_matrix (const Mapping<1> &mapping,
+ const DoFHandler<1> &dof,
+ const Quadrature<0> &q,
+ SparseMatrix<double> &matrix,
+ const FunctionMap<1>::type &boundary_functions,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int>&dof_to_boundary_mapping,
+ const Function<1> * const a = 0);
+
+
/**
* Calls the
* @p{create_boundary_mass_matrix}
* function, see above, with
* @p{mapping=MappingQ1<dim>()}.
*/
+ template <int dim>
static
void create_boundary_mass_matrix (const DoFHandler<dim> &dof,
const Quadrature<dim-1> &q,
* See the general doc of this
* class for more information.
*/
+ template <int dim>
static void create_laplace_matrix (const Mapping<dim> &mapping,
const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
* function, see above, with
* @p{mapping=MappingQ1<dim>()}.
*/
+ template <int dim>
static void create_laplace_matrix (const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
* See the general doc of this
* class for more information.
*/
+ template <int dim>
static void create_level_laplace_matrix (unsigned int level,
const MGDoFHandler<dim>& dof,
const Quadrature<dim>& q,
* See the general doc of this
* class for more information.
*/
+ template <int dim>
static void create_laplace_matrix (const Mapping<dim> &mapping,
const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
* function, see above, with
* @p{mapping=MappingQ1<dim>()}.
*/
+ template <int dim>
static void create_laplace_matrix (const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
private:
/**
* Convenience abbreviation for
- * DoF handler cell iterators.
- */
- typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
-
- /**
- * Pair of iterators denoting a
- * half-open range.
+ * pairs of DoF handler cell
+ * iterators. This type works
+ * just like a
+ * @p{std::pair<iterator,iterator>}
+ * but is templatized on the
+ * space dimension.
*/
- typedef std::pair<active_cell_iterator,active_cell_iterator> IteratorRange;
+ template <int dim>
+ struct IteratorRange
+ {
+ /**
+ * Typedef for the iterator type.
+ */
+ typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+
+ /**
+ * Abbreviation for a pair of
+ * iterators.
+ */
+ typedef typename std::pair<active_cell_iterator,active_cell_iterator> iterator_pair;
+
+ /**
+ * Constructor. Initialize
+ * the two values by the
+ * given values.
+ */
+ IteratorRange (const active_cell_iterator &first,
+ const active_cell_iterator &second);
+
+ /**
+ * Constructor taking a pair
+ * of values for
+ * initialization.
+ */
+ IteratorRange (const iterator_pair &ip);
+
+ /**
+ * Pair of iterators denoting
+ * a half-open range.
+ */
+ active_cell_iterator first, second;
+ };
+
-
/**
* Version of the same function
* (without suffix @p{_1}) with
* used to synchronise access to
* the matrix.
*/
+ template <int dim>
static
void create_mass_matrix_1 (const Mapping<dim> &mapping,
const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
const Function<dim> * const a,
- const IteratorRange &range,
+ const IteratorRange<dim> range,
Threads::ThreadMutex &mutex);
/**
* used to synchronise access to
* the matrix.
*/
+ template <int dim>
static
void create_mass_matrix_2 (const Mapping<dim> &mapping,
const DoFHandler<dim> &dof,
const Function<dim> &rhs,
Vector<double> &rhs_vector,
const Function<dim> * const a,
- const IteratorRange &range,
+ const IteratorRange<dim> range,
Threads::ThreadMutex &mutex);
/**
* used to synchronise access to
* the matrix.
*/
+ template <int dim>
static
void create_laplace_matrix_1 (const Mapping<dim> &mapping,
const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
const Function<dim> * const a,
- const IteratorRange &range,
+ const IteratorRange<dim> range,
Threads::ThreadMutex &mutex);
/**
* used to synchronise access to
* the matrix.
*/
+ template <int dim>
static
void create_laplace_matrix_2 (const Mapping<dim> &mapping,
const DoFHandler<dim> &dof,
const Function<dim> &rhs,
Vector<double> &rhs_vector,
const Function<dim> * const a,
- const IteratorRange &range,
+ const IteratorRange<dim> range,
Threads::ThreadMutex &mutex);
/**
* used to synchronise access to
* the matrix.
*/
+ template <int dim>
static
void create_boundary_mass_matrix_1 (const Mapping<dim> &mapping,
const DoFHandler<dim> &dof,
Vector<double> &rhs_vector,
std::vector<unsigned int>&dof_to_boundary_mapping,
const Function<dim> * const a,
- const IteratorRange &range,
+ const IteratorRange<dim> range,
Threads::ThreadMutex &mutex);
};
-/* -------------- declaration of explicit specializations ------------- */
-
-
-template <> void MatrixCreator<1>::create_boundary_mass_matrix (
- const Mapping<1> &,
- const DoFHandler<1> &,
- const Quadrature<0> &,
- SparseMatrix<double> &,
- const FunctionMap<1>::type&,
- Vector<double> &,
- std::vector<unsigned int> &,
- const Function<1> * const);
-
-template <> void MatrixCreator<1>::create_boundary_mass_matrix_1 (
- const Mapping<1> &,
- const DoFHandler<1> &,
- const Quadrature<0> &,
- SparseMatrix<double> &,
- const FunctionMap<1>::type&,
- Vector<double> &,
- std::vector<unsigned int> &,
- const Function<1> * const ,
- const IteratorRange &,
- Threads::ThreadMutex &);
-
/**
* Provide a collection of functions operating on matrices. These include
*
* @author Wolfgang Bangerth, 1998, 2000
*/
-template <int dim>
-class MatrixTools : public MatrixCreator<dim>
+class MatrixTools : public MatrixCreator
{
public:
/**
--- /dev/null
+//---------------------------- matrices.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- matrices.cc ---------------------------
+
+
+#include <numerics/matrices.h>
+#include <lac/vector.h>
+#include <lac/block_vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
+
+#include <algorithm>
+
+
+
+template <typename number>
+void
+MatrixTools::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ SparseMatrix<number> &matrix,
+ Vector<number> &solution,
+ Vector<number> &right_hand_side,
+ const bool preserve_symmetry)
+{
+ Assert (matrix.n() == matrix.m(),
+ ExcDimensionsDontMatch(matrix.n(), matrix.m()));
+ Assert (matrix.n() == right_hand_side.size(),
+ ExcDimensionsDontMatch(matrix.n(), right_hand_side.size()));
+ Assert (matrix.n() == solution.size(),
+ ExcDimensionsDontMatch(matrix.n(), solution.size()));
+ // if no boundary values are to be applied
+ // simply return
+ if (boundary_values.size() == 0)
+ return;
+
+
+ std::map<unsigned int,double>::const_iterator dof = boundary_values.begin(),
+ endd = boundary_values.end();
+ const unsigned int n_dofs = matrix.m();
+ const SparsityPattern &sparsity = matrix.get_sparsity_pattern();
+ const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices();
+ const unsigned int *sparsity_colnums = sparsity.get_column_numbers();
+
+ // if a diagonal entry is zero
+ // later, then we use another
+ // number instead. take it to be
+ // the first nonzero diagonal
+ // element of the matrix, or 1 if
+ // there is no such thing
+ number first_nonzero_diagonal_entry = 1;
+ for (unsigned int i=0; i<n_dofs; ++i)
+ if (matrix.diag_element(i) != 0)
+ {
+ first_nonzero_diagonal_entry = matrix.diag_element(i);
+ break;
+ };
+
+
+ for (; dof != endd; ++dof)
+ {
+ Assert (dof->first < n_dofs, ExcInternalError());
+
+ const unsigned int dof_number = dof->first;
+ // for each boundary dof:
+
+ // set entries of this line
+ // to zero except for the diagonal
+ // entry. Note that the diagonal
+ // entry is always the first one
+ // for square matrices, i.e.
+ // we shall not set
+ // matrix.global_entry(
+ // sparsity_rowstart[dof.first])
+ const unsigned int last = sparsity_rowstart[dof_number+1];
+ for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
+ matrix.global_entry(j) = 0.;
+
+
+ // set right hand side to
+ // wanted value: if main diagonal
+ // entry nonzero, don't touch it
+ // and scale rhs accordingly. If
+ // zero, take the first main
+ // diagonal entry we can find, or
+ // one if no nonzero main diagonal
+ // element exists. Normally, however,
+ // the main diagonal entry should
+ // not be zero.
+ //
+ // store the new rhs entry to make
+ // the gauss step more efficient
+ number new_rhs;
+ if (matrix.diag_element(dof_number) != 0.0)
+ new_rhs = right_hand_side(dof_number)
+ = dof->second * matrix.diag_element(dof_number);
+ else
+ {
+ // use the SparseMatrix::
+ // to work around a bug in
+ // egcs
+ matrix.SparseMatrix<number>::set (dof_number, dof_number,
+ first_nonzero_diagonal_entry);
+ new_rhs = right_hand_side(dof_number)
+ = dof->second * first_nonzero_diagonal_entry;
+ };
+
+
+ // if the user wants to have
+ // the symmetry of the matrix
+ // preserved, and if the
+ // sparsity pattern is
+ // symmetric, then do a Gauss
+ // elimination step with the
+ // present row
+ if (preserve_symmetry)
+ {
+ // store the only nonzero entry
+ // of this line for the Gauss
+ // elimination step
+ const number diagonal_entry = matrix.diag_element(dof_number);
+
+ // we have to loop over all
+ // rows of the matrix which
+ // have a nonzero entry in
+ // the column which we work
+ // in presently. if the
+ // sparsity pattern is
+ // symmetric, then we can
+ // get the positions of
+ // these rows cheaply by
+ // looking at the nonzero
+ // column numbers of the
+ // present row. we need not
+ // look at the first entry,
+ // since that is the
+ // diagonal element and
+ // thus the present row
+ for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
+ {
+ const unsigned int row = sparsity_colnums[j];
+
+ // find the position of
+ // element
+ // (row,dof_number)
+ const unsigned int *
+ p = std::lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
+ &sparsity_colnums[sparsity_rowstart[row+1]],
+ dof_number);
+
+ // check whether this line has
+ // an entry in the regarding column
+ // (check for ==dof_number and
+ // != next_row, since if
+ // row==dof_number-1, *p is a
+ // past-the-end pointer but points
+ // to dof_number anyway...)
+ //
+ // there should be such an entry!
+ Assert ((*p == dof_number) &&
+ (p != &sparsity_colnums[sparsity_rowstart[row+1]]),
+ ExcInternalError());
+
+ const unsigned int global_entry
+ = (p - &sparsity_colnums[sparsity_rowstart[0]]);
+
+ // correct right hand side
+ right_hand_side(row) -= matrix.global_entry(global_entry) /
+ diagonal_entry * new_rhs;
+
+ // set matrix entry to zero
+ matrix.global_entry(global_entry) = 0.;
+ };
+ };
+
+ // preset solution vector
+ solution(dof_number) = dof->second;
+ };
+};
+
+
+
+
+
+void
+MatrixTools::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ BlockSparseMatrix<double> &matrix,
+ BlockVector<double> &solution,
+ BlockVector<double> &right_hand_side,
+ const bool preserve_symmetry)
+{
+ const unsigned int blocks = matrix.n_block_rows();
+
+ Assert (matrix.n() == matrix.m(),
+ ExcDimensionsDontMatch(matrix.n(), matrix.m()));
+ Assert (matrix.n() == right_hand_side.size(),
+ ExcDimensionsDontMatch(matrix.n(), right_hand_side.size()));
+ Assert (matrix.n() == solution.size(),
+ ExcDimensionsDontMatch(matrix.n(), solution.size()));
+ Assert (matrix.n_block_rows() == matrix.n_block_cols(),
+ ExcMatrixNotBlockSquare());
+ Assert (matrix.get_sparsity_pattern().get_row_indices() ==
+ matrix.get_sparsity_pattern().get_column_indices(),
+ ExcMatrixNotBlockSquare());
+ Assert (matrix.get_sparsity_pattern().get_column_indices() ==
+ solution.get_block_indices (),
+ ExcBlocksDontMatch ());
+ Assert (matrix.get_sparsity_pattern().get_row_indices() ==
+ right_hand_side.get_block_indices (),
+ ExcBlocksDontMatch ());
+
+
+ // if no boundary values are to be applied
+ // simply return
+ if (boundary_values.size() == 0)
+ return;
+
+
+ std::map<unsigned int,double>::const_iterator dof = boundary_values.begin(),
+ endd = boundary_values.end();
+ const unsigned int n_dofs = matrix.m();
+ const BlockSparsityPattern &
+ sparsity_pattern = matrix.get_sparsity_pattern();
+
+ // if a diagonal entry is zero
+ // later, then we use another
+ // number instead. take it to be
+ // the first nonzero diagonal
+ // element of the matrix, or 1 if
+ // there is no such thing
+ double first_nonzero_diagonal_entry = 0;
+ for (unsigned int diag_block=0; diag_block<blocks; ++diag_block)
+ {
+ for (unsigned int i=0; i<matrix.block(diag_block,diag_block).n(); ++i)
+ if (matrix.block(diag_block,diag_block).diag_element(i) != 0)
+ {
+ first_nonzero_diagonal_entry
+ = matrix.block(diag_block,diag_block).diag_element(i);
+ break;
+ };
+ // check whether we have found
+ // something in the present
+ // block
+ if (first_nonzero_diagonal_entry != 0)
+ break;
+ };
+ // nothing found on all diagonal
+ // blocks? if so, use 1.0 instead
+ if (first_nonzero_diagonal_entry == 0)
+ first_nonzero_diagonal_entry = 1;
+
+
+ // pointer to the mapping between
+ // global and block indices. since
+ // the row and column mappings are
+ // equal, store a pointer on only
+ // one of them
+ const BlockIndices &
+ index_mapping = sparsity_pattern.get_column_indices();
+
+ // now loop over all boundary dofs
+ for (; dof != endd; ++dof)
+ {
+ Assert (dof->first < n_dofs, ExcInternalError());
+
+ // get global index and index
+ // in the block in which this
+ // dof is located
+ const unsigned int dof_number = dof->first;
+ const std::pair<unsigned int,unsigned int>
+ block_index = index_mapping.global_to_local (dof_number);
+
+ // for each boundary dof:
+
+ // set entries of this line
+ // to zero except for the diagonal
+ // entry. Note that the diagonal
+ // entry is always the first one
+ // for square matrices, i.e.
+ // we shall not set
+ // matrix.global_entry(
+ // sparsity_rowstart[dof.first])
+ // of the diagonal block
+ for (unsigned int block_col=0; block_col<blocks; ++block_col)
+ {
+ const SparsityPattern &
+ local_sparsity = sparsity_pattern.block(block_index.first,
+ block_col);
+
+ // find first and last
+ // entry in the present row
+ // of the present
+ // block. exclude the main
+ // diagonal element, which
+ // is the diagonal element
+ // of a diagonal block,
+ // which must be a square
+ // matrix so the diagonal
+ // element is the first of
+ // this row.
+ const unsigned int
+ last = local_sparsity.get_rowstart_indices()[block_index.second+1],
+ first = (block_col == block_index.first ?
+ local_sparsity.get_rowstart_indices()[block_index.second]+1 :
+ local_sparsity.get_rowstart_indices()[block_index.second]);
+
+ for (unsigned int j=first; j<last; ++j)
+ matrix.block(block_index.first,block_col).global_entry(j) = 0.;
+ };
+
+
+ // set right hand side to
+ // wanted value: if main diagonal
+ // entry nonzero, don't touch it
+ // and scale rhs accordingly. If
+ // zero, take the first main
+ // diagonal entry we can find, or
+ // one if no nonzero main diagonal
+ // element exists. Normally, however,
+ // the main diagonal entry should
+ // not be zero.
+ //
+ // store the new rhs entry to make
+ // the gauss step more efficient
+ double new_rhs;
+ if (matrix.block(block_index.first, block_index.first)
+ .diag_element(block_index.second) != 0.0)
+ new_rhs = dof->second *
+ matrix.block(block_index.first, block_index.first)
+ .diag_element(block_index.second);
+ else
+ {
+ matrix.block(block_index.first, block_index.first)
+ .diag_element(block_index.second)
+ = first_nonzero_diagonal_entry;
+ new_rhs = dof->second * first_nonzero_diagonal_entry;
+ };
+ right_hand_side.block(block_index.first)(block_index.second)
+ = new_rhs;
+
+
+ // if the user wants to have
+ // the symmetry of the matrix
+ // preserved, and if the
+ // sparsity pattern is
+ // symmetric, then do a Gauss
+ // elimination step with the
+ // present row. this is a
+ // little more complicated for
+ // block matrices.
+ if (preserve_symmetry)
+ {
+ // store the only nonzero entry
+ // of this line for the Gauss
+ // elimination step
+ const double diagonal_entry
+ = matrix.block(block_index.first,block_index.first)
+ .diag_element(block_index.second);
+
+ // we have to loop over all
+ // rows of the matrix which
+ // have a nonzero entry in
+ // the column which we work
+ // in presently. if the
+ // sparsity pattern is
+ // symmetric, then we can
+ // get the positions of
+ // these rows cheaply by
+ // looking at the nonzero
+ // column numbers of the
+ // present row.
+ //
+ // note that if we check
+ // whether row @p{row} in
+ // block (r,c) is non-zero,
+ // then we have to check
+ // for the existence of
+ // column @p{row} in block
+ // (c,r), i.e. of the
+ // transpose block
+ for (unsigned int block_row=0; block_row<blocks; ++block_row)
+ {
+ // get pointers to the
+ // sparsity patterns of
+ // this block and of
+ // the transpose one
+ const SparsityPattern &this_sparsity
+ = sparsity_pattern.block (block_row, block_index.first);
+ const SparsityPattern &transpose_sparsity
+ = sparsity_pattern.block (block_index.first, block_row);
+
+ // traverse the row of
+ // the transpose block
+ // to find the
+ // interesting rows in
+ // the present block.
+ // don't use the
+ // diagonal element of
+ // the diagonal block
+ const unsigned int
+ first = (block_index.first == block_row ?
+ transpose_sparsity.get_rowstart_indices()[block_index.second]+1 :
+ transpose_sparsity.get_rowstart_indices()[block_index.second]),
+ last = transpose_sparsity.get_rowstart_indices()[block_index.second+1];
+
+ for (unsigned int j=first; j<last; ++j)
+ {
+ // get the number
+ // of the column in
+ // this row in
+ // which a nonzero
+ // entry is. this
+ // is also the row
+ // of the transpose
+ // block which has
+ // an entry in the
+ // interesting row
+ const unsigned int row = transpose_sparsity.get_column_numbers()[j];
+
+ // find the
+ // position of
+ // element
+ // (row,dof_number)
+ // in this block
+ // (not in the
+ // transpose
+ // one). note that
+ // we have to take
+ // care of special
+ // cases with
+ // square
+ // sub-matrices
+ const unsigned int *p = 0;
+ if (this_sparsity.n_rows() == this_sparsity.n_cols())
+ {
+ if (this_sparsity.get_column_numbers()
+ [this_sparsity.get_rowstart_indices()[row]]
+ ==
+ block_index.second)
+ p = &this_sparsity.get_column_numbers()
+ [this_sparsity.get_rowstart_indices()[row]];
+ else
+ p = std::lower_bound(&this_sparsity.get_column_numbers()
+ [this_sparsity.get_rowstart_indices()[row]+1],
+ &this_sparsity.get_column_numbers()
+ [this_sparsity.get_rowstart_indices()[row+1]],
+ block_index.second);
+ }
+ else
+ p = std::lower_bound(&this_sparsity.get_column_numbers()
+ [this_sparsity.get_rowstart_indices()[row]],
+ &this_sparsity.get_column_numbers()
+ [this_sparsity.get_rowstart_indices()[row+1]],
+ block_index.second);
+
+ // check whether this line has
+ // an entry in the regarding column
+ // (check for ==dof_number and
+ // != next_row, since if
+ // row==dof_number-1, *p is a
+ // past-the-end pointer but points
+ // to dof_number anyway...)
+ //
+ // there should be
+ // such an entry!
+ // note, however,
+ // that this
+ // assertion will
+ // fail sometimes
+ // if the sparsity
+ // pattern is not
+ // symmetric!
+ Assert ((*p == block_index.second) &&
+ (p != &this_sparsity.get_column_numbers()
+ [this_sparsity.get_rowstart_indices()[row+1]]),
+ ExcInternalError());
+
+ const unsigned int global_entry
+ = (p
+ -
+ &this_sparsity.get_column_numbers()
+ [this_sparsity.get_rowstart_indices()[0]]);
+
+ // correct right hand side
+ right_hand_side.block(block_row)(row)
+ -= matrix.block(block_row,block_index.first).global_entry(global_entry) /
+ diagonal_entry * new_rhs;
+
+ // set matrix entry to zero
+ matrix.block(block_row,block_index.first).global_entry(global_entry) = 0.;
+ };
+ };
+ };
+
+ // preset solution vector
+ solution.block(block_index.first)(block_index.second) = dof->second;
+ };
+};
+
+
+
+
+// explicit instantiations
+
+template
+void
+MatrixTools::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ SparseMatrix<double> &matrix,
+ Vector<double> &solution,
+ Vector<double> &right_hand_side,
+ const bool preserve_symmetry);
+template
+void
+MatrixTools::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ SparseMatrix<float> &matrix,
+ Vector<float> &solution,
+ Vector<float> &right_hand_side,
+ const bool preserve_symmetry);
+
#include <base/function.h>
+#include <base/thread_management.h>
#include <base/multithread_info.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
template <int dim>
-void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> * const coefficient)
+inline
+MatrixCreator::IteratorRange<dim>::
+IteratorRange (const active_cell_iterator &first,
+ const active_cell_iterator &second)
+ :
+ first (first),
+ second (second)
+{};
+
+
+
+template <int dim>
+inline
+MatrixCreator::IteratorRange<dim>::IteratorRange (const iterator_pair &ip)
+ :
+ first (ip.first),
+ second (ip.second)
+{};
+
+
+
+
+template <int dim>
+void MatrixCreator::create_mass_matrix (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> * const coefficient)
{
const unsigned int n_threads = multithread_info.n_default_threads;
Threads::ThreadManager thread_manager;
// define starting and end point
// for each thread
- std::vector<IteratorRange> thread_ranges
- = Threads::split_range<active_cell_iterator> (dof.begin_active(),
- dof.end(), n_threads);
+ typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ const std::vector<std::pair<active_cell_iterator,active_cell_iterator> >
+ thread_ranges = Threads::split_range<active_cell_iterator> (dof.begin_active(),
+ dof.end(), n_threads);
// mutex to synchronise access to
// the matrix
// then assemble in parallel
for (unsigned int thread=0; thread<n_threads; ++thread)
Threads::spawn (thread_manager,
- Threads::encapsulate(&MatrixCreator<dim>::create_mass_matrix_1)
+ Threads::encapsulate(&MatrixCreator::template create_mass_matrix_1<dim>)
.collect_args (mapping, dof, q, matrix, coefficient,
thread_ranges[thread], mutex));
thread_manager.wait ();
// TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
template <int dim>
-void MatrixCreator<dim>::create_mass_matrix_1 (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> * const coefficient,
- const IteratorRange &range,
- Threads::ThreadMutex &mutex)
+void MatrixCreator::create_mass_matrix_1 (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> * const coefficient,
+ const IteratorRange<dim> range,
+ Threads::ThreadMutex &mutex)
{
UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values);
if (coefficient != 0)
template <int dim>
-void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> * const coefficient)
+void MatrixCreator::create_mass_matrix (const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> * const coefficient)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
static const MappingQ1<dim> mapping;
template <int dim>
-void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> &rhs,
- Vector<double> &rhs_vector,
- const Function<dim> * const coefficient)
+void MatrixCreator::create_mass_matrix (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<dim> * const coefficient)
{
const unsigned int n_threads = multithread_info.n_default_threads;
Threads::ThreadManager thread_manager;
// define starting and end point
// for each thread
- std::vector<IteratorRange> thread_ranges
+ typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
= Threads::split_range<active_cell_iterator> (dof.begin_active(),
dof.end(), n_threads);
// then assemble in parallel
for (unsigned int thread=0; thread<n_threads; ++thread)
Threads::spawn (thread_manager,
- Threads::encapsulate(&MatrixCreator<dim>::
- create_mass_matrix_2)
+ Threads::encapsulate(&MatrixCreator::template
+ create_mass_matrix_2<dim>)
.collect_args (mapping, dof, q, matrix, rhs,
rhs_vector, coefficient,
thread_ranges[thread], mutex));
// TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
template <int dim>
void
-MatrixCreator<dim>::create_mass_matrix_2 (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> &rhs,
- Vector<double> &rhs_vector,
- const Function<dim> * const coefficient,
- const IteratorRange &range,
- Threads::ThreadMutex &mutex)
+MatrixCreator::create_mass_matrix_2 (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<dim> * const coefficient,
+ const IteratorRange<dim> range,
+ Threads::ThreadMutex &mutex)
{
UpdateFlags update_flags = UpdateFlags(update_values |
update_q_points |
template <int dim>
-void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> &rhs,
- Vector<double> &rhs_vector,
- const Function<dim> * const coefficient)
+void MatrixCreator::create_mass_matrix (const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<dim> * const coefficient)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
static const MappingQ1<dim> mapping;
#if deal_II_dimension == 1
-template <>
-void MatrixCreator<1>::create_boundary_mass_matrix (const Mapping<1> &,
- const DoFHandler<1> &,
- const Quadrature<0> &,
- SparseMatrix<double> &,
- const FunctionMap<1>::type&,
- Vector<double> &,
- std::vector<unsigned int> &,
- const Function<1> * const)
+void MatrixCreator::create_boundary_mass_matrix (const Mapping<1> &,
+ const DoFHandler<1> &,
+ const Quadrature<0> &,
+ SparseMatrix<double> &,
+ const FunctionMap<1>::type&,
+ Vector<double> &,
+ std::vector<unsigned int> &,
+ const Function<1> * const)
{
// what would that be in 1d? the
- // identity matrix on tehe boundary
+ // identity matrix on the boundary
// dofs?
Assert (false, ExcNotImplemented());
};
// TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
template <int dim>
void
-MatrixCreator<dim>::create_boundary_mass_matrix (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim-1> &q,
- SparseMatrix<double> &matrix,
- const typename FunctionMap<dim>::type &boundary_functions,
- Vector<double> &rhs_vector,
- std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<dim> * const a)
+MatrixCreator::create_boundary_mass_matrix (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim-1> &q,
+ SparseMatrix<double> &matrix,
+ const typename FunctionMap<dim>::type &boundary_functions,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<dim> * const a)
{
const unsigned int n_threads = multithread_info.n_default_threads;
Threads::ThreadManager thread_manager;
// define starting and end point
// for each thread
- std::vector<IteratorRange> thread_ranges
+ typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
= Threads::split_range<active_cell_iterator> (dof.begin_active(),
dof.end(), n_threads);
// then assemble in parallel
for (unsigned int thread=0; thread<n_threads; ++thread)
Threads::spawn (thread_manager,
- Threads::encapsulate(&MatrixCreator<dim>::
- create_boundary_mass_matrix_1)
+ Threads::encapsulate(&MatrixCreator::template
+ create_boundary_mass_matrix_1<dim>)
.collect_args (mapping, dof, q, matrix,
boundary_functions, rhs_vector,
dof_to_boundary_mapping, a,
};
-#if deal_II_dimension == 1
-
-template <>
-void
-MatrixCreator<1>::
-create_boundary_mass_matrix_1 (const Mapping<1> &,
- const DoFHandler<1> &,
- const Quadrature<0> &,
- SparseMatrix<double> &,
- const FunctionMap<1>::type&,
- Vector<double> &,
- std::vector<unsigned int> &,
- const Function<1> * const ,
- const IteratorRange &,
- Threads::ThreadMutex &)
-{
- Assert (false, ExcNotImplemented());
-};
-
-#endif
-
template <int dim>
void
-MatrixCreator<dim>::
+MatrixCreator::
create_boundary_mass_matrix_1 (const Mapping<dim> &mapping,
const DoFHandler<dim> &dof,
const Quadrature<dim-1> &q,
Vector<double> &rhs_vector,
std::vector<unsigned int> &dof_to_boundary_mapping,
const Function<dim> * const a,
- const IteratorRange &range,
+ const IteratorRange<dim> range,
Threads::ThreadMutex &mutex)
{
const FiniteElement<dim> &fe = dof.get_fe();
rhs_vector(dof_to_boundary_mapping[dofs[j]]) += cell_vector(j);
else
{
- // compare here for relative
- // smallness
+ // compare here for relative
+ // smallness
Assert (std::fabs(cell_vector(j)) <= 1e-10 * max_diag_entry,
ExcInternalError());
};
template <int dim>
-void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim> &dof,
- const Quadrature<dim-1> &q,
- SparseMatrix<double> &matrix,
- const typename FunctionMap<dim>::type &rhs,
- Vector<double> &rhs_vector,
- std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<dim> * const a)
+void MatrixCreator::create_boundary_mass_matrix (const DoFHandler<dim> &dof,
+ const Quadrature<dim-1> &q,
+ SparseMatrix<double> &matrix,
+ const typename FunctionMap<dim>::type &rhs,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<dim> * const a)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
static const MappingQ1<dim> mapping;
template <int dim>
-void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> * const coefficient)
+void MatrixCreator::create_laplace_matrix (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> * const coefficient)
{
const unsigned int n_threads = multithread_info.n_default_threads;
Threads::ThreadManager thread_manager;
// define starting and end point
// for each thread
- std::vector<IteratorRange> thread_ranges
+ typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
= Threads::split_range<active_cell_iterator> (dof.begin_active(),
dof.end(), n_threads);
// then assemble in parallel
for (unsigned int thread=0; thread<n_threads; ++thread)
Threads::spawn (thread_manager,
- Threads::encapsulate(&MatrixCreator<dim>::create_laplace_matrix_1)
+ Threads::encapsulate(&MatrixCreator::template
+ create_laplace_matrix_1<dim>)
.collect_args (mapping, dof, q, matrix, coefficient,
thread_ranges[thread], mutex));
thread_manager.wait ();
// TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
template <int dim>
-void MatrixCreator<dim>::create_laplace_matrix_1 (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> * const coefficient,
- const IteratorRange &range,
- Threads::ThreadMutex &mutex)
+void MatrixCreator::create_laplace_matrix_1 (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> * const coefficient,
+ const IteratorRange<dim> range,
+ Threads::ThreadMutex &mutex)
{
UpdateFlags update_flags = UpdateFlags(update_JxW_values |
update_gradients);
template <int dim>
-void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> * const coefficient)
+void MatrixCreator::create_laplace_matrix (const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> * const coefficient)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
static const MappingQ1<dim> mapping;
//TODO:[GK,RH] maybe recreate this function
/*
-template <int dim>
-void MatrixCreator<dim>::create_level_laplace_matrix (unsigned int level,
- const MGDoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<float> &matrix,
- const Function<dim> * const a)
-{
+ template <int dim>
+ void MatrixCreator::create_level_laplace_matrix (unsigned int level,
+ const MGDoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<float> &matrix,
+ const Function<dim> * const a)
+ {
Vector<double> dummy; // no entries, should give an error if accessed
UpdateFlags update_flags = UpdateFlags(update_gradients |
- update_JxW_values);
+ update_JxW_values);
if (a != 0)
- update_flags = UpdateFlags(update_flags | update_q_points);
+ update_flags = UpdateFlags(update_flags | update_q_points);
const Assembler<dim>::AssemblerData data (dof,
- true, false, // assemble matrix but not rhs
- matrix, dummy,
- q, update_flags);
+ true, false, // assemble matrix but not rhs
+ matrix, dummy,
+ q, update_flags);
TriaIterator<dim, Assembler<dim> >
- assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
- dof.get_tria().begin(level)->level(),
- dof.get_tria().begin(level)->index(),
- &data);
+ assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
+ dof.get_tria().begin(level)->level(),
+ dof.get_tria().begin(level)->index(),
+ &data);
LaplaceMatrix<dim> equation (0, a);
do
- {
- assembler->assemble (equation);
- ++assembler
- }
+ {
+ assembler->assemble (equation);
+ ++assembler
+ }
while ( (assembler.state()==valid) && (assembler->level() == level) );
-};
+ };
*/
template <int dim>
-void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> &rhs,
- Vector<double> &rhs_vector,
- const Function<dim> * const coefficient)
+void MatrixCreator::create_laplace_matrix (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<dim> * const coefficient)
{
const unsigned int n_threads = multithread_info.n_default_threads;
Threads::ThreadManager thread_manager;
// define starting and end point
// for each thread
- std::vector<IteratorRange> thread_ranges
+ typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
= Threads::split_range<active_cell_iterator> (dof.begin_active(),
dof.end(), n_threads);
// then assemble in parallel
for (unsigned int thread=0; thread<n_threads; ++thread)
Threads::spawn (thread_manager,
- Threads::encapsulate(&MatrixCreator<dim>::
- create_laplace_matrix_2)
+ Threads::encapsulate(&MatrixCreator::template
+ create_laplace_matrix_2<dim>)
.collect_args (mapping, dof, q, matrix, rhs,
rhs_vector, coefficient,
thread_ranges[thread], mutex));
// TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
template <int dim>
void
-MatrixCreator<dim>::create_laplace_matrix_2 (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> &rhs,
- Vector<double> &rhs_vector,
- const Function<dim> * const coefficient,
- const IteratorRange &range,
- Threads::ThreadMutex &mutex)
+MatrixCreator::create_laplace_matrix_2 (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<dim> * const coefficient,
+ const IteratorRange<dim> range,
+ Threads::ThreadMutex &mutex)
{
UpdateFlags update_flags = UpdateFlags(update_values |
update_gradients |
template <int dim>
-void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> &rhs,
- Vector<double> &rhs_vector,
- const Function<dim> * const coefficient)
+void MatrixCreator::create_laplace_matrix (const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<dim> * const coefficient)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
static const MappingQ1<dim> mapping;
-template <int dim>
-template <typename number>
-void
-MatrixTools<dim>::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- SparseMatrix<number> &matrix,
- Vector<number> &solution,
- Vector<number> &right_hand_side,
- const bool preserve_symmetry)
-{
- Assert (matrix.n() == matrix.m(),
- ExcDimensionsDontMatch(matrix.n(), matrix.m()));
- Assert (matrix.n() == right_hand_side.size(),
- ExcDimensionsDontMatch(matrix.n(), right_hand_side.size()));
- Assert (matrix.n() == solution.size(),
- ExcDimensionsDontMatch(matrix.n(), solution.size()));
- // if no boundary values are to be applied
- // simply return
- if (boundary_values.size() == 0)
- return;
-
-
- std::map<unsigned int,double>::const_iterator dof = boundary_values.begin(),
- endd = boundary_values.end();
- const unsigned int n_dofs = matrix.m();
- const SparsityPattern &sparsity = matrix.get_sparsity_pattern();
- const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices();
- const unsigned int *sparsity_colnums = sparsity.get_column_numbers();
-
- // if a diagonal entry is zero
- // later, then we use another
- // number instead. take it to be
- // the first nonzero diagonal
- // element of the matrix, or 1 if
- // there is no such thing
- number first_nonzero_diagonal_entry = 1;
- for (unsigned int i=0; i<n_dofs; ++i)
- if (matrix.diag_element(i) != 0)
- {
- first_nonzero_diagonal_entry = matrix.diag_element(i);
- break;
- };
-
-
- for (; dof != endd; ++dof)
- {
- Assert (dof->first < n_dofs, ExcInternalError());
-
- const unsigned int dof_number = dof->first;
- // for each boundary dof:
-
- // set entries of this line
- // to zero except for the diagonal
- // entry. Note that the diagonal
- // entry is always the first one
- // for square matrices, i.e.
- // we shall not set
- // matrix.global_entry(
- // sparsity_rowstart[dof.first])
- const unsigned int last = sparsity_rowstart[dof_number+1];
- for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
- matrix.global_entry(j) = 0.;
-
-
- // set right hand side to
- // wanted value: if main diagonal
- // entry nonzero, don't touch it
- // and scale rhs accordingly. If
- // zero, take the first main
- // diagonal entry we can find, or
- // one if no nonzero main diagonal
- // element exists. Normally, however,
- // the main diagonal entry should
- // not be zero.
- //
- // store the new rhs entry to make
- // the gauss step more efficient
- number new_rhs;
- if (matrix.diag_element(dof_number) != 0.0)
- new_rhs = right_hand_side(dof_number)
- = dof->second * matrix.diag_element(dof_number);
- else
- {
- // use the SparseMatrix::
- // to work around a bug in
- // egcs
- matrix.SparseMatrix<number>::set (dof_number, dof_number,
- first_nonzero_diagonal_entry);
- new_rhs = right_hand_side(dof_number)
- = dof->second * first_nonzero_diagonal_entry;
- };
-
-
- // if the user wants to have
- // the symmetry of the matrix
- // preserved, and if the
- // sparsity pattern is
- // symmetric, then do a Gauss
- // elimination step with the
- // present row
- if (preserve_symmetry)
- {
- // store the only nonzero entry
- // of this line for the Gauss
- // elimination step
- const number diagonal_entry = matrix.diag_element(dof_number);
-
- // we have to loop over all
- // rows of the matrix which
- // have a nonzero entry in
- // the column which we work
- // in presently. if the
- // sparsity pattern is
- // symmetric, then we can
- // get the positions of
- // these rows cheaply by
- // looking at the nonzero
- // column numbers of the
- // present row. we need not
- // look at the first entry,
- // since that is the
- // diagonal element and
- // thus the present row
- for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
- {
- const unsigned int row = sparsity_colnums[j];
-
- // find the position of
- // element
- // (row,dof_number)
- const unsigned int *
- p = std::lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
- &sparsity_colnums[sparsity_rowstart[row+1]],
- dof_number);
-
- // check whether this line has
- // an entry in the regarding column
- // (check for ==dof_number and
- // != next_row, since if
- // row==dof_number-1, *p is a
- // past-the-end pointer but points
- // to dof_number anyway...)
- //
- // there should be such an entry!
- Assert ((*p == dof_number) &&
- (p != &sparsity_colnums[sparsity_rowstart[row+1]]),
- ExcInternalError());
-
- const unsigned int global_entry
- = (p - &sparsity_colnums[sparsity_rowstart[0]]);
-
- // correct right hand side
- right_hand_side(row) -= matrix.global_entry(global_entry) /
- diagonal_entry * new_rhs;
-
- // set matrix entry to zero
- matrix.global_entry(global_entry) = 0.;
- };
- };
-
- // preset solution vector
- solution(dof_number) = dof->second;
- };
-};
-
-
-
-
-
-template <int dim>
-void
-MatrixTools<dim>::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- BlockSparseMatrix<double> &matrix,
- BlockVector<double> &solution,
- BlockVector<double> &right_hand_side,
- const bool preserve_symmetry)
-{
- const unsigned int blocks = matrix.n_block_rows();
-
- Assert (matrix.n() == matrix.m(),
- ExcDimensionsDontMatch(matrix.n(), matrix.m()));
- Assert (matrix.n() == right_hand_side.size(),
- ExcDimensionsDontMatch(matrix.n(), right_hand_side.size()));
- Assert (matrix.n() == solution.size(),
- ExcDimensionsDontMatch(matrix.n(), solution.size()));
- Assert (matrix.n_block_rows() == matrix.n_block_cols(),
- ExcMatrixNotBlockSquare());
- Assert (matrix.get_sparsity_pattern().get_row_indices() ==
- matrix.get_sparsity_pattern().get_column_indices(),
- ExcMatrixNotBlockSquare());
- Assert (matrix.get_sparsity_pattern().get_column_indices() ==
- solution.get_block_indices (),
- ExcBlocksDontMatch ());
- Assert (matrix.get_sparsity_pattern().get_row_indices() ==
- right_hand_side.get_block_indices (),
- ExcBlocksDontMatch ());
-
-
- // if no boundary values are to be applied
- // simply return
- if (boundary_values.size() == 0)
- return;
-
-
- std::map<unsigned int,double>::const_iterator dof = boundary_values.begin(),
- endd = boundary_values.end();
- const unsigned int n_dofs = matrix.m();
- const BlockSparsityPattern &
- sparsity_pattern = matrix.get_sparsity_pattern();
-
- // if a diagonal entry is zero
- // later, then we use another
- // number instead. take it to be
- // the first nonzero diagonal
- // element of the matrix, or 1 if
- // there is no such thing
- double first_nonzero_diagonal_entry = 0;
- for (unsigned int diag_block=0; diag_block<blocks; ++diag_block)
- {
- for (unsigned int i=0; i<matrix.block(diag_block,diag_block).n(); ++i)
- if (matrix.block(diag_block,diag_block).diag_element(i) != 0)
- {
- first_nonzero_diagonal_entry
- = matrix.block(diag_block,diag_block).diag_element(i);
- break;
- };
- // check whether we have found
- // something in the present
- // block
- if (first_nonzero_diagonal_entry != 0)
- break;
- };
- // nothing found on all diagonal
- // blocks? if so, use 1.0 instead
- if (first_nonzero_diagonal_entry == 0)
- first_nonzero_diagonal_entry = 1;
-
-
- // pointer to the mapping between
- // global and block indices. since
- // the row and column mappings are
- // equal, store a pointer on only
- // one of them
- const BlockIndices &
- index_mapping = sparsity_pattern.get_column_indices();
-
- // now loop over all boundary dofs
- for (; dof != endd; ++dof)
- {
- Assert (dof->first < n_dofs, ExcInternalError());
-
- // get global index and index
- // in the block in which this
- // dof is located
- const unsigned int dof_number = dof->first;
- const std::pair<unsigned int,unsigned int>
- block_index = index_mapping.global_to_local (dof_number);
-
- // for each boundary dof:
-
- // set entries of this line
- // to zero except for the diagonal
- // entry. Note that the diagonal
- // entry is always the first one
- // for square matrices, i.e.
- // we shall not set
- // matrix.global_entry(
- // sparsity_rowstart[dof.first])
- // of the diagonal block
- for (unsigned int block_col=0; block_col<blocks; ++block_col)
- {
- const SparsityPattern &
- local_sparsity = sparsity_pattern.block(block_index.first,
- block_col);
-
- // find first and last
- // entry in the present row
- // of the present
- // block. exclude the main
- // diagonal element, which
- // is the diagonal element
- // of a diagonal block,
- // which must be a square
- // matrix so the diagonal
- // element is the first of
- // this row.
- const unsigned int
- last = local_sparsity.get_rowstart_indices()[block_index.second+1],
- first = (block_col == block_index.first ?
- local_sparsity.get_rowstart_indices()[block_index.second]+1 :
- local_sparsity.get_rowstart_indices()[block_index.second]);
-
- for (unsigned int j=first; j<last; ++j)
- matrix.block(block_index.first,block_col).global_entry(j) = 0.;
- };
-
-
- // set right hand side to
- // wanted value: if main diagonal
- // entry nonzero, don't touch it
- // and scale rhs accordingly. If
- // zero, take the first main
- // diagonal entry we can find, or
- // one if no nonzero main diagonal
- // element exists. Normally, however,
- // the main diagonal entry should
- // not be zero.
- //
- // store the new rhs entry to make
- // the gauss step more efficient
- double new_rhs;
- if (matrix.block(block_index.first, block_index.first)
- .diag_element(block_index.second) != 0.0)
- new_rhs = dof->second *
- matrix.block(block_index.first, block_index.first)
- .diag_element(block_index.second);
- else
- {
- matrix.block(block_index.first, block_index.first)
- .diag_element(block_index.second)
- = first_nonzero_diagonal_entry;
- new_rhs = dof->second * first_nonzero_diagonal_entry;
- };
- right_hand_side.block(block_index.first)(block_index.second)
- = new_rhs;
-
-
- // if the user wants to have
- // the symmetry of the matrix
- // preserved, and if the
- // sparsity pattern is
- // symmetric, then do a Gauss
- // elimination step with the
- // present row. this is a
- // little more complicated for
- // block matrices.
- if (preserve_symmetry)
- {
- // store the only nonzero entry
- // of this line for the Gauss
- // elimination step
- const double diagonal_entry
- = matrix.block(block_index.first,block_index.first)
- .diag_element(block_index.second);
-
- // we have to loop over all
- // rows of the matrix which
- // have a nonzero entry in
- // the column which we work
- // in presently. if the
- // sparsity pattern is
- // symmetric, then we can
- // get the positions of
- // these rows cheaply by
- // looking at the nonzero
- // column numbers of the
- // present row.
- //
- // note that if we check
- // whether row @p{row} in
- // block (r,c) is non-zero,
- // then we have to check
- // for the existence of
- // column @p{row} in block
- // (c,r), i.e. of the
- // transpose block
- for (unsigned int block_row=0; block_row<blocks; ++block_row)
- {
- // get pointers to the
- // sparsity patterns of
- // this block and of
- // the transpose one
- const SparsityPattern &this_sparsity
- = sparsity_pattern.block (block_row, block_index.first);
- const SparsityPattern &transpose_sparsity
- = sparsity_pattern.block (block_index.first, block_row);
-
- // traverse the row of
- // the transpose block
- // to find the
- // interesting rows in
- // the present block.
- // don't use the
- // diagonal element of
- // the diagonal block
- const unsigned int
- first = (block_index.first == block_row ?
- transpose_sparsity.get_rowstart_indices()[block_index.second]+1 :
- transpose_sparsity.get_rowstart_indices()[block_index.second]),
- last = transpose_sparsity.get_rowstart_indices()[block_index.second+1];
-
- for (unsigned int j=first; j<last; ++j)
- {
- // get the number
- // of the column in
- // this row in
- // which a nonzero
- // entry is. this
- // is also the row
- // of the transpose
- // block which has
- // an entry in the
- // interesting row
- const unsigned int row = transpose_sparsity.get_column_numbers()[j];
-
- // find the
- // position of
- // element
- // (row,dof_number)
- // in this block
- // (not in the
- // transpose
- // one). note that
- // we have to take
- // care of special
- // cases with
- // square
- // sub-matrices
- const unsigned int *p = 0;
- if (this_sparsity.n_rows() == this_sparsity.n_cols())
- {
- if (this_sparsity.get_column_numbers()
- [this_sparsity.get_rowstart_indices()[row]]
- ==
- block_index.second)
- p = &this_sparsity.get_column_numbers()
- [this_sparsity.get_rowstart_indices()[row]];
- else
- p = std::lower_bound(&this_sparsity.get_column_numbers()
- [this_sparsity.get_rowstart_indices()[row]+1],
- &this_sparsity.get_column_numbers()
- [this_sparsity.get_rowstart_indices()[row+1]],
- block_index.second);
- }
- else
- p = std::lower_bound(&this_sparsity.get_column_numbers()
- [this_sparsity.get_rowstart_indices()[row]],
- &this_sparsity.get_column_numbers()
- [this_sparsity.get_rowstart_indices()[row+1]],
- block_index.second);
-
- // check whether this line has
- // an entry in the regarding column
- // (check for ==dof_number and
- // != next_row, since if
- // row==dof_number-1, *p is a
- // past-the-end pointer but points
- // to dof_number anyway...)
- //
- // there should be
- // such an entry!
- // note, however,
- // that this
- // assertion will
- // fail sometimes
- // if the sparsity
- // pattern is not
- // symmetric!
- Assert ((*p == block_index.second) &&
- (p != &this_sparsity.get_column_numbers()
- [this_sparsity.get_rowstart_indices()[row+1]]),
- ExcInternalError());
-
- const unsigned int global_entry
- = (p
- -
- &this_sparsity.get_column_numbers()
- [this_sparsity.get_rowstart_indices()[0]]);
-
- // correct right hand side
- right_hand_side.block(block_row)(row)
- -= matrix.block(block_row,block_index.first).global_entry(global_entry) /
- diagonal_entry * new_rhs;
-
- // set matrix entry to zero
- matrix.block(block_row,block_index.first).global_entry(global_entry) = 0.;
- };
- };
- };
-
- // preset solution vector
- solution.block(block_index.first)(block_index.second) = dof->second;
- };
-};
-
-
// explicit instantiations
-template class MatrixCreator<deal_II_dimension>;
-template class MatrixTools<deal_II_dimension>;
-
-
+template
+void MatrixCreator::create_mass_matrix (const Mapping<deal_II_dimension> &mapping,
+ const DoFHandler<deal_II_dimension> &dof,
+ const Quadrature<deal_II_dimension> &q,
+ SparseMatrix<double> &matrix,
+ const Function<deal_II_dimension> * const coefficient);
+template
+void MatrixCreator::create_mass_matrix (const DoFHandler<deal_II_dimension> &dof,
+ const Quadrature<deal_II_dimension> &q,
+ SparseMatrix<double> &matrix,
+ const Function<deal_II_dimension> * const coefficient);
+template
+void MatrixCreator::create_mass_matrix (const Mapping<deal_II_dimension> &mapping,
+ const DoFHandler<deal_II_dimension> &dof,
+ const Quadrature<deal_II_dimension> &q,
+ SparseMatrix<double> &matrix,
+ const Function<deal_II_dimension> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<deal_II_dimension> * const coefficient);
+template
+void MatrixCreator::create_mass_matrix (const DoFHandler<deal_II_dimension> &dof,
+ const Quadrature<deal_II_dimension> &q,
+ SparseMatrix<double> &matrix,
+ const Function<deal_II_dimension> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<deal_II_dimension> * const coefficient);
+
+#if deal_II_dimension != 1
template
void
-MatrixTools<deal_II_dimension>::
-apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- SparseMatrix<double> &matrix,
- Vector<double> &solution,
- Vector<double> &right_hand_side,
- const bool preserve_symmetry);
+MatrixCreator::create_boundary_mass_matrix (const Mapping<deal_II_dimension> &mapping,
+ const DoFHandler<deal_II_dimension> &dof,
+ const Quadrature<deal_II_dimension-1> &q,
+ SparseMatrix<double> &matrix,
+ const FunctionMap<deal_II_dimension>::type &boundary_functions,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<deal_II_dimension> * const a);
+#endif
template
-void
-MatrixTools<deal_II_dimension>::
-apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- SparseMatrix<float> &matrix,
- Vector<float> &solution,
- Vector<float> &right_hand_side,
- const bool preserve_symmetry);
+void MatrixCreator::create_boundary_mass_matrix (const DoFHandler<deal_II_dimension> &dof,
+ const Quadrature<deal_II_dimension-1> &q,
+ SparseMatrix<double> &matrix,
+ const FunctionMap<deal_II_dimension>::type &rhs,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<deal_II_dimension> * const a);
+template
+void MatrixCreator::create_laplace_matrix (const DoFHandler<deal_II_dimension> &dof,
+ const Quadrature<deal_II_dimension> &q,
+ SparseMatrix<double> &matrix,
+ const Function<deal_II_dimension> * const coefficient);
+
+template
+void MatrixCreator::create_laplace_matrix (const Mapping<deal_II_dimension> &mapping,
+ const DoFHandler<deal_II_dimension> &dof,
+ const Quadrature<deal_II_dimension> &q,
+ SparseMatrix<double> &matrix,
+ const Function<deal_II_dimension> * const coefficient);
+template
+void MatrixCreator::create_laplace_matrix (const Mapping<deal_II_dimension> &mapping,
+ const DoFHandler<deal_II_dimension> &dof,
+ const Quadrature<deal_II_dimension> &q,
+ SparseMatrix<double> &matrix,
+ const Function<deal_II_dimension> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<deal_II_dimension> * const coefficient);
+template
+void MatrixCreator::create_laplace_matrix (const DoFHandler<deal_II_dimension> &dof,
+ const Quadrature<deal_II_dimension> &q,
+ SparseMatrix<double> &matrix,
+ const Function<deal_II_dimension> &rhs,
+ Vector<double> &rhs_vector,
+ const Function<deal_II_dimension> * const coefficient);
projects / dealii-svn.git / commitdiff