From: Guido Kanschat Date: Sun, 21 Mar 2010 19:22:22 +0000 (+0000) Subject: working on constraints for matrix blocks X-Git-Tag: v8.0.0~6289 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ae1409b2a1dfac6d044f20d362b7662c3dcd69e2;p=dealii.git working on constraints for matrix blocks git-svn-id: https://svn.dealii.org/trunk@20871 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/multigrid/mg_level_object.h b/deal.II/base/include/base/mg_level_object.h similarity index 97% rename from deal.II/deal.II/include/multigrid/mg_level_object.h rename to deal.II/base/include/base/mg_level_object.h index 1b97b476a8..7b75ed4154 100644 --- a/deal.II/deal.II/include/multigrid/mg_level_object.h +++ b/deal.II/base/include/base/mg_level_object.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2003, 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 2003, 2004, 2005, 2006, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -21,9 +21,6 @@ DEAL_II_NAMESPACE_OPEN -/*!@addtogroup mg */ -/*@{*/ - /** * An array with an object for each level. The purpose of this class * is mostly to store objects and allow access by level number, even @@ -34,7 +31,9 @@ DEAL_II_NAMESPACE_OPEN * In most cases, the objects which are stored on each levels, are * either matrices or vectors. * - * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2005 + * @ingroup mg + * @ingroup data + * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2005, 2010 */ template class MGLevelObject : public Subscriptor @@ -116,7 +115,6 @@ class MGLevelObject : public Subscriptor std::vector > objects; }; -/*@}*/ /* ------------------------------------------------------------------- */ diff --git a/deal.II/deal.II/include/multigrid/mg_block_smoother.h b/deal.II/deal.II/include/multigrid/mg_block_smoother.h index 3065fc7e90..b51afecf49 100644 --- a/deal.II/deal.II/include/multigrid/mg_block_smoother.h +++ b/deal.II/deal.II/include/multigrid/mg_block_smoother.h @@ -20,7 +20,7 @@ #include #include #include -#include +#include #include DEAL_II_NAMESPACE_OPEN diff --git a/deal.II/deal.II/include/multigrid/mg_matrix.h b/deal.II/deal.II/include/multigrid/mg_matrix.h index d66cf2c3de..274b1eb3ea 100644 --- a/deal.II/deal.II/include/multigrid/mg_matrix.h +++ b/deal.II/deal.II/include/multigrid/mg_matrix.h @@ -16,7 +16,7 @@ #include #include #include -#include +#include DEAL_II_NAMESPACE_OPEN diff --git a/deal.II/deal.II/include/multigrid/mg_smoother.h b/deal.II/deal.II/include/multigrid/mg_smoother.h index d298470a95..e28f7cd486 100644 --- a/deal.II/deal.II/include/multigrid/mg_smoother.h +++ b/deal.II/deal.II/include/multigrid/mg_smoother.h @@ -19,7 +19,7 @@ #include #include #include -#include +#include #include DEAL_II_NAMESPACE_OPEN diff --git a/deal.II/deal.II/include/multigrid/mg_transfer.h b/deal.II/deal.II/include/multigrid/mg_transfer.h index 8c63dc237a..7db0779cc5 100644 --- a/deal.II/deal.II/include/multigrid/mg_transfer.h +++ b/deal.II/deal.II/include/multigrid/mg_transfer.h @@ -27,7 +27,7 @@ #include #include -#include +#include diff --git a/deal.II/deal.II/include/multigrid/mg_transfer_block.h b/deal.II/deal.II/include/multigrid/mg_transfer_block.h index bc895bc3ec..a9f7451ef2 100644 --- a/deal.II/deal.II/include/multigrid/mg_transfer_block.h +++ b/deal.II/deal.II/include/multigrid/mg_transfer_block.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -27,7 +27,7 @@ #include #include -#include +#include diff --git a/deal.II/deal.II/include/multigrid/mg_transfer_component.h b/deal.II/deal.II/include/multigrid/mg_transfer_component.h index 2d1e1de352..3a209635f4 100644 --- a/deal.II/deal.II/include/multigrid/mg_transfer_component.h +++ b/deal.II/deal.II/include/multigrid/mg_transfer_component.h @@ -27,7 +27,7 @@ #include #include -#include +#include diff --git a/deal.II/deal.II/include/multigrid/multigrid.h b/deal.II/deal.II/include/multigrid/multigrid.h index 8442d467ef..fbcccdce6b 100644 --- a/deal.II/deal.II/include/multigrid/multigrid.h +++ b/deal.II/deal.II/include/multigrid/multigrid.h @@ -21,7 +21,7 @@ #include #include #include -#include +#include #include #include diff --git a/deal.II/deal.II/include/numerics/mesh_worker.h b/deal.II/deal.II/include/numerics/mesh_worker.h index 3602b4331f..6f9eda3ae4 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker.h +++ b/deal.II/deal.II/include/numerics/mesh_worker.h @@ -18,7 +18,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h index 01371224c7..e1af235b0b 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h @@ -17,7 +17,7 @@ #include #include #include -#include +#include DEAL_II_NAMESPACE_OPEN diff --git a/deal.II/deal.II/include/numerics/mesh_worker_output.h b/deal.II/deal.II/include/numerics/mesh_worker_output.h index 0e557c67af..1ae30f1ec8 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_output.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_output.h @@ -17,7 +17,7 @@ #include #include #include -#include +#include DEAL_II_NAMESPACE_OPEN diff --git a/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h index 80312e8f88..41f709fd65 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h @@ -16,7 +16,7 @@ #include #include #include -#include +#include DEAL_II_NAMESPACE_OPEN diff --git a/deal.II/deal.II/source/multigrid/mg_base.cc b/deal.II/deal.II/source/multigrid/mg_base.cc index 90476e759b..e7088cdfd7 100644 --- a/deal.II/deal.II/source/multigrid/mg_base.cc +++ b/deal.II/deal.II/source/multigrid/mg_base.cc @@ -17,7 +17,7 @@ #include #include #include -#include +#include #include diff --git a/deal.II/deal.II/source/multigrid/mg_tools.all_dimensions.cc b/deal.II/deal.II/source/multigrid/mg_tools.all_dimensions.cc index e85e42f8fd..8a7fbba6f1 100644 --- a/deal.II/deal.II/source/multigrid/mg_tools.all_dimensions.cc +++ b/deal.II/deal.II/source/multigrid/mg_tools.all_dimensions.cc @@ -18,7 +18,7 @@ #include #include #include -#include +#include #include #include diff --git a/deal.II/deal.II/source/multigrid/mg_tools.cc b/deal.II/deal.II/source/multigrid/mg_tools.cc index 92eaf98f41..035e80a42c 100644 --- a/deal.II/deal.II/source/multigrid/mg_tools.cc +++ b/deal.II/deal.II/source/multigrid/mg_tools.cc @@ -26,7 +26,7 @@ #include #include #include -#include +#include #include #include diff --git a/deal.II/examples/step-12/step-12.cc b/deal.II/examples/step-12/step-12.cc index 8eff89bddf..64322d2285 100644 --- a/deal.II/examples/step-12/step-12.cc +++ b/deal.II/examples/step-12/step-12.cc @@ -456,14 +456,14 @@ void Step12::integrate_cell_term (DoFInfo& dinfo, CellInfo& info) } // Now the same for the boundary terms. Note - // that now we use FEFaceValuesBase, the base + // that now we use FEValuesBase, the base // class for both FEFaceValues and // FESubfaceValues, in order to get access to // normal vectors. template void Step12::integrate_boundary_term (DoFInfo& dinfo, FaceInfo& info) { - const FEFaceValuesBase& fe_v = info.fe_values(); + const FEValuesBase& fe_v = info.fe_values(); FullMatrix& local_matrix = dinfo.matrix(0).matrix; Vector& local_vector = dinfo.vector(0).block(0); @@ -512,14 +512,14 @@ void Step12::integrate_face_term (DoFInfo& dinfo1, DoFInfo& dinfo2, { // For quadrature points, weights, // etc., we use the - // FEFaceValuesBase object of the + // FEValuesBase object of the // first argument. - const FEFaceValuesBase& fe_v = info1.fe_values(); + const FEValuesBase& fe_v = info1.fe_values(); // For additional shape functions, // we have to ask the neighbors - // FEFaceValuesBase. - const FEFaceValuesBase& fe_v_neighbor = info2.fe_values(); + // FEValuesBase. + const FEValuesBase& fe_v_neighbor = info2.fe_values(); // Then we get references to the // four local matrices. The letters diff --git a/deal.II/lac/include/lac/constraint_matrix.h b/deal.II/lac/include/lac/constraint_matrix.h index c4c21c4727..79ffcfe2b3 100644 --- a/deal.II/lac/include/lac/constraint_matrix.h +++ b/deal.II/lac/include/lac/constraint_matrix.h @@ -1825,6 +1825,18 @@ class ConstraintMatrix : public Subscriptor const Table<2,bool> &dof_mask, internal::bool2type) const; + /** + * Internal helper function for + * distribute_local_to_global + * function. + * + * Creates a list of affected + * rows for sistribution. + */ + void + make_sorted_row_list (const std::vector &local_dof_indices, + internals::GlobalRowsFromLocal &global_rows) const; + /** * Internal helper function for * distribute_local_to_global function. diff --git a/deal.II/lac/include/lac/constraint_matrix.templates.h b/deal.II/lac/include/lac/constraint_matrix.templates.h index 5f20f9465f..d9a2516ad3 100644 --- a/deal.II/lac/include/lac/constraint_matrix.templates.h +++ b/deal.II/lac/include/lac/constraint_matrix.templates.h @@ -1722,6 +1722,85 @@ namespace internals } // end of namespace internals +void +ConstraintMatrix:: +make_sorted_row_list (const std::vector &local_dof_indices, + internals::GlobalRowsFromLocal &global_rows) const +{ + const unsigned int n_local_dofs = local_dof_indices.size(); + // when distributing the local data to + // the global matrix, we can quite + // cheaply sort the indices (obviously, + // this introduces the need for + // allocating some memory on the way, but + // we need to do this only for rows, + // whereas the distribution process + // itself goes over rows and + // columns). This has the advantage that + // when writing into the global matrix, + // we can make use of the sortedness. + + // so the first step is to create a + // sorted list of all row values that are + // possible. these values are either the + // rows from unconstrained dofs, or some + // indices introduced by dofs constrained + // to a combination of some other + // dofs. regarding the data type, choose + // an STL vector of a pair of unsigned + // ints (for global columns) and internal + // data (containing local columns + + // possible jumps from + // constraints). Choosing an STL map or + // anything else M.K. knows of would be + // much more expensive here! + + // cache whether we have to resolve any + // indirect rows generated from resolving + // constrained dofs. + unsigned int added_rows = 0; + + // first add the indices in an unsorted + // way and only keep track of the + // constraints that appear. They are + // resolved in a second step. + for (unsigned int i = 0; iM-1FA-1 with M the + * pressure mass matrix, A the pressure Laplacian and F + * the advection diffusion operator applied to the pressure + * space. Since only a single matrix is needed for the other blocks, + * using BlockSparseMatrix or similar would be a waste of memory. + * + * While the add() functions make a MatrixBlock appear like a block + * matrix for assembling, the functions vmult(), Tvmult(), + * vmult_add(), and Tvmult_add() make it behave like a MATRIX, when it + * comes to applying it to a vector. This behavior allows us to store + * MatrixBlock objects in vectors, for instance in MGLevelObject + * without extracting the #matrix first. * * MatrixBlock comes handy when using BlockMatrixArray. Once the * MatrixBlock has been properly initalized and filled, it can be used - * as: + * in the simplest case as: * @code - * std::vector > > blocks; + * MatrixBlockVector > > blocks; * * ... * * BlockMatrixArray matrix (n_blocks, n_blocks); * * for (unsigned int i=0;i @@ -79,33 +101,34 @@ class MatrixBlock /** * Add value to the - * element (i,j). Throws + * element (i,j). Throws * an error if the entry does not - * exist or if value is - * not a finite number. Still, it - * is allowed to store zero - * values in non-existent fields. + * exist or if it is in a + * different block. */ void add (const unsigned int i, const unsigned int j, const typename MATRIX::value_type value); /** - * Add all elements given in a - * FullMatrix into sparse + * Add all elements in a + * FullMatrix into sparse * matrix locations given by - * indices. In other words, - * this function adds the elements in - * full_matrix to the - * respective entries in calling - * matrix, using the local-to-global - * indexing specified by - * indices for both the rows - * and the columns of the - * matrix. This function assumes a - * quadratic sparse matrix and a - * quadratic full_matrix, the usual - * situation in FE calculations. + * indices. This function + * assumes a quadratic sparse + * matrix and a quadratic + * full_matrix. The global + * locations are translated + * into locations in this block + * and ExcBlockIndexMismatch is + * thrown, if the global index + * does not point into the + * block refered to by #row and + * #column. + * + * @todo + * elide_zero_values is + * currently ignored. * * The optional parameter * elide_zero_values can be @@ -123,11 +146,33 @@ class MatrixBlock const bool elide_zero_values = true); /** - * Same function as before, but now - * including the possibility to use - * rectangular full_matrices and - * different local-to-global indexing - * on rows and columns, respectively. + * Add all elements in a + * FullMatrix into global + * locations given by + * row_indices and + * col_indices, + * respectively. The global + * locations are translated + * into locations in this block + * and ExcBlockIndexMismatch is + * thrown, if the global index + * does not point into the + * block refered to by #row and + * #column. + * + * @todo + * elide_zero_values is + * currently ignored. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be added anyway or + * these should be filtered away and + * only non-zero data is added. The + * default value is true, + * i.e., zero values won't be added + * into the matrix. */ template void add (const std::vector& row_indices, @@ -137,10 +182,27 @@ class MatrixBlock /** * Set several elements in the - * specified row of the matrix with - * column indices as given by - * col_indices to the - * respective value. + * specified row of the matrix + * with column indices as given + * by col_indices to + * the respective value. This + * is the function doing thye + * actual work for the ones + * adding full matrices. The + * global locations + * row_index and + * col_indices are + * translated into locations in + * this block and + * ExcBlockIndexMismatch is + * thrown, if the global index + * does not point into the + * block refered to by #row and + * #column. + * + * @todo + * elide_zero_values is + * currently ignored. * * The optional parameter * elide_zero_values can be @@ -153,7 +215,7 @@ class MatrixBlock * into the matrix. */ template - void add (const unsigned int row, + void add (const unsigned int row_index, const std::vector& col_indices, const std::vector& values, const bool elide_zero_values = true); @@ -183,6 +245,50 @@ class MatrixBlock const bool elide_zero_values = true, const bool col_indices_are_sorted = false); + /** + * Matrix-vector-multiplication, + * forwarding to the same + * function in MATRIX. No index + * computations are done, thus, + * the vectors need to have sizes + * matching #matrix. + */ + template + void vmult (VECTOR& w, const VECTOR& v) const; + + /** + * Matrix-vector-multiplication, + * forwarding to the same + * function in MATRIX. No index + * computations are done, thus, + * the vectors need to have sizes + * matching #matrix. + */ + template + void vmult_add (VECTOR& w, const VECTOR& v) const; + + /** + * Matrix-vector-multiplication, + * forwarding to the same + * function in MATRIX. No index + * computations are done, thus, + * the vectors need to have sizes + * matching #matrix. + */ + template + void Tvmult (VECTOR& w, const VECTOR& v) const; + + /** + * Matrix-vector-multiplication, + * forwarding to the same + * function in MATRIX. No index + * computations are done, thus, + * the vectors need to have sizes + * matching #matrix. + */ + template + void Tvmult_add (VECTOR& w, const VECTOR& v) const; + /** * The block number computed from * an index by using @@ -229,6 +335,8 @@ class MatrixBlock * pointers, in order to allow for copying and rearranging. Each * matrix block can be identified by name. * + * @relates MatrixBlock + * @ingroup vector_valued * @author Baerbel Janssen, Guido Kanschat, 2010 */ template @@ -264,6 +372,48 @@ class MatrixBlockVector : public NamedData }; +/** + * A vector of MGLevelObject, which is implemented using shared + * pointers, in order to allow for copying and rearranging. Each + * matrix block can be identified by name. + * + * @relates MatrixBlock + * @ingroup vector_valued + * @author Baerbel Janssen, Guido Kanschat, 2010 + */ +template +class MGMatrixBlockVector : public NamedData > > +{ + public: + /** + * Add a new matrix block at the + * position + * in the block system. + */ + void add(unsigned int row, unsigned int column, + const std::string& name, + const BlockIndices* block_indices); + /** + * Access a constant reference to + * the block at position i. + */ + const MatrixBlock& block(unsigned int i) const; + /** + * Access a reference to + * the block at position i. + */ + MatrixBlock& block(unsigned int i); + /** + * Access the matrix at position + * i for read and write + * access. + */ + MATRIX& matrix(unsigned int i); + + +}; + + //----------------------------------------------------------------------// template @@ -336,6 +486,43 @@ MatrixBlock::add (const std::vector& row_indices, } +template +template +inline +void +MatrixBlock::add (const unsigned int b_row, + const unsigned int n_cols, + const unsigned int *col_indices, + const number *values, + const bool, + const bool) +{ + Assert(block_indices != 0, ExcNotInitialized()); + const std::pair bi + = block_indices->global_to_local(b_row); + + // In debug mode, we check whether + // all indices are in the correct + // block. + + // Actually, for the time being, we + // leave it at this. While it may + // not be the most efficient way, + // it is at least thread safe. +//#ifdef DEBUG + Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row)); + + for (unsigned int j=0;j bj + = block_indices->global_to_local(col_indices[j]); + Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column)); + + matrix.add(bi.second, bj.second, values[j]); + } +//#endif +} + template template @@ -374,45 +561,46 @@ MatrixBlock::add (const unsigned int row, template -template +template inline void -MatrixBlock::add (const unsigned int b_row, - const unsigned int n_cols, - const unsigned int *col_indices, - const number *values, - const bool, - const bool) +MatrixBlock::vmult (VECTOR& w, const VECTOR& v) const { - Assert(block_indices != 0, ExcNotInitialized()); - const std::pair bi - = block_indices->global_to_local(b_row); + matrix.vmult(w,v); +} - // In debug mode, we check whether - // all indices are in the correct - // block. - // Actually, for the time being, we - // leave it at this. While it may - // not be the most efficient way, - // it is at least thread safe. -//#ifdef DEBUG - Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row)); +template +template +inline +void +MatrixBlock::vmult_add (VECTOR& w, const VECTOR& v) const +{ + matrix.vmult_add(w,v); +} - for (unsigned int j=0;j bj - = block_indices->global_to_local(col_indices[j]); - Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column)); - - matrix.add(bi.second, bj.second, values[j]); - } - -//#endif +template +template +inline +void +MatrixBlock::Tvmult (VECTOR& w, const VECTOR& v) const +{ + matrix.Tvmult(w,v); +} + + +template +template +inline +void +MatrixBlock::Tvmult_add (VECTOR& w, const VECTOR& v) const +{ + matrix.Tvmult_add(w,v); } + //----------------------------------------------------------------------// template