From e31c3c1731632fa8f85bcde9713844f1ae6edea5 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 6 Jan 2010 09:32:48 +0000 Subject: [PATCH] Move functions into .cc files that are only used in these very files. git-svn-id: https://svn.dealii.org/trunk@20302 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/multigrid/mg_base.h | 148 -------- .../include/multigrid/mg_transfer.templates.h | 31 -- .../multigrid/mg_transfer_block.templates.h | 85 ----- .../mg_transfer_component.templates.h | 100 ------ deal.II/deal.II/source/multigrid/mg_base.cc | 311 ----------------- .../source/multigrid/mg_transfer_block.cc | 212 ++++++++++-- .../source/multigrid/mg_transfer_component.cc | 318 ++++++++++++++++-- .../source/multigrid/mg_transfer_prebuilt.cc | 110 ++++++ 8 files changed, 599 insertions(+), 716 deletions(-) diff --git a/deal.II/deal.II/include/multigrid/mg_base.h b/deal.II/deal.II/include/multigrid/mg_base.h index 8bea78bf4e..de6ea59612 100644 --- a/deal.II/deal.II/include/multigrid/mg_base.h +++ b/deal.II/deal.II/include/multigrid/mg_base.h @@ -31,154 +31,6 @@ 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. - * For compatibility reasons with - * the next function - * the target_component is added - * here but is not used. - */ - template - static void - reinit_vector (const dealii::MGDoFHandler &mg_dof, - std::vector target_component, - 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. The target_component - * is handed to MGTools::count_dofs_per_block. - * See for documentation there. - */ - template - static void - reinit_vector (const dealii::MGDoFHandler &mg_dof, - std::vector target_component, - 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 */ /*@{*/ 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 a756eff5fb..6ddff57548 100644 --- a/deal.II/deal.II/include/multigrid/mg_transfer.templates.h +++ b/deal.II/deal.II/include/multigrid/mg_transfer.templates.h @@ -29,37 +29,6 @@ DEAL_II_NAMESPACE_OPEN /* --------------------- MGTransferPrebuilt -------------- */ -template -template -void -MGTransferPrebuilt::copy_to_mg ( - const MGDoFHandler& mg_dof_handler, - MGLevelObject& dst, - const InVector& src) const -{ - internal::mg::reinit_vector(mg_dof_handler, component_to_block_map, dst); - bool first = true; - for (unsigned int level=mg_dof_handler.get_tria().n_levels();level != 0;) - { - --level; - - typedef std::map::const_iterator IT; - for (IT i= copy_indices[level].begin(); - i != copy_indices[level].end();++i) - dst[level](i->second) = src(i->first); - // For non-DG: degrees of - // freedom in the refinement - // face may need special - // attention, since they belong - // to the coarse level, but - // have fine level basis - // functions - if (!first) - restrict_and_add (level+1, dst[level], dst[level+1]); - first = false; - } -} - template 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 f9a903efac..b80d4d6528 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 @@ -34,63 +34,6 @@ typedef std::map::const_iterator IT; -template -template -void -MGTransferBlockSelect::copy_to_mg ( - const MGDoFHandler &mg_dof_handler, - MGLevelObject > &dst, - const BlockVector &src) const -{ - 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. - bool first = true; - for (unsigned int level=mg_dof_handler.get_tria().n_levels();level != 0;) - { - --level; - for (IT i= copy_indices[selected_block][level].begin(); - i != copy_indices[selected_block][level].end();++i) - dst[level](i->second) = src.block(selected_block)(i->first); - if (!first) - restrict_and_add (level+1, dst[level], dst[level+1]); - first = false; - } -} - - - -template -template -void -MGTransferBlockSelect::copy_to_mg ( - const MGDoFHandler &mg_dof_handler, - MGLevelObject > &dst, - const Vector &src) const -{ - 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. - bool first = true; - for (unsigned int level=mg_dof_handler.get_tria().n_levels();level != 0;) - { - --level; - for (IT i= copy_indices[selected_block][level].begin(); - i != copy_indices[selected_block][level].end();++i) - dst[level](i->second) = src(i->first); - if (!first) - restrict_and_add (level+1, dst[level], dst[level+1]); - first = false; - } -} - - - template template void @@ -167,34 +110,6 @@ MGTransferBlockSelect::memory_consumption () const -template -template -void -MGTransferBlock::copy_to_mg ( - const MGDoFHandler& mg_dof_handler, - MGLevelObject >& dst, - const BlockVector& src) const -{ - 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;) - { - --level; - for (unsigned int block=0;blocksecond) = src.block(block)(i->first); - if (!first) - restrict_and_add (level+1, dst[level], dst[level+1]); - first = false; - } -} - - - - template template void 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 65a2b9c0f2..4765f323c9 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 @@ -57,106 +57,6 @@ MGTransferSelect::copy_to_mg ( -template -template -void -MGTransferSelect::do_copy_to_mg ( - const MGDoFHandler& mg_dof_handler, - MGLevelObject >& dst, - const InVector& src, - 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()); - - 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 - // slow. - std::vector selected_indices; - selected_indices.reserve(dofs_per_cell); - for (unsigned int i=0; i=static_cast(minlevel); --level) - { - - typename MGDoFHandler::active_cell_iterator - level_cell = mg_dof_handler.begin_active(level); - const typename MGDoFHandler::active_cell_iterator - level_end = mg_dof_handler.end_active(level); - // Compute coarse level right hand side - // by restricting from fine level. - for (; level_cell!=level_end; ++level_cell) - { - // get the dof numbers of - // this cell for the global - // and the level-wise - // numbering - level_cell->get_dof_indices(global_dof_indices); - level_cell->get_mg_dof_indices (level_dof_indices); - - // transfer the global - // defect in the vector - // into the level-wise one - const unsigned int level_start - = 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) - { - dst[level](level_dof_indices[*i] - level_start) - = src(global_dof_indices[*i] - offset); - } - } - - // for that part of the level - // which is further refined: - // get the defect by - // restriction of the defect on - // one level higher - if (static_cast(level) < maxlevel) - { - restrict_and_add (level+1, dst[level], dst[level+1]); - } - }; -} - - - template template void diff --git a/deal.II/deal.II/source/multigrid/mg_base.cc b/deal.II/deal.II/source/multigrid/mg_base.cc index 77e6ea3a0e..90476e759b 100644 --- a/deal.II/deal.II/source/multigrid/mg_base.cc +++ b/deal.II/deal.II/source/multigrid/mg_base.cc @@ -20,260 +20,10 @@ #include #include -#include DEAL_II_NAMESPACE_OPEN -namespace internal -{ - namespace mg - { - template - void - reinit_vector (const dealii::MGDoFHandler &mg_dof, - std::vector , - 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, - std::vector target_component, - MGLevelObject > &v) - { - const unsigned int n_blocks = mg_dof.get_fe().n_blocks(); - if (target_component.size()==0) - { - target_component.resize(n_blocks); - for (unsigned int i=0;i > - ndofs(mg_dof.get_tria().n_levels(), - std::vector(n_blocks)); - MGTools::count_dofs_per_block (mg_dof, ndofs, target_component); - - for (unsigned int level=v.get_minlevel(); - level<=v.get_maxlevel();++level) - { - v[level].reinit(n_target_blocks); - for (unsigned int b=0; b - 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() {} @@ -297,67 +47,6 @@ MGCoarseGridBase::~MGCoarseGridBase() // Explicit instantiations //TODO: Use the template expander script for this -namespace internal -{ - namespace mg - { - template void reinit_vector ( - const dealii::MGDoFHandler&, - std::vector, - MGLevelObject >&); - template void reinit_vector ( - const dealii::MGDoFHandler&, - std::vector, - MGLevelObject >&); - - template void reinit_vector ( - const dealii::MGDoFHandler&, - std::vector, - MGLevelObject >&); - template void reinit_vector ( - const dealii::MGDoFHandler&, - std::vector, - 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 >; diff --git a/deal.II/deal.II/source/multigrid/mg_transfer_block.cc b/deal.II/deal.II/source/multigrid/mg_transfer_block.cc index 4ac5aaa0d7..9684b3668c 100644 --- a/deal.II/deal.II/source/multigrid/mg_transfer_block.cc +++ b/deal.II/deal.II/source/multigrid/mg_transfer_block.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors +// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -28,10 +28,178 @@ #include #include +#include #include DEAL_II_NAMESPACE_OPEN +namespace +{ + /** + * 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 + 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 +template +void +MGTransferBlockSelect::copy_to_mg ( + const MGDoFHandler &mg_dof_handler, + MGLevelObject > &dst, + const BlockVector &src) const +{ + 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. + bool first = true; + for (unsigned int level=mg_dof_handler.get_tria().n_levels();level != 0;) + { + --level; + for (IT i= copy_indices[selected_block][level].begin(); + i != copy_indices[selected_block][level].end();++i) + dst[level](i->second) = src.block(selected_block)(i->first); + if (!first) + restrict_and_add (level+1, dst[level], dst[level+1]); + first = false; + } +} + + + +template +template +void +MGTransferBlockSelect::copy_to_mg ( + const MGDoFHandler &mg_dof_handler, + MGLevelObject > &dst, + const Vector &src) const +{ + 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. + bool first = true; + for (unsigned int level=mg_dof_handler.get_tria().n_levels();level != 0;) + { + --level; + for (IT i= copy_indices[selected_block][level].begin(); + i != copy_indices[selected_block][level].end();++i) + dst[level](i->second) = src(i->first); + if (!first) + restrict_and_add (level+1, dst[level], dst[level+1]); + first = false; + } +} + + + +template +template +void +MGTransferBlock::copy_to_mg ( + const MGDoFHandler& mg_dof_handler, + MGLevelObject >& dst, + const BlockVector& src) const +{ + 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;) + { + --level; + for (unsigned int block=0;blocksecond) = src.block(block)(i->first); + if (!first) + restrict_and_add (level+1, dst[level], dst[level+1]); + first = false; + } +} + template @@ -43,10 +211,10 @@ void MGTransferBlockBase::build_matrices ( const unsigned int n_blocks = fe.n_blocks(); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_levels = mg_dof.get_tria().n_levels(); - + Assert (selected.size() == n_blocks, ExcDimensionMismatch(selected.size(), n_blocks)); - + // Compute the mapping between real // blocks and blocks used for // multigrid computations. @@ -57,11 +225,11 @@ void MGTransferBlockBase::build_matrices ( mg_block[i] = n_mg_blocks++; else mg_block[i] = numbers::invalid_unsigned_int; - + // Compute the lengths of all blocks sizes.resize(n_levels); MGTools::count_dofs_per_block(mg_dof, sizes); - + // Fill some index vectors // for later use. mg_block_start = sizes; @@ -76,7 +244,7 @@ void MGTransferBlockBase::build_matrices ( k += t; } } - + block_start.resize(n_blocks); DoFTools::count_dofs_per_block (static_cast&>(mg_dof), block_start); @@ -97,9 +265,9 @@ void MGTransferBlockBase::build_matrices ( for (unsigned int block=0;block dof_indices_parent (dofs_per_cell); std::vector dof_indices_child (dofs_per_cell); - + // for each level: first build the // sparsity pattern of the matrices // and then build the matrices @@ -159,7 +327,7 @@ void MGTransferBlockBase::build_matrices ( 0, false); prolongation_sparsities[level]->collect_sizes(); - + for (typename MGDoFHandler::cell_iterator cell=mg_dof.begin(level); cell != mg_dof.end(level); ++cell) if (cell->has_children()) @@ -175,7 +343,7 @@ void MGTransferBlockBase::build_matrices ( // this child const FullMatrix &prolongation = mg_dof.get_fe().get_prolongation_matrix (child, cell->refinement_case()); - + cell->child(child)->get_mg_dof_indices (dof_indices_child); // now tag the entries in the @@ -214,7 +382,7 @@ void MGTransferBlockBase::build_matrices ( // this child const FullMatrix &prolongation = mg_dof.get_fe().get_prolongation_matrix (child, cell->refinement_case()); - + cell->child(child)->get_mg_dof_indices (dof_indices_child); // now set the entries in the @@ -246,13 +414,13 @@ void MGTransferBlockSelect::build_matrices ( { const FiniteElement& fe = mg_dof.get_fe(); unsigned int n_blocks = mg_dof.get_fe().n_blocks(); - + selected_block = select; selected.resize(n_blocks, false); selected[select] = true; - + MGTransferBlockBase::build_matrices (dof, mg_dof); - + std::vector global_dof_indices (fe.dofs_per_cell); std::vector level_dof_indices (fe.dofs_per_cell); for (int level=dof.get_tria().n_levels()-1; level>=0; --level) @@ -261,7 +429,7 @@ void MGTransferBlockSelect::build_matrices ( level_cell = mg_dof.begin_active(level); const typename MGDoFHandler::active_cell_iterator level_end = mg_dof.end_active(level); - + // Compute coarse level right hand side // by restricting from fine level. for (; level_cell!=level_end; ++level_cell) @@ -273,7 +441,7 @@ void MGTransferBlockSelect::build_matrices ( // numbering global_cell.get_dof_indices(global_dof_indices); level_cell->get_mg_dof_indices (level_dof_indices); - + for (unsigned int i=0; i::build_matrices ( { const FiniteElement& fe = mg_dof.get_fe(); unsigned int n_blocks = mg_dof.get_fe().n_blocks(); - + if (sel.size() != 0) { Assert(sel.size() == n_blocks, ExcDimensionMismatch(sel.size(), n_blocks)); selected = sel; - } + } if (selected.size() == 0) selected = std::vector (n_blocks, true); MGTransferBlockBase::build_matrices (dof, mg_dof); - + std::vector global_dof_indices (fe.dofs_per_cell); std::vector level_dof_indices (fe.dofs_per_cell); for (int level=dof.get_tria().n_levels()-1; level>=0; --level) @@ -318,7 +486,7 @@ void MGTransferBlock::build_matrices ( level_cell = mg_dof.begin_active(level); const typename MGDoFHandler::active_cell_iterator level_end = mg_dof.end_active(level); - + // Compute coarse level right hand side // by restricting from fine level. for (; level_cell!=level_end; ++level_cell) @@ -330,7 +498,7 @@ void MGTransferBlock::build_matrices ( // numbering global_cell.get_dof_indices(global_dof_indices); level_cell->get_mg_dof_indices (level_dof_indices); - + for (unsigned int i=0; i #include +#include DEAL_II_NAMESPACE_OPEN +namespace +{ + /** + * 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 + 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 +template +void +MGTransferSelect::do_copy_to_mg ( + const MGDoFHandler& mg_dof_handler, + MGLevelObject >& dst, + const InVector& src, + 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()); + + 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 + // slow. + std::vector selected_indices; + selected_indices.reserve(dofs_per_cell); + for (unsigned int i=0; i=static_cast(minlevel); --level) + { + + typename MGDoFHandler::active_cell_iterator + level_cell = mg_dof_handler.begin_active(level); + const typename MGDoFHandler::active_cell_iterator + level_end = mg_dof_handler.end_active(level); + // Compute coarse level right hand side + // by restricting from fine level. + for (; level_cell!=level_end; ++level_cell) + { + // get the dof numbers of + // this cell for the global + // and the level-wise + // numbering + level_cell->get_dof_indices(global_dof_indices); + level_cell->get_mg_dof_indices (level_dof_indices); + + // transfer the global + // defect in the vector + // into the level-wise one + const unsigned int level_start + = 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) + { + dst[level](level_dof_indices[*i] - level_start) + = src(global_dof_indices[*i] - offset); + } + } + + // for that part of the level + // which is further refined: + // get the defect by + // restriction of the defect on + // one level higher + if (static_cast(level) < maxlevel) + { + restrict_and_add (level+1, dst[level], dst[level+1]); + } + } +} + template void MGTransferComponentBase::build_matrices ( @@ -51,7 +331,7 @@ void MGTransferComponentBase::build_matrices ( Assert (target_component.size() == mg_dof.get_fe().n_components(), ExcDimensionMismatch(target_component.size(), mg_dof.get_fe().n_components())); - + for (unsigned int i=0;i& fe = mg_dof.get_fe(); // Effective number of components @@ -89,10 +369,10 @@ void MGTransferComponentBase::build_matrices ( *std::max_element(mg_target_component.begin(), mg_target_component.end()) + 1; const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_levels = mg_dof.get_tria().n_levels(); - + Assert (mg_selected.size() == fe.n_components(), ExcDimensionMismatch(mg_selected.size(), fe.n_components())); - + // Compute the lengths of all blocks sizes.resize(n_levels); MGTools::count_dofs_per_component(mg_dof, sizes, true, mg_target_component); @@ -134,9 +414,9 @@ void MGTransferComponentBase::build_matrices ( // too slow copy_to_and_from_indices.resize(n_levels); - + // Building the prolongation matrices starts here! - + // reset the size of the array of // matrices. call resize(0) first, // in order to delete all elements @@ -163,7 +443,7 @@ void MGTransferComponentBase::build_matrices ( // for a cell and one of its children std::vector dof_indices_parent (dofs_per_cell); std::vector dof_indices_child (dofs_per_cell); - + // for each level: first build the // sparsity pattern of the matrices // and then build the matrices @@ -196,7 +476,7 @@ void MGTransferComponentBase::build_matrices ( 0, false); prolongation_sparsities[level]->collect_sizes(); - + for (typename MGDoFHandler::cell_iterator cell=mg_dof.begin(level); cell != mg_dof.end(level); ++cell) if (cell->has_children()) @@ -210,7 +490,7 @@ void MGTransferComponentBase::build_matrices ( // this child const FullMatrix &prolongation = mg_dof.get_fe().get_prolongation_matrix (child, cell->refinement_case()); - + cell->child(child)->get_mg_dof_indices (dof_indices_child); // now tag the entries in the @@ -247,7 +527,7 @@ void MGTransferComponentBase::build_matrices ( // this child const FullMatrix &prolongation = mg_dof.get_fe().get_prolongation_matrix (child, cell->refinement_case()); - + cell->child(child)->get_mg_dof_indices (dof_indices_child); // now set the entries in the @@ -281,10 +561,10 @@ void MGTransferSelect::build_matrices ( { const FiniteElement& fe = mg_dof.get_fe(); unsigned int ncomp = mg_dof.get_fe().n_components(); - + target_component = t_component; mg_target_component = mg_t_component; - + selected_component = select; mg_selected_component = mg_select; selected.resize(ncomp, false); @@ -303,7 +583,7 @@ void MGTransferSelect::build_matrices ( break; } } - + for (unsigned int i=0;i::build_matrices ( mg_selected_component = i; break; } - } + } MGTransferComponentBase::build_matrices (dof, mg_dof); - + std::vector global_dof_indices (fe.dofs_per_cell); std::vector level_dof_indices (fe.dofs_per_cell); for (int level=dof.get_tria().n_levels()-1; level>=0; --level) @@ -322,7 +602,7 @@ void MGTransferSelect::build_matrices ( level_cell = mg_dof.begin_active(level); const typename MGDoFHandler::active_cell_iterator level_end = mg_dof.end_active(level); - + // Compute coarse level right hand side // by restricting from fine level. for (; level_cell!=level_end; ++level_cell) @@ -334,7 +614,7 @@ void MGTransferSelect::build_matrices ( // numbering global_cell.get_dof_indices(global_dof_indices); level_cell->get_mg_dof_indices (level_dof_indices); - + for (unsigned int i=0; i + void + reinit_vector (const dealii::MGDoFHandler &mg_dof, + std::vector , + 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); + } + + } + + + /** + * Adjust 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. The target_component + * is handed to MGTools::count_dofs_per_block. + * See for documentation there. + */ + template + void + reinit_vector (const dealii::MGDoFHandler &mg_dof, + std::vector target_component, + MGLevelObject > &v) + { + const unsigned int n_blocks = mg_dof.get_fe().n_blocks(); + if (target_component.size()==0) + { + target_component.resize(n_blocks); + for (unsigned int i=0;i > + ndofs(mg_dof.get_tria().n_levels(), + std::vector(n_blocks)); + MGTools::count_dofs_per_block (mg_dof, ndofs, target_component); + + for (unsigned int level=v.get_minlevel(); + level<=v.get_maxlevel();++level) + { + v[level].reinit(n_target_blocks); + for (unsigned int b=0; b +template +void +MGTransferPrebuilt::copy_to_mg ( + const MGDoFHandler& mg_dof_handler, + MGLevelObject& dst, + const InVector& src) const +{ + reinit_vector(mg_dof_handler, component_to_block_map, dst); + bool first = true; + for (unsigned int level=mg_dof_handler.get_tria().n_levels();level != 0;) + { + --level; + + typedef std::map::const_iterator IT; + for (IT i= copy_indices[level].begin(); + i != copy_indices[level].end();++i) + dst[level](i->second) = src(i->first); + // For non-DG: degrees of + // freedom in the refinement + // face may need special + // attention, since they belong + // to the coarse level, but + // have fine level basis + // functions + if (!first) + restrict_and_add (level+1, dst[level], dst[level+1]); + first = false; + } +} + template @@ -233,6 +342,7 @@ void MGTransferPrebuilt::build_matrices ( } +//TODO: Use template expander script template void MGTransferPrebuilt >::build_matrices -- 2.39.5