From b3709838106ad07d817e493eed3f34d5f01a1313 Mon Sep 17 00:00:00 2001 From: bangerth Date: Tue, 5 Jan 2010 17:07:42 +0000 Subject: [PATCH] Move the reinit_vector functions into an internal namespace. These are not supposed to be called from user code. git-svn-id: https://svn.dealii.org/trunk@20295 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/multigrid/mg_base.h | 153 ++++++- deal.II/deal.II/include/multigrid/mg_tools.h | 132 ------ .../include/multigrid/mg_transfer.templates.h | 10 +- .../multigrid/mg_transfer_block.templates.h | 17 +- .../mg_transfer_component.templates.h | 22 +- deal.II/deal.II/source/multigrid/mg_base.cc | 315 ++++++++++++- .../deal.II/source/multigrid/mg_dof_tools.cc | 418 +++--------------- 7 files changed, 554 insertions(+), 513 deletions(-) diff --git a/deal.II/deal.II/include/multigrid/mg_base.h b/deal.II/deal.II/include/multigrid/mg_base.h index 87a2185ddd..0c8a42d5e9 100644 --- a/deal.II/deal.II/include/multigrid/mg_base.h +++ b/deal.II/deal.II/include/multigrid/mg_base.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -23,8 +23,155 @@ #include #include + DEAL_II_NAMESPACE_OPEN + +template class MGDoFHandler; +template class MGLevelObject; + + +namespace internal +{ + namespace mg + { + /** + * Ajust vectors on all levels to + * correct size. Here, we just + * count the numbers of degrees + * of freedom on each level and + * @p reinit each level vector + * to this length. + */ + template + static void + reinit_vector (const dealii::MGDoFHandler &mg_dof, + MGLevelObject > &vector); + + /** + * Ajust vectors on all levels to + * correct size. Here, we just + * count the numbers of degrees + * of freedom on each level and + * @p reinit each level vector + * to this length. + */ + template + static void + reinit_vector (const dealii::MGDoFHandler &mg_dof, + MGLevelObject > &vector); + + + /** + * Adjust vectors on all levels + * to correct size. The degrees + * of freedom on each level are + * counted by block and only the + * block selected is used. + */ + template + static void + reinit_vector_by_blocks ( + const dealii::MGDoFHandler &mg_dof, + MGLevelObject > &v, + const unsigned int selected, + std::vector >& cached_sizes); + + /** + * Adjust block vectors on all + * levels to correct size. The + * degrees of freedom on each + * level are counted by block. + */ + template + static void + reinit_vector_by_blocks ( + const dealii::MGDoFHandler &mg_dof, + MGLevelObject > &v, + const std::vector &selected, + std::vector >& cached_sizes); + + /** + * Adjust block-vectors on all + * levels to correct size. Count + * the numbers of degrees of + * freedom on each level + * component-wise. Then, assign + * each block of @p vector the + * corresponding size. + * + * The boolean field @p selected + * allows restricting this + * operation to certain + * components. In this case, @p + * vector will only have as many + * blocks as there are true + * values in @p selected (no + * blocks of length zero are + * padded in). If this argument + * is omitted, all blocks will be + * considered. + * + * Degrees of freedom must be + * sorted by component in order + * to obtain reasonable results + * from this function. + * + * The argument + * @p target_component allows to + * re-sort and group components + * as in + * DoFRenumbering::component_wise. + * + * + */ + template + static void + reinit_vector_by_components (const dealii::MGDoFHandler& mg_dof, + MGLevelObject >& v, + const std::vector& selected, + const std::vector& target_component, + std::vector >& cached_sizes); + /** + * Adjust vectors on all levels + * to correct size. Count the + * numbers of degrees of freedom + * on each level component-wise + * in a single component. Then, + * assign @p vector the + * corresponding size. + * + * The boolean field @p selected + * may be nonzero in a single + * component, indicating the + * block of a block vector the + * argument @p v corresponds to. + * + * Degrees of freedom must be + * sorted by component in order + * to obtain reasonable results + * from this function. + * + * The argument + * @p target_component allows to + * re-sort and group components + * as in + * DoFRenumbering::component_wise. + */ + template + static void + reinit_vector_by_components ( + const dealii::MGDoFHandler &mg_dof, + MGLevelObject > &v, + const std::vector &selected, + const std::vector &target_component, + std::vector >& cached_sizes); + + } +} + + + /*!@addtogroup mg */ /*@{*/ @@ -239,7 +386,7 @@ class MGSmootherBase : public Subscriptor * Release matrices. */ virtual void clear() = 0; - + /** * Smoothing function. This is the * function used in multigrid @@ -247,7 +394,7 @@ class MGSmootherBase : public Subscriptor */ virtual void smooth (const unsigned int level, VECTOR& u, - const VECTOR& rhs) const = 0; + const VECTOR& rhs) const = 0; }; /*@}*/ diff --git a/deal.II/deal.II/include/multigrid/mg_tools.h b/deal.II/deal.II/include/multigrid/mg_tools.h index ecd4c83b83..d421dffcdf 100644 --- a/deal.II/deal.II/include/multigrid/mg_tools.h +++ b/deal.II/deal.II/include/multigrid/mg_tools.h @@ -258,138 +258,6 @@ class MGTools BlockSparseMatrix& matrix, const bool preserve_symmetry); - /** - * Ajust vectors on all levels to - * correct size. Here, we just - * count the numbers of degrees - * of freedom on each level and - * @p reinit each level vector - * to this length. - */ - template - static void - reinit_vector (const MGDoFHandler &mg_dof, - MGLevelObject > &vector); - - /** - * Ajust vectors on all levels to - * correct size. Here, we just - * count the numbers of degrees - * of freedom on each level and - * @p reinit each level vector - * to this length. - */ - template - static void - reinit_vector (const MGDoFHandler &mg_dof, - MGLevelObject > &vector); - - - /** - * Adjust vectors on all levels - * to correct size. The degrees - * of freedom on each level are - * counted by block and only the - * block selected is used. - */ - template - static void - reinit_vector_by_blocks ( - const MGDoFHandler &mg_dof, - MGLevelObject > &v, - const unsigned int selected, - std::vector >& cached_sizes); - - /** - * Adjust block vectors on all - * levels to correct size. The - * degrees of freedom on each - * level are counted by block. - */ - template - static void - reinit_vector_by_blocks ( - const MGDoFHandler &mg_dof, - MGLevelObject > &v, - const std::vector &selected, - std::vector >& cached_sizes); - - /** - * Adjust block-vectors on all - * levels to correct size. Count - * the numbers of degrees of - * freedom on each level - * component-wise. Then, assign - * each block of @p vector the - * corresponding size. - * - * The boolean field @p selected - * allows restricting this - * operation to certain - * components. In this case, @p - * vector will only have as many - * blocks as there are true - * values in @p selected (no - * blocks of length zero are - * padded in). If this argument - * is omitted, all blocks will be - * considered. - * - * Degrees of freedom must be - * sorted by component in order - * to obtain reasonable results - * from this function. - * - * The argument - * @p target_component allows to - * re-sort and group components - * as in - * DoFRenumbering::component_wise. - * - * - */ - template - static void - reinit_vector_by_components (const MGDoFHandler& mg_dof, - MGLevelObject >& v, - const std::vector& selected, - const std::vector& target_component, - std::vector >& cached_sizes); - /** - * Adjust vectors on all levels - * to correct size. Count the - * numbers of degrees of freedom - * on each level component-wise - * in a single component. Then, - * assign @p vector the - * corresponding size. - * - * The boolean field @p selected - * may be nonzero in a single - * component, indicating the - * block of a block vector the - * argument @p v corresponds to. - * - * Degrees of freedom must be - * sorted by component in order - * to obtain reasonable results - * from this function. - * - * The argument - * @p target_component allows to - * re-sort and groupt components - * as in - * DoFRenumbering::component_wise. - */ - template - static void - reinit_vector_by_components ( - const MGDoFHandler &mg_dof, - MGLevelObject > &v, - const std::vector &selected, - const std::vector &target_component, - std::vector >& cached_sizes); - /** * For each level in a multigrid diff --git a/deal.II/deal.II/include/multigrid/mg_transfer.templates.h b/deal.II/deal.II/include/multigrid/mg_transfer.templates.h index 804421ac47..e22a2981d7 100644 --- a/deal.II/deal.II/include/multigrid/mg_transfer.templates.h +++ b/deal.II/deal.II/include/multigrid/mg_transfer.templates.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors +// Copyright (C) 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 @@ -40,7 +40,7 @@ MGTransferPrebuilt::copy_to_mg ( MGLevelObject& dst, const InVector& src) const { - MGTools::reinit_vector(mg_dof_handler, dst); + internal::mg::reinit_vector(mg_dof_handler, dst); bool first = true; for (unsigned int level=mg_dof_handler.get_tria().n_levels();level != 0;) { @@ -115,7 +115,7 @@ MGTransferPrebuilt::copy_from_mg_add ( template template void MGTransferPrebuilt::find_dofs_on_refinement_edges ( - const MGDoFHandler& mg_dof) + const MGDoFHandler& mg_dof) { for (unsigned int level = 0; level::find_dofs_on_refinement_edges ( std::vector boundary_cell_dofs(dofs_per_cell); const unsigned int level = cell->level(); cell->get_mg_dof_indices (local_dof_indices); - for (unsigned int face_nr=0; + for (unsigned int face_nr=0; face_nr::faces_per_cell; ++face_nr) { typename DoFHandler::face_iterator face = cell->face(face_nr); if(!cell->at_boundary(face_nr)) { //interior face - typename MGDoFHandler::cell_iterator neighbor + typename MGDoFHandler::cell_iterator neighbor = cell->neighbor(face_nr); // Do refinement face // from the coarse side diff --git a/deal.II/deal.II/include/multigrid/mg_transfer_block.templates.h b/deal.II/deal.II/include/multigrid/mg_transfer_block.templates.h index 5ead161738..f9a903efac 100644 --- a/deal.II/deal.II/include/multigrid/mg_transfer_block.templates.h +++ b/deal.II/deal.II/include/multigrid/mg_transfer_block.templates.h @@ -1,8 +1,8 @@ //--------------------------------------------------------------------------- // mg_transfer.templates.h,v 1.22 2006/01/29 15:03:55 guido Exp -// Version: +// Version: // -// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors +// Copyright (C) 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 @@ -31,7 +31,7 @@ DEAL_II_NAMESPACE_OPEN // Simplify some things below typedef std::map::const_iterator IT; - + template @@ -42,7 +42,8 @@ MGTransferBlockSelect::copy_to_mg ( MGLevelObject > &dst, const BlockVector &src) const { - MGTools::reinit_vector_by_blocks(mg_dof_handler, dst, selected_block, sizes); + internal::mg::reinit_vector_by_blocks(mg_dof_handler, dst, + selected_block, sizes); // For MGTransferBlockSelect, the // multilevel block is always the // first, since only one block is @@ -70,7 +71,8 @@ MGTransferBlockSelect::copy_to_mg ( MGLevelObject > &dst, const Vector &src) const { - MGTools::reinit_vector_by_blocks(mg_dof_handler, dst, selected_block, sizes); + internal::mg::reinit_vector_by_blocks(mg_dof_handler, dst, + selected_block, sizes); // For MGTransferBlockSelect, the // multilevel block is always the // first, since only one block is selected. @@ -84,7 +86,7 @@ MGTransferBlockSelect::copy_to_mg ( if (!first) restrict_and_add (level+1, dst[level], dst[level+1]); first = false; - } + } } @@ -173,7 +175,8 @@ MGTransferBlock::copy_to_mg ( MGLevelObject >& dst, const BlockVector& src) const { - MGTools::reinit_vector_by_blocks(mg_dof_handler, dst, selected, sizes); + internal::mg::reinit_vector_by_blocks(mg_dof_handler, dst, + selected, sizes); bool first = true; for (unsigned int level=mg_dof_handler.get_tria().n_levels();level != 0;) { diff --git a/deal.II/deal.II/include/multigrid/mg_transfer_component.templates.h b/deal.II/deal.II/include/multigrid/mg_transfer_component.templates.h index ce7bc81818..65a2b9c0f2 100644 --- a/deal.II/deal.II/include/multigrid/mg_transfer_component.templates.h +++ b/deal.II/deal.II/include/multigrid/mg_transfer_component.templates.h @@ -1,8 +1,8 @@ //--------------------------------------------------------------------------- // mg_transfer.templates.h,v 1.22 2006/01/29 15:03:55 guido Exp -// Version: +// Version: // -// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors +// Copyright (C) 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 @@ -67,25 +67,25 @@ MGTransferSelect::do_copy_to_mg ( const unsigned int offset) const { const FiniteElement& fe = mg_dof_handler.get_fe(); - + const unsigned int dofs_per_cell = fe.dofs_per_cell; - + // set the elements of the vectors // on all levels to zero unsigned int minlevel = dst.get_minlevel(); unsigned int maxlevel = dst.get_maxlevel(); - + dst=0; - + Assert(sizes.size()==mg_dof_handler.get_tria().n_levels(), ExcMatricesNotBuilt()); - MGTools::reinit_vector_by_components( + internal::mg::reinit_vector_by_components( mg_dof_handler, dst, mg_selected, mg_target_component, sizes); - + std::vector global_dof_indices (dofs_per_cell); std::vector level_dof_indices (dofs_per_cell); - + // Build a vector of the selected // indices, since traversing all // indices on each cell is too @@ -133,7 +133,7 @@ MGTransferSelect::do_copy_to_mg ( = mg_component_start[level][selected_component]; const typename std::vector::const_iterator end = selected_indices.end(); - + for (typename std::vector::const_iterator i=selected_indices.begin(); i != end; ++i) @@ -221,7 +221,7 @@ MGTransferSelect::do_copy_from_mg ( const unsigned int offset) const { const FiniteElement& fe = mg_dof_handler.get_fe(); - + const unsigned int dofs_per_cell = fe.dofs_per_cell; std::vector global_dof_indices (dofs_per_cell); std::vector level_dof_indices (dofs_per_cell); diff --git a/deal.II/deal.II/source/multigrid/mg_base.cc b/deal.II/deal.II/source/multigrid/mg_base.cc index 8b011beb66..ab2cc9bec1 100644 --- a/deal.II/deal.II/source/multigrid/mg_base.cc +++ b/deal.II/deal.II/source/multigrid/mg_base.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2005, 2006 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 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 @@ -12,13 +12,250 @@ //--------------------------------------------------------------------------- -#include #include #include +#include +#include +#include +#include +#include + +#include DEAL_II_NAMESPACE_OPEN +namespace internal +{ + namespace mg + { + template + void + reinit_vector (const dealii::MGDoFHandler &mg_dof, + MGLevelObject > &v) + { + for (unsigned int level=v.get_minlevel(); + level<=v.get_maxlevel();++level) + { + unsigned int n = mg_dof.n_dofs (level); + v[level].reinit(n); + } + + } + + + template + void + reinit_vector (const dealii::MGDoFHandler &mg_dof, + MGLevelObject > &v) + { + const unsigned int n_blocks = mg_dof.get_fe().n_blocks(); + std::vector > + ndofs(mg_dof.get_tria().n_levels(), + std::vector(n_blocks)); + MGTools::count_dofs_per_block (mg_dof, ndofs); + + for (unsigned int level=v.get_minlevel(); + level<=v.get_maxlevel();++level) + { + v[level].reinit(ndofs[level]); + } + } + + + template + void + reinit_vector_by_components ( + const dealii::MGDoFHandler &mg_dof, + MGLevelObject > &v, + const std::vector &sel, + const std::vector &target_comp, + std::vector >& ndofs) + { + std::vector selected=sel; + std::vector target_component=target_comp; + const unsigned int ncomp = mg_dof.get_fe().n_components(); + + // If the selected and + // target_component have size 0, + // they must be replaced by default + // values. + // + // Since we already made copies + // directly after this function was + // called, we use the arguments + // directly. + if (target_component.size() == 0) + { + target_component.resize(ncomp); + for (unsigned int i=0;i > + new_dofs(mg_dof.get_tria().n_levels(), + std::vector(target_component.size())); + std::swap(ndofs, new_dofs); + MGTools::count_dofs_per_component (mg_dof, ndofs, + true, target_component); + } + + for (unsigned int level=v.get_minlevel(); + level<=v.get_maxlevel();++level) + { + v[level].reinit(n_selected, 0); + unsigned int k=0; + for (unsigned int i=0;i + void + reinit_vector_by_components ( + const dealii::MGDoFHandler &mg_dof, + MGLevelObject > &v, + const std::vector &selected, + const std::vector &target_component, + std::vector >& ndofs) + { + Assert (selected.size() == target_component.size(), + ExcDimensionMismatch(selected.size(), target_component.size())); + + // Compute the number of blocks needed +#ifdef DEBUG + const unsigned int n_selected + = std::accumulate(selected.begin(), + selected.end(), + 0U); + Assert(n_selected == 1, ExcDimensionMismatch(n_selected, 1)); +#endif + + unsigned int selected_block = 0; + while (!selected[selected_block]) + ++selected_block; + + if (ndofs.size() == 0) + { + std::vector > + new_dofs(mg_dof.get_tria().n_levels(), + std::vector(target_component.size())); + std::swap(ndofs, new_dofs); + MGTools::count_dofs_per_component (mg_dof, ndofs, + true, target_component); + } + + for (unsigned int level=v.get_minlevel(); + level<=v.get_maxlevel();++level) + { + v[level].reinit(ndofs[level][selected_block]); + } + } + + + template + void + reinit_vector_by_blocks ( + const dealii::MGDoFHandler &mg_dof, + MGLevelObject > &v, + const std::vector &sel, + std::vector >& ndofs) + { + std::vector selected=sel; + // Compute the number of blocks needed + const unsigned int n_selected + = std::accumulate(selected.begin(), + selected.end(), + 0U); + + if (ndofs.size() == 0) + { + std::vector > + new_dofs(mg_dof.get_tria().n_levels(), + std::vector(selected.size())); + std::swap(ndofs, new_dofs); + MGTools::count_dofs_per_block (mg_dof, ndofs); + } + + for (unsigned int level=v.get_minlevel(); + level<=v.get_maxlevel();++level) + { + v[level].reinit(n_selected, 0); + unsigned int k=0; + for (unsigned int i=0;i + void + reinit_vector_by_blocks ( + const dealii::MGDoFHandler &mg_dof, + MGLevelObject > &v, + const unsigned int selected_block, + std::vector >& ndofs) + { + const unsigned int n_blocks = mg_dof.get_fe().n_blocks(); + Assert(selected_block < n_blocks, ExcIndexRange(selected_block, 0, n_blocks)); + + std::vector selected(n_blocks, false); + selected[selected_block] = true; + + if (ndofs.size() == 0) + { + std::vector > + new_dofs(mg_dof.get_tria().n_levels(), + std::vector(selected.size())); + std::swap(ndofs, new_dofs); + MGTools::count_dofs_per_block (mg_dof, ndofs); + } + + for (unsigned int level=v.get_minlevel(); + level<=v.get_maxlevel();++level) + { + v[level].reinit(ndofs[level][selected_block]); + } + } + } +} + + + template MGTransferBase::~MGTransferBase() {} @@ -41,23 +278,81 @@ MGCoarseGridBase::~MGCoarseGridBase() // Explicit instantiations -template class MGTransferBase >; -template class MGTransferBase >; +//TODO: Use the template expander script for this +namespace internal +{ + namespace mg + { + template void reinit_vector ( + const dealii::MGDoFHandler&, + MGLevelObject >&); + template void reinit_vector ( + const dealii::MGDoFHandler&, + MGLevelObject >&); + + template void reinit_vector ( + const dealii::MGDoFHandler&, + MGLevelObject >&); + template void reinit_vector ( + const dealii::MGDoFHandler&, + MGLevelObject >&); + + template void reinit_vector_by_components ( + const dealii::MGDoFHandler&, + MGLevelObject >&, + const std::vector &, + const std::vector &, + std::vector >&); + template void reinit_vector_by_components ( + const dealii::MGDoFHandler&, + MGLevelObject >&, + const std::vector &, + const std::vector &, + std::vector >&); + + template void reinit_vector_by_components ( + const dealii::MGDoFHandler&, MGLevelObject >&, + const std::vector&, const std::vector&, + std::vector >&); + template void reinit_vector_by_components ( + const dealii::MGDoFHandler&, MGLevelObject >&, + const std::vector&, const std::vector&, + std::vector >&); + + template void reinit_vector_by_blocks ( + const dealii::MGDoFHandler&, MGLevelObject >&, + const std::vector &, std::vector >&); + template void reinit_vector_by_blocks ( + const dealii::MGDoFHandler&, MGLevelObject >&, + const std::vector &, std::vector >&); + + template void reinit_vector_by_blocks ( + const dealii::MGDoFHandler&, MGLevelObject >&, + unsigned int, std::vector >&); + template void reinit_vector_by_blocks ( + const dealii::MGDoFHandler&, MGLevelObject >&, + unsigned int, std::vector >&); + } +} + + +template class MGTransferBase >; +template class MGTransferBase >; template class MGTransferBase >; template class MGTransferBase >; -template class MGMatrixBase >; -template class MGMatrixBase >; +template class MGMatrixBase >; +template class MGMatrixBase >; template class MGMatrixBase >; template class MGMatrixBase >; -template class MGSmootherBase >; -template class MGSmootherBase >; +template class MGSmootherBase >; +template class MGSmootherBase >; template class MGSmootherBase >; template class MGSmootherBase >; -template class MGCoarseGridBase >; -template class MGCoarseGridBase >; +template class MGCoarseGridBase >; +template class MGCoarseGridBase >; template class MGCoarseGridBase >; template class MGCoarseGridBase >; diff --git a/deal.II/deal.II/source/multigrid/mg_dof_tools.cc b/deal.II/deal.II/source/multigrid/mg_dof_tools.cc index 1ebfe6fb2c..62ccc8ae2f 100644 --- a/deal.II/deal.II/source/multigrid/mg_dof_tools.cc +++ b/deal.II/deal.II/source/multigrid/mg_dof_tools.cc @@ -93,7 +93,7 @@ MGTools::compute_row_length_vector( const_cast&> (dofs.get_tria()); user_flags_triangulation.save_user_flags(old_flags); user_flags_triangulation.clear_user_flags(); - + const typename MGDoFHandler::cell_iterator end = dofs.end(level); typename MGDoFHandler::active_cell_iterator cell; std::vector cell_indices; @@ -133,7 +133,7 @@ MGTools::compute_row_length_vector( // dofs_per_vertex is // arbitrary, not the other way // round. -//TODO: This assumes that even in hp context, the dofs per face coincide! +//TODO: This assumes that even in hp context, the dofs per face coincide! unsigned int increment = fe.dofs_per_cell - dim * fe.dofs_per_face; while (i < fe.first_line_index) row_lengths[cell_indices[i++]] += increment; @@ -143,7 +143,7 @@ MGTools::compute_row_length_vector( // faces are handled below, we // have to subtract ALL faces // in this case. - + // In all other cases we // subtract adjacent faces to be // added in the loop below. @@ -152,7 +152,7 @@ MGTools::compute_row_length_vector( : fe.dofs_per_cell - GeometryInfo::faces_per_cell * fe.dofs_per_face; while (i < fe.first_quad_index) row_lengths[cell_indices[i++]] += increment; - + // Now quads in 2D and 3D increment = (dim>2) ? fe.dofs_per_cell - (dim-2) * fe.dofs_per_face @@ -170,7 +170,7 @@ MGTools::compute_row_length_vector( // coupled topologically to the // adjacent cells, but we // subtracted some faces. - + // Now, let's go by the faces // and add the missing // contribution as well as the @@ -185,21 +185,21 @@ MGTools::compute_row_length_vector( if (static_cast(neighbor->level()) != level) level_boundary = true; } - + if (level_boundary) { for (unsigned int i=0;i& nfe = neighbor->get_fe(); typename MGDoFHandler::face_iterator face = cell->face(iface); - + // Flux couplings are // computed from both sides // for simplicity. - + // The dofs on the common face // will be handled below, // therefore, we subtract them @@ -210,7 +210,7 @@ MGTools::compute_row_length_vector( for (unsigned int i=0;iuser_flag_set()) @@ -223,7 +223,7 @@ MGTools::compute_row_length_vector( // face dofs. Since we // subtracted two faces, we // have to re-add one. - + // If one side of the face // is refined, all the fine // face dofs couple with @@ -252,7 +252,7 @@ MGTools::compute_row_length_vector( { Assert (row_lengths.size() == dofs.n_dofs(), ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs())); - + // Function starts here by // resetting the counters. std::fill(row_lengths.begin(), row_lengths.end(), 0); @@ -267,12 +267,12 @@ MGTools::compute_row_length_vector( const_cast&> (dofs.get_tria()); user_flags_triangulation.save_user_flags(old_flags); user_flags_triangulation.clear_user_flags(); - + const typename MGDoFHandler::cell_iterator end = dofs.end(level); typename MGDoFHandler::active_cell_iterator cell; std::vector cell_indices; std::vector neighbor_indices; - + // We have to translate the // couplings from components to // blocks, so this works for @@ -281,7 +281,7 @@ MGTools::compute_row_length_vector( std::vector > couple_face; DoFTools::convert_couplings_to_blocks(dofs, couplings, couple_cell); DoFTools::convert_couplings_to_blocks(dofs, flux_couplings, couple_face); - + // We loop over cells and go from // cells to lower dimensional // objects. This is the only way to @@ -293,7 +293,7 @@ MGTools::compute_row_length_vector( { const FiniteElement& fe = cell->get_fe(); const unsigned int fe_index = cell->active_fe_index(); - + Assert (couplings.n_rows()==fe.n_components(), ExcDimensionMismatch(couplings.n_rows(), fe.n_components())); Assert (couplings.n_cols()==fe.n_components(), @@ -302,7 +302,7 @@ MGTools::compute_row_length_vector( ExcDimensionMismatch(flux_couplings.n_rows(), fe.n_components())); Assert (flux_couplings.n_cols()==fe.n_components(), ExcDimensionMismatch(flux_couplings.n_cols(), fe.n_components())); - + cell_indices.resize(fe.dofs_per_cell); cell->get_mg_dof_indices(cell_indices); unsigned int i = 0; @@ -347,7 +347,7 @@ MGTools::compute_row_length_vector( // faces are handled below, we // have to subtract ALL faces // in this case. - + // In all other cases we // subtract adjacent faces to be // added in the loop below. @@ -367,7 +367,7 @@ MGTools::compute_row_length_vector( } ++i; } - + // Now quads in 2D and 3D while (i < fe.first_hex_index) { @@ -385,7 +385,7 @@ MGTools::compute_row_length_vector( } ++i; } - + // Finally, cells in 3D while (i < fe.dofs_per_cell) { @@ -401,14 +401,14 @@ MGTools::compute_row_length_vector( } ++i; } - + // At this point, we have // counted all dofs // contributiong from cells // coupled topologically to the // adjacent cells, but we // subtracted some faces. - + // Now, let's go by the faces // and add the missing // contribution as well as the @@ -423,21 +423,21 @@ MGTools::compute_row_length_vector( if (static_cast(neighbor->level()) != level) level_boundary = true; } - + if (level_boundary) { for (unsigned int i=0;i& nfe = neighbor->get_fe(); typename MGDoFHandler::face_iterator face = cell->face(iface); - + // Flux couplings are // computed from both sides // for simplicity. - + // The dofs on the common face // will be handled below, // therefore, we subtract them @@ -452,7 +452,7 @@ MGTools::compute_row_length_vector( - nfe.base_element(base).dofs_per_face; row_lengths[cell_indices[i]] += increment; } - + // Do this only once per // face and not on the // hanging faces. @@ -466,12 +466,12 @@ MGTools::compute_row_length_vector( // face dofs. Since we // subtracted two faces, we // have to re-add one. - + // If one side of the face // is refined, all the fine // face dofs couple with // the coarse one. - + // Wolfgang, do they couple // with each other by // constraints? @@ -519,7 +519,7 @@ void MGTools::make_sparsity_pattern ( std::vector dofs_on_this_cell(dofs_per_cell); typename MGDoFHandler::cell_iterator cell = dof.begin(level), endc = dof.end(level); - for (; cell!=endc; ++cell) + for (; cell!=endc; ++cell) { cell->get_mg_dof_indices (dofs_on_this_cell); // make sparsity pattern for this cell @@ -540,7 +540,7 @@ MGTools::make_flux_sparsity_pattern ( const unsigned int level) { const unsigned int n_dofs = dof.n_dofs(level); - + Assert (sparsity.n_rows() == n_dofs, ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); Assert (sparsity.n_cols() == n_dofs, @@ -583,7 +583,7 @@ MGTools::make_flux_sparsity_pattern ( } } } - } + } } } @@ -601,7 +601,7 @@ MGTools::make_flux_sparsity_pattern_edge ( const unsigned int fine_dofs = dof.n_dofs(level); const unsigned int coarse_dofs = dof.n_dofs(level-1); - + // Matrix maps from fine level to coarse level Assert (sparsity.n_rows() == coarse_dofs, @@ -642,7 +642,7 @@ MGTools::make_flux_sparsity_pattern_edge ( } } } - } + } } } @@ -663,7 +663,7 @@ MGTools::make_flux_sparsity_pattern ( const FiniteElement& fe = dof.get_fe(); const unsigned int n_dofs = dof.n_dofs(level); const unsigned int n_comp = fe.n_components(); - + Assert (sparsity.n_rows() == n_dofs, ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); Assert (sparsity.n_cols() == n_dofs, @@ -676,23 +676,23 @@ MGTools::make_flux_sparsity_pattern ( ExcDimensionMismatch (flux_mask.n_rows(), n_comp)); Assert (flux_mask.n_cols() == n_comp, ExcDimensionMismatch (flux_mask.n_cols(), n_comp)); - + const unsigned int total_dofs = fe.dofs_per_cell; std::vector dofs_on_this_cell(total_dofs); std::vector dofs_on_other_cell(total_dofs); Table<2,bool> support_on_face(total_dofs, GeometryInfo::faces_per_cell); - + typename MGDoFHandler::cell_iterator cell = dof.begin(level), endc = dof.end(level); const Table<2,DoFTools::Coupling> int_dof_mask = DoFTools::dof_couplings_from_component_couplings(fe, int_mask), flux_dof_mask = DoFTools::dof_couplings_from_component_couplings(fe, flux_mask); - + for (unsigned int i=0; i::faces_per_cell;++f) support_on_face(i,f) = fe.has_support_on_face(i,f); - + // Clear user flags because we will // need them. But first we save // them and make sure that we @@ -704,7 +704,7 @@ MGTools::make_flux_sparsity_pattern ( std::vector user_flags; dof.get_tria().save_user_flags(user_flags); const_cast &>(dof.get_tria()).clear_user_flags (); - + for (; cell!=endc; ++cell) { cell->get_mg_dof_indices (dofs_on_this_cell); @@ -732,7 +732,7 @@ MGTools::make_flux_sparsity_pattern ( for (unsigned int j=0; j::cell_iterator neighbor = cell->neighbor(face); - + if (neighbor->level() < cell->level()) continue; - + unsigned int neighbor_face = cell->neighbor_of_neighbor(face); - + neighbor->get_mg_dof_indices (dofs_on_other_cell); for (unsigned int i=0; iface(neighbor_face)->set_user_flag (); + neighbor->face(neighbor_face)->set_user_flag (); } } } - + // finally restore the user flags const_cast &>(dof.get_tria()).load_user_flags(user_flags); } @@ -838,13 +838,13 @@ MGTools::make_flux_sparsity_pattern_edge ( { const FiniteElement& fe = dof.get_fe(); const unsigned int n_comp = fe.n_components(); - + Assert ((level>=1) && (level dofs_on_this_cell(dofs_per_cell); std::vector dofs_on_other_cell(dofs_per_cell); Table<2,bool> support_on_face(dofs_per_cell, GeometryInfo::faces_per_cell); - + typename MGDoFHandler::cell_iterator cell = dof.begin(level), endc = dof.end(level); const Table<2,DoFTools::Coupling> flux_dof_mask = DoFTools::dof_couplings_from_component_couplings(fe, flux_mask); - + for (unsigned int i=0; i::faces_per_cell;++f) support_on_face(i,f) = fe.has_support_on_face(i,f); - + for (; cell!=endc; ++cell) { cell->get_mg_dof_indices (dofs_on_this_cell); @@ -902,7 +902,7 @@ MGTools::make_flux_sparsity_pattern_edge ( } } } - } + } } } @@ -918,10 +918,10 @@ count_dofs_per_component (const MGDoFHandler &dof_handler, const FiniteElement& fe = dof_handler.get_fe(); const unsigned int n_components = fe.n_components(); const unsigned int nlevels = dof_handler.get_tria().n_levels(); - + Assert (result.size() == nlevels, ExcDimensionMismatch(result.size(), nlevels)); - + if (target_component.size() == 0) { target_component.resize(n_components); @@ -936,7 +936,7 @@ count_dofs_per_component (const MGDoFHandler &dof_handler, { result[l].resize (n_components); std::fill (result[l].begin(),result[l].end(), 0U); - + // special case for only one // component. treat this first // since it does not require any @@ -979,7 +979,7 @@ count_dofs_per_component (const MGDoFHandler &dof_handler, const FiniteElement& base = fe.base_element(b); // Dimension of base element unsigned int d = base.n_components(); - + for (unsigned int m=0;m& fe = dof_handler.get_fe(); const unsigned int n_blocks = fe.n_blocks(); const unsigned int n_levels = dof_handler.get_tria().n_levels(); - + dofs_per_block.resize (n_levels); - + for (unsigned int l=0;l::type::const_iterator i=function_map.begin(); // i!=function_map.end(); ++i) // Assert (n_components == i->second->n_components, @@ -1154,7 +1153,7 @@ MGTools::make_boundary_list( std::vector face_dofs; face_dofs.reserve (DoFTools::max_dofs_per_face(dof)); std::fill (face_dofs.begin (), face_dofs.end (), DoFHandler::invalid_dof_index); - + typename MGDoFHandler::cell_iterator cell = dof.begin(), endc = dof.end(); for (; cell!=endc; ++cell) @@ -1301,227 +1300,6 @@ MGTools::make_boundary_list( #endif - -template -void -MGTools::reinit_vector (const MGDoFHandler &mg_dof, - MGLevelObject > &v) -{ - for (unsigned int level=v.get_minlevel(); - level<=v.get_maxlevel();++level) - { - unsigned int n = mg_dof.n_dofs (level); - v[level].reinit(n); - } - -} - - -template -void -MGTools::reinit_vector (const MGDoFHandler &mg_dof, - MGLevelObject > &v) -{ - const unsigned int n_blocks = mg_dof.get_fe().n_blocks(); - std::vector > - ndofs(mg_dof.get_tria().n_levels(), - std::vector(n_blocks)); - count_dofs_per_block (mg_dof, ndofs); - - for (unsigned int level=v.get_minlevel(); - level<=v.get_maxlevel();++level) - { - v[level].reinit(ndofs[level]); - } -} - - -template -void -MGTools::reinit_vector_by_components ( - const MGDoFHandler &mg_dof, - MGLevelObject > &v, - const std::vector &sel, - const std::vector &target_comp, - std::vector >& ndofs) -{ - std::vector selected=sel; - std::vector target_component=target_comp; - const unsigned int ncomp = mg_dof.get_fe().n_components(); - - // If the selected and - // target_component have size 0, - // they must be replaced by default - // values. - // - // Since we already made copies - // directly after this function was - // called, we use the arguments - // directly. - if (target_component.size() == 0) - { - target_component.resize(ncomp); - for (unsigned int i=0;i > - new_dofs(mg_dof.get_tria().n_levels(), - std::vector(target_component.size())); - std::swap(ndofs, new_dofs); - count_dofs_per_component (mg_dof, ndofs, true, target_component); - } - - for (unsigned int level=v.get_minlevel(); - level<=v.get_maxlevel();++level) - { - v[level].reinit(n_selected, 0); - unsigned int k=0; - for (unsigned int i=0;i -void -MGTools::reinit_vector_by_components ( - const MGDoFHandler &mg_dof, - MGLevelObject > &v, - const std::vector &selected, - const std::vector &target_component, - std::vector >& ndofs) -{ - Assert (selected.size() == target_component.size(), - ExcDimensionMismatch(selected.size(), target_component.size())); - - // Compute the number of blocks needed -#ifdef DEBUG - const unsigned int n_selected - = std::accumulate(selected.begin(), - selected.end(), - 0U); - Assert(n_selected == 1, ExcDimensionMismatch(n_selected, 1)); -#endif - - unsigned int selected_block = 0; - while (!selected[selected_block]) - ++selected_block; - - if (ndofs.size() == 0) - { - std::vector > - new_dofs(mg_dof.get_tria().n_levels(), - std::vector(target_component.size())); - std::swap(ndofs, new_dofs); - count_dofs_per_component (mg_dof, ndofs, true, target_component); - } - - for (unsigned int level=v.get_minlevel(); - level<=v.get_maxlevel();++level) - { - v[level].reinit(ndofs[level][selected_block]); - } -} - - -template -void -MGTools::reinit_vector_by_blocks ( - const MGDoFHandler &mg_dof, - MGLevelObject > &v, - const std::vector &sel, - std::vector >& ndofs) -{ - std::vector selected=sel; - // Compute the number of blocks needed - const unsigned int n_selected - = std::accumulate(selected.begin(), - selected.end(), - 0U); - - if (ndofs.size() == 0) - { - std::vector > - new_dofs(mg_dof.get_tria().n_levels(), - std::vector(selected.size())); - std::swap(ndofs, new_dofs); - count_dofs_per_block (mg_dof, ndofs); - } - - for (unsigned int level=v.get_minlevel(); - level<=v.get_maxlevel();++level) - { - v[level].reinit(n_selected, 0); - unsigned int k=0; - for (unsigned int i=0;i -void -MGTools::reinit_vector_by_blocks ( - const MGDoFHandler &mg_dof, - MGLevelObject > &v, - const unsigned int selected_block, - std::vector >& ndofs) -{ - const unsigned int n_blocks = mg_dof.get_fe().n_blocks(); - Assert(selected_block < n_blocks, ExcIndexRange(selected_block, 0, n_blocks)); - - std::vector selected(n_blocks, false); - selected[selected_block] = true; - - if (ndofs.size() == 0) - { - std::vector > - new_dofs(mg_dof.get_tria().n_levels(), - std::vector(selected.size())); - std::swap(ndofs, new_dofs); - count_dofs_per_block (mg_dof, ndofs); - } - - for (unsigned int level=v.get_minlevel(); - level<=v.get_maxlevel();++level) - { - v[level].reinit(ndofs[level][selected_block]); - } -} @@ -1546,7 +1324,7 @@ extract_inner_interface_dofs (const MGDoFHandler &mg_dof_handler, mg_dof_handler.n_dofs(l))); Assert (boundary_interface_dofs[l].size() == mg_dof_handler.n_dofs(l), ExcDimensionMismatch (boundary_interface_dofs[l].size(), - mg_dof_handler.n_dofs(l))); + mg_dof_handler.n_dofs(l))); std::fill (interface_dofs[l].begin(), interface_dofs[l].end(), @@ -1555,9 +1333,9 @@ extract_inner_interface_dofs (const MGDoFHandler &mg_dof_handler, boundary_interface_dofs[l].end(), false); } - + const FiniteElement &fe = mg_dof_handler.get_fe(); - + const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int dofs_per_face = fe.dofs_per_face; @@ -1566,17 +1344,17 @@ extract_inner_interface_dofs (const MGDoFHandler &mg_dof_handler, std::vector cell_dofs(dofs_per_cell, false); std::vector boundary_cell_dofs(dofs_per_cell, false); - + typename MGDoFHandler::cell_iterator cell = mg_dof_handler.begin(), endc = mg_dof_handler.end(); for (; cell!=endc; ++cell) { bool has_coarser_neighbor = false; - + std::fill (cell_dofs.begin(), cell_dofs.end(), false); std::fill (boundary_cell_dofs.begin(), boundary_cell_dofs.end(), false); - + for (unsigned int face_nr=0; face_nr::faces_per_cell; ++face_nr) { const typename DoFHandler::face_iterator face = cell->face(face_nr); @@ -1605,7 +1383,7 @@ extract_inner_interface_dofs (const MGDoFHandler &mg_dof_handler, if (cell_dofs[fe.face_to_cell_index(j,face_nr)] == true) boundary_cell_dofs[fe.face_to_cell_index(j,face_nr)] = true; - + const unsigned int level = cell->level(); cell->get_mg_dof_indices (local_dof_indices); @@ -1675,56 +1453,6 @@ MGTools::compute_row_length_vector( std::vector&, const Table<2,DoFTools::Coupling>&, const Table<2,DoFTools::Coupling>&); -template void MGTools::reinit_vector ( - const MGDoFHandler&, - MGLevelObject >&); -template void MGTools::reinit_vector ( - const MGDoFHandler&, - MGLevelObject >&); - -template void MGTools::reinit_vector ( - const MGDoFHandler&, - MGLevelObject >&); -template void MGTools::reinit_vector ( - const MGDoFHandler&, - MGLevelObject >&); - -template void MGTools::reinit_vector_by_components ( - const MGDoFHandler&, - MGLevelObject >&, - const std::vector &, - const std::vector &, - std::vector >&); -template void MGTools::reinit_vector_by_components ( - const MGDoFHandler&, - MGLevelObject >&, - const std::vector &, - const std::vector &, - std::vector >&); - -template void MGTools::reinit_vector_by_components ( - const MGDoFHandler&, MGLevelObject >&, - const std::vector&, const std::vector&, - std::vector >&); -template void MGTools::reinit_vector_by_components ( - const MGDoFHandler&, MGLevelObject >&, - const std::vector&, const std::vector&, - std::vector >&); - -template void MGTools::reinit_vector_by_blocks ( - const MGDoFHandler&, MGLevelObject >&, - const std::vector &, std::vector >&); -template void MGTools::reinit_vector_by_blocks ( - const MGDoFHandler&, MGLevelObject >&, - const std::vector &, std::vector >&); - -template void MGTools::reinit_vector_by_blocks ( - const MGDoFHandler&, MGLevelObject >&, - unsigned int, std::vector >&); -template void MGTools::reinit_vector_by_blocks ( - const MGDoFHandler&, MGLevelObject >&, - unsigned int, std::vector >&); - template void MGTools::count_dofs_per_component ( const MGDoFHandler&, std::vector >&, bool, std::vector); -- 2.39.5