From f49d41e35dc0ed2771f8608e48ae75f3b63007cc Mon Sep 17 00:00:00 2001 From: Karl Ljungkvist Date: Thu, 17 Nov 2016 17:26:57 +0100 Subject: [PATCH] Extract some common MG transfer setup code --- .../deal.II/multigrid/mg_transfer_internal.h | 74 ++- source/multigrid/mg_transfer_internal.cc | 473 +++++++++++++++++- source/multigrid/mg_transfer_internal.inst.in | 52 ++ source/multigrid/mg_transfer_matrix_free.cc | 464 +---------------- 4 files changed, 613 insertions(+), 450 deletions(-) diff --git a/include/deal.II/multigrid/mg_transfer_internal.h b/include/deal.II/multigrid/mg_transfer_internal.h index 27532a3ed6..dec3ea3b11 100644 --- a/include/deal.II/multigrid/mg_transfer_internal.h +++ b/include/deal.II/multigrid/mg_transfer_internal.h @@ -17,10 +17,12 @@ #ifndef dealii__mg_transfer_internal_h #define dealii__mg_transfer_internal_h +#include #include +#include +#include #include - DEAL_II_NAMESPACE_OPEN namespace internal @@ -38,6 +40,76 @@ namespace internal std::vector > > ©_indices, std::vector > > ©_indices_global_mine, std::vector > > ©_indices_level_mine); + + + + /** + * Given the collection of child cells in lexicographic ordering as seen + * from the parent, this function computes the first index of the given + * child + */ + template + unsigned int + compute_shift_within_children(const unsigned int child, + const unsigned int fe_shift_1d, + const unsigned int fe_degree); + + /** + * Stores data related to the finite element contained in the + * DoFHandler. Used only for the initialization using + * setup_transfer. + */ + template + struct ElementInfo + { + + /** + * Stores the degree of the finite element. The selection of the + * computational kernel is based on this number. + */ + unsigned int fe_degree; + + /** + * Stores whether the element is continuous and there is a joint degree of + * freedom in the center of the 1D line. + */ + bool element_is_continuous; + + /** + * Stores the number of components in the finite element. + */ + unsigned int n_components; + + /** + * Stores the number of degrees of freedom on all child cells. It is + * 2dim*fe.dofs_per_cell for DG elements and somewhat + * less for continuous elements. + */ + unsigned int n_child_cell_dofs; + + /** + * Holds the one-dimensional embedding (prolongation) matrix from mother + * element to the children. + */ + internal::MatrixFreeFunctions::ShapeInfo shape_info; + + }; + + /** + * Sets up most of the internal data structures of MGTransferMatrixFree + */ + template + void setup_transfer(const dealii::DoFHandler &mg_dof, + const MGConstrainedDoFs *mg_constrained_dofs, + ElementInfo &elem_info, + std::vector > &level_dof_indices, + std::vector > > &parent_child_connect, + std::vector &n_owned_level_cells, + std::vector > > &dirichlet_indices, + std::vector > > &weights_on_refined, + std::vector > > ©_indices_global_mine, + MGLevelObject > &ghosted_level_vector); + } } diff --git a/source/multigrid/mg_transfer_internal.cc b/source/multigrid/mg_transfer_internal.cc index 1c3696e078..6561d033a0 100644 --- a/source/multigrid/mg_transfer_internal.cc +++ b/source/multigrid/mg_transfer_internal.cc @@ -16,6 +16,7 @@ #include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -24,11 +25,10 @@ namespace internal { namespace MGTransfer { - /** - * Internal data structure that is used in the MPI communication in - * fill_copy_indices(). It represents an entry in the copy_indices* map, - * that associates a level dof index with a global dof index. - */ + + // Internal data structure that is used in the MPI communication in + // fill_copy_indices(). It represents an entry in the copy_indices* map, + // that associates a level dof index with a global dof index. struct DoFPair { unsigned int level; @@ -52,6 +52,8 @@ namespace internal {} }; + // Internal function for filling the copy indices from global to level + // indices template void fill_copy_indices(const dealii::DoFHandler &mg_dof, const MGConstrainedDoFs *mg_constrained_dofs, @@ -272,6 +274,467 @@ namespace internal for (unsigned int level=0; level + void + reinit_ghosted_vector(const IndexSet &locally_owned, + std::vector &ghosted_level_dofs, + const MPI_Comm &communicator, + LinearAlgebra::distributed::Vector &ghosted_level_vector, + std::vector > ©_indices_global_mine) + { + std::sort(ghosted_level_dofs.begin(), ghosted_level_dofs.end()); + IndexSet ghosted_dofs(locally_owned.size()); + ghosted_dofs.add_indices(ghosted_level_dofs.begin(), + std::unique(ghosted_level_dofs.begin(), + ghosted_level_dofs.end())); + ghosted_dofs.compress(); + + // Add possible ghosts from the previous content in the vector + if (ghosted_level_vector.size() == locally_owned.size()) + { + // shift the local number of the copy indices according to the new + // partitioner that we are going to use for the vector + const std_cxx11::shared_ptr part + = ghosted_level_vector.get_partitioner(); + ghosted_dofs.add_indices(part->ghost_indices()); + for (unsigned int i=0; ilocal_to_global(copy_indices_global_mine[i].second)); + } + ghosted_level_vector.reinit(locally_owned, ghosted_dofs, communicator); + } + + // Transform the ghost indices to local index space for the vector + inline void + copy_indices_to_mpi_local_numbers(const Utilities::MPI::Partitioner &part, + const std::vector &mine, + const std::vector &remote, + std::vector &localized_indices) + { + localized_indices.resize(mine.size()+remote.size(), + numbers::invalid_unsigned_int); + for (unsigned int i=0; i + unsigned int + compute_shift_within_children(const unsigned int child, + const unsigned int fe_shift_1d, + const unsigned int fe_degree) + { + // we put the degrees of freedom of all child cells in lexicographic + // ordering + unsigned int c_tensor_index[dim]; + unsigned int tmp = child; + for (unsigned int d=0; d + void add_child_indices(const unsigned int child, + const unsigned int fe_shift_1d, + const unsigned int fe_degree, + const std::vector &lexicographic_numbering, + const std::vector &local_dof_indices, + types::global_dof_index *target_indices) + { + const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d; + const unsigned int shift = + compute_shift_within_children(child, fe_shift_1d, fe_degree); + const unsigned int n_components = + local_dof_indices.size()/Utilities::fixed_power(fe_degree+1); + types::global_dof_index *indices = target_indices + shift; + const unsigned int n_scalar_cell_dofs = Utilities::fixed_power(n_child_dofs_1d); + for (unsigned int c=0, m=0; c2 ? (fe_degree+1) : 1); ++k) + for (unsigned int j=0; j<(dim>1 ? (fe_degree+1) : 1); ++j) + for (unsigned int i=0; i<(fe_degree+1); ++i, ++m) + { + const unsigned int index = c*n_scalar_cell_dofs+k*n_child_dofs_1d* + n_child_dofs_1d+j*n_child_dofs_1d+i; + Assert(indices[index] == numbers::invalid_dof_index || + indices[index] == local_dof_indices[lexicographic_numbering[m]], + ExcInternalError()); + indices[index] = local_dof_indices[lexicographic_numbering[m]]; + } + } + + template + void setup_element_info(ElementInfo &elem_info,const FiniteElement<1> &fe, + const dealii::DoFHandler &mg_dof) + { + // currently, we have only FE_Q and FE_DGQ type elements implemented + elem_info.n_components = mg_dof.get_fe().element_multiplicity(0); + AssertDimension(Utilities::fixed_power(fe.dofs_per_cell)*elem_info.n_components, + mg_dof.get_fe().dofs_per_cell); + AssertDimension(fe.degree, mg_dof.get_fe().degree); + elem_info.fe_degree = fe.degree; + elem_info.element_is_continuous = fe.dofs_per_vertex > 0; + Assert(fe.dofs_per_vertex < 2, ExcNotImplemented()); + + // step 1.2: get renumbering of 1D basis functions to lexicographic + // numbers. The distinction according to fe.dofs_per_vertex is to support + // both continuous and discontinuous bases. + std::vector renumbering(fe.dofs_per_cell); + { + AssertIndexRange(fe.dofs_per_vertex, 2); + renumbering[0] = 0; + for (unsigned int i=0; i::vertices_per_cell*fe.dofs_per_vertex + i; + if (fe.dofs_per_vertex > 0) + renumbering[fe.dofs_per_cell-fe.dofs_per_vertex] = fe.dofs_per_vertex; + } + + // step 1.3: create a 1D quadrature formula from the finite element that + // collects the support points of the basis functions on the two children. + std::vector > basic_support_points = fe.get_unit_support_points(); + Assert(fe.dofs_per_vertex == 0 || fe.dofs_per_vertex == 1, + ExcNotImplemented()); + std::vector > points_refined(fe.dofs_per_vertex > 0 ? + (2 * fe.dofs_per_cell - 1) : + (2 * fe.dofs_per_cell)); + const unsigned int shift = fe.dofs_per_cell - fe.dofs_per_vertex; + for (unsigned int c=0; c::max_children_per_cell; ++c) + for (unsigned int j=0; j(n_child_dofs_1d); + + // step 1.4: evaluate the polynomials and store the data in ShapeInfo + const Quadrature<1> quadrature(points_refined); + elem_info.shape_info.reinit(quadrature, mg_dof.get_fe(), 0); + + for (unsigned int c=0; c::max_children_per_cell; ++c) + for (unsigned int i=0; i::epsilon(),1e-12), + ExcInternalError()); + } + + + // Sets up most of the internal data structures of the MGTransferMatrixFree + // class + template + void setup_transfer(const dealii::DoFHandler &mg_dof, + const MGConstrainedDoFs *mg_constrained_dofs, + ElementInfo &elem_info, + std::vector > &level_dof_indices, + std::vector > > &parent_child_connect, + std::vector &n_owned_level_cells, + std::vector > > &dirichlet_indices, + std::vector > > &weights_on_refined, + std::vector > > ©_indices_global_mine, + MGLevelObject > &ghosted_level_vector) + { + + // we collect all child DoFs of a mother cell together. For faster + // tensorized operations, we align the degrees of freedom + // lexicographically. We distinguish FE_Q elements and FE_DGQ elements + + const dealii::Triangulation &tria = mg_dof.get_triangulation(); + + // ---------------------------- 1. Extract 1D info about the finite element + // step 1.1: create a 1D copy of the finite element from FETools where we + // substitute the template argument + AssertDimension(mg_dof.get_fe().n_base_elements(), 1); + std::string fe_name = mg_dof.get_fe().base_element(0).get_name(); + { + const std::size_t template_starts = fe_name.find_first_of('<'); + Assert (fe_name[template_starts+1] == (dim==1?'1':(dim==2?'2':'3')), + ExcInternalError()); + fe_name[template_starts+1] = '1'; + } + std_cxx11::shared_ptr > fe_1d + (FETools::get_fe_by_name<1,1>(fe_name)); + const FiniteElement<1> &fe = *fe_1d; + + setup_element_info(elem_info,fe,mg_dof); + + unsigned int n_child_dofs_1d = elem_info.shape_info.n_q_points; + + // -------------- 2. Extract and match dof indices between child and parent + const unsigned int n_levels = tria.n_global_levels(); + level_dof_indices.resize(n_levels); + parent_child_connect.resize(n_levels-1); + n_owned_level_cells.resize(n_levels-1); + std::vector > coarse_level_indices(n_levels-1); + for (unsigned int level=0; level local_dof_indices(mg_dof.get_fe().dofs_per_cell); + dirichlet_indices.resize(n_levels-1); + + // We use the vectors stored ghosted_level_vector in the base class for + // keeping ghosted transfer indices. To avoid keeping two very similar + // vectors, we merge them here. + if (ghosted_level_vector.max_level() != n_levels-1) + ghosted_level_vector.resize(0, n_levels-1); + + for (unsigned int level=n_levels-1; level > 0; --level) + { + unsigned int counter = 0; + std::vector global_level_dof_indices; + std::vector global_level_dof_indices_remote; + std::vector ghosted_level_dofs; + std::vector global_level_dof_indices_l0; + std::vector ghosted_level_dofs_l0; + + // step 2.1: loop over the cells on the coarse side + typename dealii::DoFHandler::cell_iterator cell, endc = mg_dof.end(level-1); + for (cell = mg_dof.begin(level-1); cell != endc; ++cell) + { + // need to look into a cell if it has children and it is locally + // owned + if (!cell->has_children()) + continue; + + bool consider_cell = false; + if (tria.locally_owned_subdomain()==numbers::invalid_subdomain_id + || cell->level_subdomain_id()==tria.locally_owned_subdomain() + ) + consider_cell = true; + + // due to the particular way we store DoF indices (via children), + // we also need to add the DoF indices for coarse cells where we + // own at least one child + bool cell_is_remote = !consider_cell; + for (unsigned int c=0; c::max_children_per_cell; ++c) + if (cell->child(c)->level_subdomain_id()==tria.locally_owned_subdomain()) + { + consider_cell = true; + break; + } + + if (!consider_cell) + continue; + + // step 2.2: loop through children and append the dof indices to + // the appropriate list. We need separate lists for the owned + // coarse cell case (which will be part of + // restriction/prolongation between level-1 and level) and the + // remote case (which needs to store DoF indices for the + // operations between level and level+1). + AssertDimension(cell->n_children(), + GeometryInfo::max_children_per_cell); + std::vector &next_indices = + cell_is_remote ? global_level_dof_indices_remote : global_level_dof_indices; + const std::size_t start_index = next_indices.size(); + next_indices.resize(start_index + elem_info.n_child_cell_dofs, + numbers::invalid_dof_index); + for (unsigned int c=0; c::max_children_per_cell; ++c) + { + if (cell_is_remote && cell->child(c)->level_subdomain_id() != + tria.locally_owned_subdomain()) + continue; + cell->child(c)->get_mg_dof_indices(local_dof_indices); + + const IndexSet &owned_level_dofs = mg_dof.locally_owned_mg_dofs(level); + for (unsigned int i=0; i(c, fe.dofs_per_cell - fe.dofs_per_vertex, + fe.degree, elem_info.shape_info.lexicographic_numbering, + local_dof_indices, + &next_indices[start_index]); + + // step 2.3 store the connectivity to the parent + if (cell->child(c)->has_children() && + (tria.locally_owned_subdomain()==numbers::invalid_subdomain_id + || cell->child(c)->level_subdomain_id()==tria.locally_owned_subdomain() + )) + { + const unsigned int child_index = coarse_level_indices[level][cell->child(c)->index()]; + AssertIndexRange(child_index, parent_child_connect[level].size()); + unsigned int parent_index = counter; + // remote cells, i.e., cells where we work on a further + // level but are not treated on the current level, need to + // be placed at the end of the list; however, we do not + // yet know the exact position in the array, so shift + // their parent index by the number of cells so we can set + // the correct number after the end of this loop + if (cell_is_remote) + parent_index = start_index/elem_info.n_child_cell_dofs + tria.n_cells(level); + parent_child_connect[level][child_index] = + std::make_pair(parent_index, c); + AssertIndexRange(mg_dof.get_fe().dofs_per_cell, + static_cast(-1)); + + // set Dirichlet boundary conditions (as a list of + // constrained DoFs) for the child + if (mg_constrained_dofs != 0) + for (unsigned int i=0; iis_boundary_index(level, + local_dof_indices[elem_info.shape_info.lexicographic_numbering[i]])) + dirichlet_indices[level][child_index].push_back(i); + } + } + if (!cell_is_remote) + { + AssertIndexRange(static_cast(cell->index()), + coarse_level_indices[level-1].size()); + coarse_level_indices[level-1][cell->index()] = counter++; + } + + // step 2.4: include indices for the coarsest cells. we still + // insert the indices as if they were from a child in order to use + // the same code (the coarsest level does not matter much in terms + // of memory, so we gain in code simplicity) + if (level == 1 && !cell_is_remote) + { + cell->get_mg_dof_indices(local_dof_indices); + + const IndexSet &owned_level_dofs_l0 = mg_dof.locally_owned_mg_dofs(0); + for (unsigned int i=0; i(0, fe.dofs_per_cell - fe.dofs_per_vertex, + fe.degree, elem_info.shape_info.lexicographic_numbering, + local_dof_indices, + &global_level_dof_indices_l0[start_index]); + + dirichlet_indices[0].push_back(std::vector()); + if (mg_constrained_dofs != 0) + for (unsigned int i=0; iis_boundary_index(0, local_dof_indices[elem_info.shape_info.lexicographic_numbering[i]])) + dirichlet_indices[0].back().push_back(i); + } + } + + // step 2.5: store information about the current level and prepare the + // Dirichlet indices and parent-child relationship for the next + // coarser level + AssertDimension(counter*elem_info.n_child_cell_dofs, global_level_dof_indices.size()); + n_owned_level_cells[level-1] = counter; + dirichlet_indices[level-1].resize(counter); + parent_child_connect[level-1]. + resize(counter, std::make_pair(numbers::invalid_unsigned_int, + numbers::invalid_unsigned_int)); + + // step 2.6: put the cells with remotely owned parent to the end of + // the list (these are needed for the transfer from level to level+1 + // but not for the transfer from level-1 to level). + if (level < n_levels-1) + for (std::vector >::iterator + i=parent_child_connect[level].begin(); i!=parent_child_connect[level].end(); ++i) + if (i->first >= tria.n_cells(level)) + { + i->first -= tria.n_cells(level); + i->first += counter; + } + + // step 2.7: Initialize the ghosted vector + const parallel::Triangulation *ptria = + (dynamic_cast*> (&tria)); + const MPI_Comm communicator = + ptria != 0 ? ptria->get_communicator() : MPI_COMM_SELF; + + reinit_ghosted_vector (mg_dof.locally_owned_mg_dofs(level), + ghosted_level_dofs, communicator, + ghosted_level_vector[level], + copy_indices_global_mine[level]); + + copy_indices_to_mpi_local_numbers(*ghosted_level_vector[level].get_partitioner(), + global_level_dof_indices, + global_level_dof_indices_remote, + level_dof_indices[level]); + // step 2.8: Initialize the ghosted vector for level 0 + if (level == 1) + { + for (unsigned int i = 0; i(), + level_dof_indices[0]); + } + } + + // ---------------------- 3. compute weights to make restriction additive + + // get the valence of the individual components and compute the weights as + // the inverse of the valence + weights_on_refined.resize(n_levels-1); + for (unsigned int level = 1; level::n_array_elements; + std::vector degree_to_3 (n_child_dofs_1d); + degree_to_3[0] = 0; + for (unsigned int i=1; i(3)); + for (unsigned int c=0; c2 ? n_child_dofs_1d : 1); ++k) + for (unsigned int j=0; j<(dim>1 ? n_child_dofs_1d : 1); ++j) + { + unsigned int shift = 9*degree_to_3[k] + 3*degree_to_3[j]; + for (unsigned int i=0; i(3)+shift+degree_to_3[i]][v] = Number(1.)/ + ghosted_level_vector[level].local_element(level_dof_indices[level][elem_info.n_child_cell_dofs*c+m]); + } + } + } + + } + } } diff --git a/source/multigrid/mg_transfer_internal.inst.in b/source/multigrid/mg_transfer_internal.inst.in index bab949b153..96c231acbd 100644 --- a/source/multigrid/mg_transfer_internal.inst.in +++ b/source/multigrid/mg_transfer_internal.inst.in @@ -35,3 +35,55 @@ for (deal_II_dimension : DIMENSIONS; \} \} } + + +for (deal_II_dimension : DIMENSIONS) +{ + namespace internal + \{ + namespace MGTransfer + \{ + + template + unsigned int + compute_shift_within_children(const unsigned int, + const unsigned int, + const unsigned int); + \} + \} +} + +for (S : REAL_SCALARS) +{ + namespace internal + \{ + namespace MGTransfer + \{ + + template + struct ElementInfo; + \} + \} +} + +for (deal_II_dimension : DIMENSIONS; S : REAL_SCALARS) +{ + namespace internal + \{ + namespace MGTransfer + \{ + + template + void setup_transfer(const dealii::DoFHandler&, + const MGConstrainedDoFs*, + ElementInfo&, + std::vector >&, + std::vector > >&, + std::vector&, + std::vector > >&, + std::vector > >&, + std::vector > >&, + MGLevelObject >&); + \} + \} +} diff --git a/source/multigrid/mg_transfer_matrix_free.cc b/source/multigrid/mg_transfer_matrix_free.cc index f8117424b3..e23881f1b3 100644 --- a/source/multigrid/mg_transfer_matrix_free.cc +++ b/source/multigrid/mg_transfer_matrix_free.cc @@ -22,10 +22,10 @@ #include #include #include -#include #include #include #include +#include #include #include @@ -91,453 +91,29 @@ void MGTransferMatrixFree::clear () } - -namespace -{ - // given the collection of child cells in lexicographic ordering as seen - // from the parent, compute the first index of the given child - template - unsigned int - compute_shift_within_children(const unsigned int child, - const unsigned int fe_shift_1d, - const unsigned int fe_degree) - { - // we put the degrees of freedom of all child cells in - // lexicographic ordering - unsigned int c_tensor_index[dim]; - unsigned int tmp = child; - for (unsigned int d=0; d - void add_child_indices(const unsigned int child, - const unsigned int fe_shift_1d, - const unsigned int fe_degree, - const std::vector &lexicographic_numbering, - const std::vector &local_dof_indices, - types::global_dof_index *target_indices) - { - const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d; - const unsigned int shift = - compute_shift_within_children(child, fe_shift_1d, fe_degree); - const unsigned int n_components = - local_dof_indices.size()/Utilities::fixed_power(fe_degree+1); - types::global_dof_index *indices = target_indices + shift; - const unsigned int n_scalar_cell_dofs = Utilities::fixed_power(n_child_dofs_1d); - for (unsigned int c=0, m=0; c2 ? (fe_degree+1) : 1); ++k) - for (unsigned int j=0; j<(dim>1 ? (fe_degree+1) : 1); ++j) - for (unsigned int i=0; i<(fe_degree+1); ++i, ++m) - { - const unsigned int index = c*n_scalar_cell_dofs+k*n_child_dofs_1d* - n_child_dofs_1d+j*n_child_dofs_1d+i; - Assert(indices[index] == numbers::invalid_dof_index || - indices[index] == local_dof_indices[lexicographic_numbering[m]], - ExcInternalError()); - indices[index] = local_dof_indices[lexicographic_numbering[m]]; - } - } - - - - // initialize the vectors needed for the transfer (and merge with the - // content in copy_indices_global_mine) - template - void - reinit_ghosted_vector(const IndexSet &locally_owned, - std::vector &ghosted_level_dofs, - const MPI_Comm &communicator, - LinearAlgebra::distributed::Vector &ghosted_level_vector, - std::vector > ©_indices_global_mine) - { - std::sort(ghosted_level_dofs.begin(), ghosted_level_dofs.end()); - IndexSet ghosted_dofs(locally_owned.size()); - ghosted_dofs.add_indices(ghosted_level_dofs.begin(), - std::unique(ghosted_level_dofs.begin(), - ghosted_level_dofs.end())); - ghosted_dofs.compress(); - - // Add possible ghosts from the previous content in the vector - if (ghosted_level_vector.size() == locally_owned.size()) - { - // shift the local number of the copy indices according to the new - // partitioner that we are going to use for the vector - const std_cxx11::shared_ptr part - = ghosted_level_vector.get_partitioner(); - ghosted_dofs.add_indices(part->ghost_indices()); - for (unsigned int i=0; ilocal_to_global(copy_indices_global_mine[i].second)); - } - ghosted_level_vector.reinit(locally_owned, ghosted_dofs, communicator); - } - - // Transform the ghost indices to local index space for the vector - void - copy_indices_to_mpi_local_numbers(const Utilities::MPI::Partitioner &part, - const std::vector &mine, - const std::vector &remote, - std::vector &localized_indices) - { - localized_indices.resize(mine.size()+remote.size(), - numbers::invalid_unsigned_int); - for (unsigned int i=0; i void MGTransferMatrixFree::build (const DoFHandler &mg_dof) { this->fill_and_communicate_copy_indices(mg_dof); - // we collect all child DoFs of a mother cell together. For faster - // tensorized operations, we align the degrees of freedom - // lexicographically. We distinguish FE_Q elements and FE_DGQ elements - - const Triangulation &tria = mg_dof.get_triangulation(); - - // ---------------------------- 1. Extract 1D info about the finite element - // step 1.1: create a 1D copy of the finite element from FETools where we - // substitute the template argument - AssertDimension(mg_dof.get_fe().n_base_elements(), 1); - std::string fe_name = mg_dof.get_fe().base_element(0).get_name(); - { - const std::size_t template_starts = fe_name.find_first_of('<'); - Assert (fe_name[template_starts+1] == (dim==1?'1':(dim==2?'2':'3')), - ExcInternalError()); - fe_name[template_starts+1] = '1'; - } - std_cxx11::shared_ptr > fe_1d - (FETools::get_fe_by_name<1,1>(fe_name)); - const FiniteElement<1> &fe = *fe_1d; - unsigned int n_child_dofs_1d = numbers::invalid_unsigned_int; - - { - // currently, we have only FE_Q and FE_DGQ type elements implemented - n_components = mg_dof.get_fe().element_multiplicity(0); - AssertDimension(Utilities::fixed_power(fe.dofs_per_cell)*n_components, - mg_dof.get_fe().dofs_per_cell); - AssertDimension(fe.degree, mg_dof.get_fe().degree); - fe_degree = fe.degree; - element_is_continuous = fe.dofs_per_vertex > 0; - Assert(fe.dofs_per_vertex < 2, ExcNotImplemented()); - - // step 1.2: get renumbering of 1D basis functions to lexicographic - // numbers. The distinction according to fe.dofs_per_vertex is to support - // both continuous and discontinuous bases. - std::vector renumbering(fe.dofs_per_cell); - { - AssertIndexRange(fe.dofs_per_vertex, 2); - renumbering[0] = 0; - for (unsigned int i=0; i::vertices_per_cell*fe.dofs_per_vertex + i; - if (fe.dofs_per_vertex > 0) - renumbering[fe.dofs_per_cell-fe.dofs_per_vertex] = fe.dofs_per_vertex; - } - - // step 1.3: create a 1D quadrature formula from the finite element that - // collects the support points of the basis functions on the two children. - std::vector > basic_support_points = fe.get_unit_support_points(); - Assert(fe.dofs_per_vertex == 0 || fe.dofs_per_vertex == 1, - ExcNotImplemented()); - std::vector > points_refined(fe.dofs_per_vertex > 0 ? - (2 * fe.dofs_per_cell - 1) : - (2 * fe.dofs_per_cell)); - const unsigned int shift = fe.dofs_per_cell - fe.dofs_per_vertex; - for (unsigned int c=0; c::max_children_per_cell; ++c) - for (unsigned int j=0; j(n_child_dofs_1d); - - // step 1.4: evaluate the polynomials and store the data in ShapeInfo - const Quadrature<1> quadrature(points_refined); - shape_info.reinit(quadrature, mg_dof.get_fe(), 0); - - for (unsigned int c=0; c::max_children_per_cell; ++c) - for (unsigned int i=0; i::epsilon(),1e-12), - ExcInternalError()); - } - - // -------------- 2. Extract and match dof indices between child and parent - const unsigned int n_levels = tria.n_global_levels(); - level_dof_indices.resize(n_levels); - parent_child_connect.resize(n_levels-1); - n_owned_level_cells.resize(n_levels-1); - std::vector > coarse_level_indices(n_levels-1); - for (unsigned int level=0; level local_dof_indices(mg_dof.get_fe().dofs_per_cell); - dirichlet_indices.resize(n_levels-1); - - // We use the vectors stored ghosted_level_vector in the base class for - // keeping ghosted transfer indices. To avoid keeping two very similar - // vectors, we merge them here. - if (this->ghosted_level_vector.max_level() != n_levels-1) - this->ghosted_level_vector.resize(0, n_levels-1); - - for (unsigned int level=n_levels-1; level > 0; --level) - { - unsigned int counter = 0; - std::vector global_level_dof_indices; - std::vector global_level_dof_indices_remote; - std::vector ghosted_level_dofs; - std::vector global_level_dof_indices_l0; - std::vector ghosted_level_dofs_l0; - - // step 2.1: loop over the cells on the coarse side - typename DoFHandler::cell_iterator cell, endc = mg_dof.end(level-1); - for (cell = mg_dof.begin(level-1); cell != endc; ++cell) - { - // need to look into a cell if it has children and it is locally owned - if (!cell->has_children()) - continue; - - bool consider_cell = false; - if (tria.locally_owned_subdomain()==numbers::invalid_subdomain_id - || cell->level_subdomain_id()==tria.locally_owned_subdomain() - ) - consider_cell = true; - - // due to the particular way we store DoF indices (via children), we - // also need to add the DoF indices for coarse cells where we own at - // least one child - bool cell_is_remote = !consider_cell; - for (unsigned int c=0; c::max_children_per_cell; ++c) - if (cell->child(c)->level_subdomain_id()==tria.locally_owned_subdomain()) - { - consider_cell = true; - break; - } - - if (!consider_cell) - continue; - - // step 2.2: loop through children and append the dof indices to the - // appropriate list. We need separate lists for the owned coarse - // cell case (which will be part of restriction/prolongation between - // level-1 and level) and the remote case (which needs to store DoF - // indices for the operations between level and level+1). - AssertDimension(cell->n_children(), - GeometryInfo::max_children_per_cell); - std::vector &next_indices = - cell_is_remote ? global_level_dof_indices_remote : global_level_dof_indices; - const std::size_t start_index = next_indices.size(); - next_indices.resize(start_index + n_child_cell_dofs, - numbers::invalid_dof_index); - for (unsigned int c=0; c::max_children_per_cell; ++c) - { - if (cell_is_remote && cell->child(c)->level_subdomain_id() != - tria.locally_owned_subdomain()) - continue; - cell->child(c)->get_mg_dof_indices(local_dof_indices); - - const IndexSet &owned_level_dofs = mg_dof.locally_owned_mg_dofs(level); - for (unsigned int i=0; i(c, fe.dofs_per_cell - fe.dofs_per_vertex, - fe.degree, shape_info.lexicographic_numbering, - local_dof_indices, - &next_indices[start_index]); - - // step 2.3 store the connectivity to the parent - if (cell->child(c)->has_children() && - (tria.locally_owned_subdomain()==numbers::invalid_subdomain_id - || cell->child(c)->level_subdomain_id()==tria.locally_owned_subdomain() - )) - { - const unsigned int child_index = coarse_level_indices[level][cell->child(c)->index()]; - AssertIndexRange(child_index, parent_child_connect[level].size()); - unsigned int parent_index = counter; - // remote cells, i.e., cells where we work on a further - // level but are not treated on the current level, need to - // be placed at the end of the list; however, we do not yet - // know the exact position in the array, so shift their - // parent index by the number of cells so we can set the - // correct number after the end of this loop - if (cell_is_remote) - parent_index = start_index/n_child_cell_dofs + tria.n_cells(level); - parent_child_connect[level][child_index] = - std::make_pair(parent_index, c); - AssertIndexRange(mg_dof.get_fe().dofs_per_cell, - static_cast(-1)); - - // set Dirichlet boundary conditions (as a list of - // constrained DoFs) for the child - if (this->mg_constrained_dofs != 0) - for (unsigned int i=0; img_constrained_dofs->is_boundary_index(level, local_dof_indices[shape_info.lexicographic_numbering[i]])) - dirichlet_indices[level][child_index].push_back(i); - } - } - if (!cell_is_remote) - { - AssertIndexRange(static_cast(cell->index()), - coarse_level_indices[level-1].size()); - coarse_level_indices[level-1][cell->index()] = counter++; - } - - // step 2.4: include indices for the coarsest cells. we still insert - // the indices as if they were from a child in order to use the same - // code (the coarsest level does not matter much in terms of memory, - // so we gain in code simplicity) - if (level == 1 && !cell_is_remote) - { - cell->get_mg_dof_indices(local_dof_indices); - - const IndexSet &owned_level_dofs_l0 = mg_dof.locally_owned_mg_dofs(0); - for (unsigned int i=0; i(0, fe.dofs_per_cell - fe.dofs_per_vertex, - fe.degree, shape_info.lexicographic_numbering, - local_dof_indices, - &global_level_dof_indices_l0[start_index]); - - dirichlet_indices[0].push_back(std::vector()); - if (this->mg_constrained_dofs != 0) - for (unsigned int i=0; img_constrained_dofs->is_boundary_index(0, local_dof_indices[shape_info.lexicographic_numbering[i]])) - dirichlet_indices[0].back().push_back(i); - } - } - - // step 2.5: store information about the current level and prepare the - // Dirichlet indices and parent-child relationship for the next coarser - // level - AssertDimension(counter*n_child_cell_dofs, global_level_dof_indices.size()); - n_owned_level_cells[level-1] = counter; - dirichlet_indices[level-1].resize(counter); - parent_child_connect[level-1]. - resize(counter, std::make_pair(numbers::invalid_unsigned_int, - numbers::invalid_unsigned_int)); - - // step 2.6: put the cells with remotely owned parent to the end of the - // list (these are needed for the transfer from level to level+1 but not - // for the transfer from level-1 to level). - if (level < n_levels-1) - for (std::vector >::iterator - i=parent_child_connect[level].begin(); i!=parent_child_connect[level].end(); ++i) - if (i->first >= tria.n_cells(level)) - { - i->first -= tria.n_cells(level); - i->first += counter; - } - - // step 2.7: Initialize the ghosted vector - const parallel::Triangulation *ptria = - (dynamic_cast*> (&tria)); - const MPI_Comm communicator = - ptria != 0 ? ptria->get_communicator() : MPI_COMM_SELF; - - reinit_ghosted_vector(mg_dof.locally_owned_mg_dofs(level), - ghosted_level_dofs, communicator, - this->ghosted_level_vector[level], - this->copy_indices_global_mine[level]); - - copy_indices_to_mpi_local_numbers(*this->ghosted_level_vector[level].get_partitioner(), - global_level_dof_indices, - global_level_dof_indices_remote, - level_dof_indices[level]); - - // step 2.8: Initialize the ghosted vector for level 0 - if (level == 1) - { - for (unsigned int i = 0; ighosted_level_vector[0], - this->copy_indices_global_mine[0]); - - copy_indices_to_mpi_local_numbers(*this->ghosted_level_vector[0].get_partitioner(), - global_level_dof_indices_l0, - std::vector(), - level_dof_indices[0]); - } - } - - // ------------------------ 3. compute weights to make restriction additive - // - // get the valence of the individual components and compute the weights as - // the inverse of the valence - weights_on_refined.resize(n_levels-1); - for (unsigned int level = 1; levelghosted_level_vector[level] = 0; - for (unsigned int c=0; cghosted_level_vector[level].local_element(level_dof_indices[level][n_child_cell_dofs*c+j]) += Number(1.); - this->ghosted_level_vector[level].compress(VectorOperation::add); - this->ghosted_level_vector[level].update_ghost_values(); - - const unsigned int vec_size = VectorizedArray::n_array_elements; - std::vector degree_to_3 (n_child_dofs_1d); - degree_to_3[0] = 0; - for (unsigned int i=1; i(3)); - for (unsigned int c=0; c2 ? n_child_dofs_1d : 1); ++k) - for (unsigned int j=0; j<(dim>1 ? n_child_dofs_1d : 1); ++j) - { - unsigned int shift = 9*degree_to_3[k] + 3*degree_to_3[j]; - for (unsigned int i=0; i(3)+shift+degree_to_3[i]][v] = Number(1.)/ - this->ghosted_level_vector[level].local_element(level_dof_indices[level][n_child_cell_dofs*c+m]); - } - } - } + internal::MGTransfer::ElementInfo elem_info; + internal::MGTransfer::setup_transfer(mg_dof, + (const MGConstrainedDoFs *)this->mg_constrained_dofs, + elem_info, + level_dof_indices, + parent_child_connect, + n_owned_level_cells, + dirichlet_indices, + weights_on_refined, + this->copy_indices_global_mine, + this->ghosted_level_vector); + // unpack element info data + fe_degree = elem_info.fe_degree; + element_is_continuous = elem_info.element_is_continuous; + n_components = elem_info.n_components; + n_child_cell_dofs = elem_info.n_child_cell_dofs; + shape_info = elem_info.shape_info; evaluation_data.resize(3*n_child_cell_dofs); } @@ -764,7 +340,7 @@ void MGTransferMatrixFree // read from source vector for (unsigned int v=0; v + const unsigned int shift = internal::MGTransfer::compute_shift_within_children (parent_child_connect[to_level-1][cell+v].second, degree+1-element_is_continuous, degree); const unsigned int *indices = &level_dof_indices[to_level-1][parent_child_connect[to_level-1][cell+v].first*n_child_cell_dofs+shift]; @@ -897,7 +473,7 @@ void MGTransferMatrixFree // write into dst vector for (unsigned int v=0; v + const unsigned int shift = internal::MGTransfer::compute_shift_within_children (parent_child_connect[from_level-1][cell+v].second, degree+1-element_is_continuous, degree); AssertIndexRange(parent_child_connect[from_level-1][cell+v].first* -- 2.39.5