From 8ef8ae2d2cd35f61ab37cf5ab33947c2e88f6a4d Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 1 Feb 2010 15:55:59 +0000 Subject: [PATCH] Use std::vector > instead of std::map for describing the transfer from and to multigrid. This is considerably faster. Also improve imposition of boundary conditions in mg_transfer_prebuilt git-svn-id: https://svn.dealii.org/trunk@20485 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/multigrid/mg_transfer.h | 4 +- .../include/multigrid/mg_transfer.templates.h | 3 +- .../include/multigrid/mg_transfer_block.h | 4 +- .../multigrid/mg_transfer_block.templates.h | 2 +- .../source/multigrid/mg_transfer_block.cc | 41 ++++- .../source/multigrid/mg_transfer_prebuilt.cc | 161 +++++++++++------- 6 files changed, 143 insertions(+), 72 deletions(-) diff --git a/deal.II/deal.II/include/multigrid/mg_transfer.h b/deal.II/deal.II/include/multigrid/mg_transfer.h index 630cd295f7..670ed5f0b3 100644 --- a/deal.II/deal.II/include/multigrid/mg_transfer.h +++ b/deal.II/deal.II/include/multigrid/mg_transfer.h @@ -220,8 +220,8 @@ class MGTransferPrebuilt : public MGTransferBase * The data is first the global * index, then the level index. */ - std::vector > - copy_indices; + std::vector > > + copy_indices; /** * The vector that stores what 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 6ddff57548..96cdc1c4f7 100644 --- a/deal.II/deal.II/include/multigrid/mg_transfer.templates.h +++ b/deal.II/deal.II/include/multigrid/mg_transfer.templates.h @@ -28,6 +28,7 @@ DEAL_II_NAMESPACE_OPEN /* --------------------- MGTransferPrebuilt -------------- */ +typedef std::vector >::const_iterator IT; @@ -49,7 +50,6 @@ MGTransferPrebuilt::copy_from_mg( dst = 0; for (unsigned int level=0;level::const_iterator IT; for (IT i= copy_indices[level].begin(); i != copy_indices[level].end();++i) dst(i->first) = src[level](i->second); @@ -75,7 +75,6 @@ MGTransferPrebuilt::copy_from_mg_add ( // to the coarse level, but // have fine level basis // functions - typedef std::map::const_iterator IT; for (unsigned int level=0;level > > - copy_indices; + std::vector > > > + copy_indices; }; /** 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 b80d4d6528..95028b8220 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 @@ -30,7 +30,7 @@ DEAL_II_NAMESPACE_OPEN /* --------------------- MGTransferBlockSelect -------------- */ // Simplify some things below -typedef std::map::const_iterator IT; +typedef std::vector >::const_iterator IT; 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 c2a149d0bb..119432c140 100644 --- a/deal.II/deal.II/source/multigrid/mg_transfer_block.cc +++ b/deal.II/deal.II/source/multigrid/mg_transfer_block.cc @@ -447,10 +447,29 @@ void MGTransferBlockSelect::build_matrices ( { const unsigned int block = fe.system_to_block_index(i).first; if (selected[block]) - copy_indices[block][level].insert( + copy_indices[block][level].push_back( std::make_pair(global_dof_indices[i] - block_start[block], level_dof_indices[i] - mg_block_start[level][block])); } + } + // sort the list of dof numbers and compress + // out duplicates + for (unsigned int block=0; block 0) + { + std::sort (copy_indices[block][level].begin(), + copy_indices[block][level].end()); + + std::vector >::iterator + it1 = copy_indices[block][level].begin(), it2 = it1, it3 = it1; + it2++, it1++; + for ( ; it2!=copy_indices[block][level].end(); ++it2, ++it3) + if (*it2 > *it3) + *it1++ = *it2; + copy_indices[block][level].erase(it1,copy_indices[block][level].end()); + } } } } @@ -504,11 +523,29 @@ void MGTransferBlock::build_matrices ( { const unsigned int block = fe.system_to_block_index(i).first; if (selected[block]) - copy_indices[block][level].insert( + copy_indices[block][level].push_back( std::make_pair(global_dof_indices[i] - block_start[block], level_dof_indices[i] - mg_block_start[level][block])); } } + + for (unsigned int block=0; block 0) + { + std::sort (copy_indices[block][level].begin(), + copy_indices[block][level].end()); + + std::vector >::iterator + it1 = copy_indices[block][level].begin(), it2 = it1, it3 = it1; + it2++, it1++; + for ( ; it2!=copy_indices[block][level].end(); ++it2, ++it3) + if (*it2 > *it3) + *it1++ = *it2; + copy_indices[block][level].erase(it1,copy_indices[block][level].end()); + } + } } } diff --git a/deal.II/deal.II/source/multigrid/mg_transfer_prebuilt.cc b/deal.II/deal.II/source/multigrid/mg_transfer_prebuilt.cc index 51263fcf10..8e3972bf72 100644 --- a/deal.II/deal.II/source/multigrid/mg_transfer_prebuilt.cc +++ b/deal.II/deal.II/source/multigrid/mg_transfer_prebuilt.cc @@ -122,7 +122,7 @@ MGTransferPrebuilt::copy_to_mg ( { --level; - typedef std::map::const_iterator IT; + typedef std::vector >::const_iterator IT; for (IT i= copy_indices[level].begin(); i != copy_indices[level].end();++i) dst[level](i->second) = src(i->first); @@ -181,6 +181,9 @@ void MGTransferPrebuilt::build_matrices ( std::vector dof_indices_parent (dofs_per_cell); std::vector dof_indices_child (dofs_per_cell); + copy_indices.resize(mg_dof.get_tria().n_levels()); + find_dofs_on_refinement_edges (mg_dof); + // for each level: first build the sparsity // pattern of the matrices and then build the // matrices themselves. note that we only @@ -238,6 +241,11 @@ void MGTransferPrebuilt::build_matrices ( prolongation_matrices[level]->reinit (*prolongation_sparsities[level]); + copy_indices[level+1].clear(); + if (level == 0) + copy_indices[0].clear(); + std::vector global_dof_indices (dofs_per_cell); + // now actually build the matrices for (typename MGDoFHandler::cell_iterator cell=mg_dof.begin(level); cell != mg_dof.end(level); ++cell) @@ -260,82 +268,109 @@ void MGTransferPrebuilt::build_matrices ( // now set the entries in the // matrix for (unsigned int i=0; iset (dof_indices_child[i], - dof_indices_parent[j], - prolongation(i,j)); + prolongation_matrices[level]->set (dof_indices_child[i], + dofs_per_cell, + &dof_indices_parent[0], + &prolongation(i,0), + true); + + if (cell->child(child)->has_children()==false) + { + DoFCellAccessor >& global_cell + = *cell->child(child); + global_cell.get_dof_indices(global_dof_indices); + for (unsigned int i=0; i >& global_cell = *cell; + + // get the dof numbers of + // this cell for the global + // and the level-wise + // numbering + global_cell.get_dof_indices(global_dof_indices); + cell->get_mg_dof_indices (dof_indices_parent); + + for (unsigned int i=0; i=0; --level) + { + if (copy_indices[level].size() > 0) + { + std::sort (copy_indices[level].begin(), copy_indices[level].end()); + + std::vector >::iterator + it1 = copy_indices[level].begin(), it2 = it1, it3 = it1; + it2++, it1++; + for ( ; it2!=copy_indices[level].end(); ++it2, ++it3) + if (*it2 > *it3) + *it1++ = *it2; + copy_indices[level].erase(it1,copy_indices[level].end()); + } + } + // if we have only one single level, we did + // not enter the above code. Need to create + // the respective vector manually + if (n_levels == 1) + { + copy_indices[0].resize (sizes[0]); + for (unsigned int i=0; i(i,i); } + // impose boundary conditions // but only in the column of // the prolongation matrix - //TODO: this way is not very efficient - std::vector > boundary_indices (mg_dof.get_tria().n_levels()); + std::vector > boundary_indices (n_levels); typename FunctionMap::type boundary; ZeroFunction homogeneous_dirichlet_bc (1); boundary[0] = &homogeneous_dirichlet_bc; MGTools::make_boundary_list(mg_dof, boundary, boundary_indices); - for (unsigned int level=0; levelm(); - for (unsigned int i=0; i constrain_indices; + for (int level=n_levels-2; level>=0; --level) { - SparseMatrix::iterator anfang = prolongation_matrices[level]->begin(i), - ende = prolongation_matrices[level]->end(i); - for(; anfang != ende; ++anfang) - { - std::set::const_iterator dof = boundary_indices[level].begin(), - endd = boundary_indices[level].end(); - for (; dof != endd; ++dof) - { - const unsigned int column_number = anfang->column(); - const unsigned int dof_number = *dof; - if(dof_number == column_number) - prolongation_matrices[level]->set(i,dof_number,0); - } - } - } - } - - find_dofs_on_refinement_edges (mg_dof); - // Prepare indices to copy from a - // global vector to the parts of an - // MGVector corresponding to active - // cells. - copy_indices.resize(mg_dof.get_tria().n_levels()); - std::vector global_dof_indices (dofs_per_cell); - std::vector level_dof_indices (dofs_per_cell); - for (int level=mg_dof.get_tria().n_levels()-1; level>=0; --level) - { - copy_indices[level].clear(); - typename MGDoFHandler::active_cell_iterator - 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) + if (boundary_indices[level].size() == 0) + continue; + + // need to delete all the columns in the + // matrix that are on the boundary. to fix + // this, create an array as long as there are + // matrix columns, and find which columns we + // need to filter away. + constrain_indices.resize (0); + constrain_indices.resize (prolongation_matrices[level]->n(), 0); + std::set::const_iterator dof = boundary_indices[level].begin(), + endd = boundary_indices[level].end(); + for (; dof != endd; ++dof) + constrain_indices[*dof] = 1; + + const unsigned int n_dofs = prolongation_matrices[level]->m(); + for (unsigned int i=0; i >& global_cell = *level_cell; - // get the dof numbers of - // this cell for the global - // and the level-wise - // numbering - global_cell.get_dof_indices(global_dof_indices); - level_cell->get_mg_dof_indices (level_dof_indices); - - for (unsigned int i=0; i::iterator start_row = prolongation_matrices[level]->begin(i), + end_row = prolongation_matrices[level]->end(i); + for(; start_row != end_row; ++start_row) { - if(!dofs_on_refinement_edge[level][level_dof_indices[i]]) - copy_indices[level].insert( - std::make_pair(global_dof_indices[i], level_dof_indices[i])); + if (constrain_indices[start_row->column()] == 1) + start_row->value() = 0; } } } -- 2.39.5