From: Baerbel Jannsen Date: Fri, 5 Feb 2010 13:56:31 +0000 (+0000) Subject: prepared MGTransferSelect for adaptive refinement, cleaned up mg_dof_tools.cc X-Git-Tag: v8.0.0~6541 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c3b888fc307a02881591c5ea5add58c313924089;p=dealii.git prepared MGTransferSelect for adaptive refinement, cleaned up mg_dof_tools.cc git-svn-id: https://svn.dealii.org/trunk@20506 0785d39b-7218-0410-832d-ea1e28bc413d --- 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 46c2d25134..aa11f8b5d6 100644 --- a/deal.II/deal.II/source/multigrid/mg_dof_tools.cc +++ b/deal.II/deal.II/source/multigrid/mg_dof_tools.cc @@ -30,7 +30,6 @@ #include #include -#include #include #include #include @@ -1072,7 +1071,10 @@ MGTools::count_dofs_per_block ( std::vector > dofs_in_block (n_blocks, std::vector(dof_handler.n_dofs(l), false)); std::vector > - block_select (n_blocks, std::vector(n_blocks, false)); + block_select; + block_select.resize(n_blocks); + for (int iii=0;iii tasks; for (unsigned int i=0; i +void +MGTools:: +extract_inner_interface_dofs (const MGDoFHandler &mg_dof_handler, + std::vector > &interface_dofs) +{ + Assert (interface_dofs.size() == mg_dof_handler.get_tria().n_levels(), + ExcDimensionMismatch (interface_dofs.size(), + mg_dof_handler.get_tria().n_levels())); + + for (unsigned int l=0; l &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; + + std::vector local_dof_indices (dofs_per_cell); + std::vector cell_dofs(dofs_per_cell, false); + + typename MGDoFHandler::cell_iterator cell = mg_dof_handler.begin(), + endc = mg_dof_handler.end(); + + for (; cell!=endc; ++cell) + { + std::fill (cell_dofs.begin(), 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); + if (!face->at_boundary()) + { + //interior face + const typename MGDoFHandler::cell_iterator + neighbor = cell->neighbor(face_nr); + + if (neighbor->level() < cell->level()) + { + for (unsigned int j=0; jlevel(); + cell->get_mg_dof_indices (local_dof_indices); + + for(unsigned int i=0; i @@ -1380,7 +1444,7 @@ extract_inner_interface_dofs (const MGDoFHandler &mg_dof_handler, for (unsigned int face_nr=0; face_nr::faces_per_cell; ++face_nr) if(cell->at_boundary(face_nr)) for(unsigned int j=0; j &mg_dof_handler, std::vector > &interface_dofs, std::vector > &boundary_interface_dofs); +template +void +MGTools:: +extract_inner_interface_dofs (const MGDoFHandler &mg_dof_handler, + std::vector > &interface_dofs); #endif DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/multigrid/mg_transfer_component.cc b/deal.II/deal.II/source/multigrid/mg_transfer_component.cc index 7d2ae176eb..26d2891c86 100644 --- a/deal.II/deal.II/source/multigrid/mg_transfer_component.cc +++ b/deal.II/deal.II/source/multigrid/mg_transfer_component.cc @@ -12,9 +12,11 @@ //--------------------------------------------------------------------------- #include +#include #include #include +#include #include #include #include @@ -29,7 +31,7 @@ #include #include - +#include DEAL_II_NAMESPACE_OPEN @@ -123,8 +125,7 @@ namespace 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); + MGTools::count_dofs_per_block (mg_dof, ndofs, target_component); } for (unsigned int level=v.get_minlevel(); @@ -184,25 +185,26 @@ namespace // 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)); +// 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; - + 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); + MGTools::count_dofs_per_block (mg_dof, ndofs, + target_component); } for (unsigned int level=v.get_minlevel(); @@ -220,18 +222,8 @@ void MGTransferSelect::do_copy_to_mg ( const MGDoFHandler& mg_dof_handler, MGLevelObject >& dst, - const InVector& src, - const unsigned int offset) const + const InVector& src) 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(), @@ -240,21 +232,6 @@ MGTransferSelect::do_copy_to_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::do_copy_to_mg ( // the respective vector on the // next finer level, which we then // already have built. - for (int level=maxlevel; level>=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); - } - } + 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_to_and_from_indices[level].begin(); + i != copy_to_and_from_indices[level].end(); ++i) + dst[level](i->second) = src(i->first); // 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]); - } + if(!first) + restrict_and_add (level+1, dst[level], dst[level+1]); + first = false; } } @@ -375,7 +326,10 @@ void MGTransferComponentBase::build_matrices ( // Compute the lengths of all blocks sizes.resize(n_levels); - MGTools::count_dofs_per_component(mg_dof, sizes, true, mg_target_component); + for(unsigned int l=0; l&>(mg_dof), - component_start, true, target_component); + count_dofs_per_block (mg_dof, component_start, target_component); unsigned int k=0; for (unsigned int i=0;iadd(dof_indices_child[i], dof_indices_parent[j]); }; @@ -538,7 +491,7 @@ void MGTransferComponentBase::build_matrices ( { const unsigned int icomp = fe.system_to_component_index(i).first; const unsigned int jcomp = fe.system_to_component_index(j).first; - if ((icomp==jcomp) && mg_selected[mg_target_component[icomp]]) + if ((icomp==jcomp) && mg_selected[icomp]) prolongation_matrices[level]->set(dof_indices_child[i], dof_indices_parent[j], prolongation(i,j)); @@ -546,6 +499,56 @@ void MGTransferComponentBase::build_matrices ( } } } + // impose boundary conditions + // but only in the column of + // the prolongation matrix + //TODO: this way is not very efficient + + if(boundary_indices.size() != 0) + { + std::vector > + dofs_per_component(mg_dof.get_tria().n_levels(), + std::vector(n_components)); + + MGTools::count_dofs_per_block (mg_dof, dofs_per_component, mg_target_component); + for (unsigned int level=0; levelblock(iblock,jblock).m(); + for (unsigned int i=0; i::iterator anfang = prolongation_matrices[level]->block(iblock,jblock).begin(i), + ende = prolongation_matrices[level]->block(iblock,jblock).end(i); + for(; anfang != ende; ++anfang) + { + const unsigned int column_number = anfang->column(); + + //convert global indices into local ones + const BlockIndices block_indices_coarse (dofs_per_component[level]); + const unsigned int global_j = block_indices_coarse.local_to_global(iblock, column_number); + + std::set::const_iterator found_dof = + boundary_indices[level].find(global_j); + + const bool is_boundary_index = + (found_dof != boundary_indices[level].end()); + + if(is_boundary_index) + { + prolongation_matrices[level]->block(iblock,jblock) + .set(i,column_number,0); + } + } + } + } + } + } } @@ -557,20 +560,27 @@ void MGTransferSelect::build_matrices ( unsigned int select, unsigned int mg_select, const std::vector& t_component, - const std::vector& mg_t_component) + const std::vector& mg_t_component, + const std::vector >& bdry_indices) { 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; + boundary_indices = bdry_indices; selected_component = select; mg_selected_component = mg_select; selected.resize(ncomp, false); - selected[select] = true; + for(unsigned int c=0; c::build_matrices ( break; } } + MGTransferComponentBase::build_matrices (dof, mg_dof); + interface_dofs.resize(mg_dof.get_tria().n_levels()); + for(unsigned int l=0; l 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) { + copy_to_and_from_indices[level].clear(); typename MGDoFHandler::active_cell_iterator level_cell = mg_dof.begin_active(level); const typename MGDoFHandler::active_cell_iterator @@ -618,14 +635,15 @@ void MGTransferSelect::build_matrices ( for (unsigned int i=0; i::build_matrices const MGDoFHandler &, unsigned int, unsigned int, const std::vector&, - const std::vector&); + const std::vector&, + const std::vector >&); template void MGTransferSelect::build_matrices @@ -651,7 +670,8 @@ void MGTransferSelect::build_matrices const MGDoFHandler &, unsigned int, unsigned int, const std::vector&, - const std::vector&); + const std::vector&, + const std::vector >&); template void MGTransferSelect::copy_to_mg ( 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 0e676f06c7..fa7bcc0bab 100644 --- a/deal.II/deal.II/source/multigrid/mg_transfer_prebuilt.cc +++ b/deal.II/deal.II/source/multigrid/mg_transfer_prebuilt.cc @@ -123,7 +123,7 @@ MGTransferPrebuilt::copy_to_mg ( --level; VECTOR& dst_level = dst[level]; - typedef std::vector >::const_iterator IT; + 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); @@ -145,7 +145,9 @@ MGTransferPrebuilt::copy_to_mg ( template template void MGTransferPrebuilt::build_matrices ( - const MGDoFHandler &mg_dof) + const MGDoFHandler &mg_dof, + const std::vector >&boundary_indices + ) { //start building the matrices here const unsigned int n_levels = mg_dof.get_tria().n_levels(); @@ -274,42 +276,48 @@ void MGTransferPrebuilt::build_matrices ( // impose boundary conditions // but only in the column of // the prolongation matrix - 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); - - std::vector constrain_indices; - for (int level=n_levels-2; level>=0; --level) + + if(boundary_indices.size() != 0) + { + std::vector constrain_indices; + for (int level=n_levels-2; level>=0; --level) +// for (unsigned int level=0; leveln(), 0); std::set::const_iterator dof = boundary_indices[level].begin(), - endd = boundary_indices[level].end(); + endd = boundary_indices[level].end(); for (; dof != endd; ++dof) - constrain_indices[*dof] = 1; + constrain_indices[*dof] = 1; const unsigned int n_dofs = prolongation_matrices[level]->m(); 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 (constrain_indices[start_row->column()] == 1) - start_row->value() = 0; - } - } + { + SparseMatrix::iterator start_row = prolongation_matrices[level]->begin(i), + end_row = prolongation_matrices[level]->end(i); + for(; start_row != end_row; ++start_row) + { + std::set::const_iterator dof = boundary_indices[level].begin(), + endd = boundary_indices[level].end(); + for (; dof != endd; ++dof) + { + const unsigned int column_number = start_row->column(); + const unsigned int dof_number = *dof; + if(dof_number == column_number) + prolongation_matrices[level]->set(i,dof_number,0); + } + } + } } + } // to find the indices that describe the // relation between global dofs and local @@ -321,11 +329,13 @@ void MGTransferPrebuilt::build_matrices ( // the std::vector > that // only contains the active dofs on the // levels. - - find_dofs_on_refinement_edges (mg_dof); + interface_dofs.resize(mg_dof.get_tria().n_levels()); + for(unsigned int l=0; l temp_copy_indices; +// std::vector temp_copy_indices; 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) @@ -336,8 +346,8 @@ void MGTransferPrebuilt::build_matrices ( const typename MGDoFHandler::active_cell_iterator level_end = mg_dof.end_active(level); - temp_copy_indices.resize (0); - temp_copy_indices.resize (mg_dof.n_dofs(level), numbers::invalid_unsigned_int); +// temp_copy_indices.resize (0); +// temp_copy_indices.resize (mg_dof.n_dofs(level), numbers::invalid_unsigned_int); // Compute coarse level right hand side // by restricting from fine level. @@ -352,8 +362,14 @@ void MGTransferPrebuilt::build_matrices ( level_cell->get_mg_dof_indices (level_dof_indices); for (unsigned int i=0; i::build_matrices ( // copy_indices object. Then, insert the pairs // of global index and level index into // copy_indices. - const unsigned int n_active_dofs = - std::count_if (temp_copy_indices.begin(), temp_copy_indices.end(), - std::bind2nd(std::not_equal_to(), - numbers::invalid_unsigned_int)); - copy_indices[level].resize (n_active_dofs); - unsigned int counter = 0; - for (unsigned int i=0; i (temp_copy_indices[i], i); - Assert (counter == n_active_dofs, ExcInternalError()); +// const unsigned int n_active_dofs = +// std::count_if (temp_copy_indices.begin(), temp_copy_indices.end(), +// std::bind2nd(std::not_equal_to(), +// numbers::invalid_unsigned_int)); +// copy_indices[level].resize (n_active_dofs); +// unsigned int counter = 0; +// for (unsigned int i=0; i (temp_copy_indices[i], i); +// Assert (counter == n_active_dofs, ExcInternalError()); } } @@ -381,19 +397,23 @@ void MGTransferPrebuilt::build_matrices ( template void MGTransferPrebuilt >::build_matrices -(const MGDoFHandler &mg_dof); +(const MGDoFHandler &mg_dof, + const std::vector >&boundary_indices); template void MGTransferPrebuilt >::build_matrices -(const MGDoFHandler &mg_dof); +(const MGDoFHandler &mg_dof, + const std::vector >&boundary_indices); template void MGTransferPrebuilt >::build_matrices -(const MGDoFHandler &mg_dof); +(const MGDoFHandler &mg_dof, + const std::vector >&boundary_indices); template void MGTransferPrebuilt >::build_matrices -(const MGDoFHandler &mg_dof); +(const MGDoFHandler &mg_dof, + const std::vector >&boundary_indices); template void MGTransferPrebuilt >::copy_to_mg (