#ifndef dealii__mg_transfer_internal_h
#define dealii__mg_transfer_internal_h
+#include <deal.II/base/mg_level_object.h>
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>
-
DEAL_II_NAMESPACE_OPEN
namespace internal
std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > > ©_indices,
std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > > ©_indices_global_mine,
std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > > ©_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 <int dim>
+ 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
+ * <tt>setup_transfer</tt>.
+ */
+ template <typename Number>
+ 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
+ * <tt>2<sup>dim</sup>*fe.dofs_per_cell</tt> 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<Number> shape_info;
+
+ };
+
+ /**
+ * Sets up most of the internal data structures of MGTransferMatrixFree
+ */
+ template <int dim, typename Number>
+ void setup_transfer(const dealii::DoFHandler<dim> &mg_dof,
+ const MGConstrainedDoFs *mg_constrained_dofs,
+ ElementInfo<Number> &elem_info,
+ std::vector<std::vector<unsigned int> > &level_dof_indices,
+ std::vector<std::vector<std::pair<unsigned int,unsigned int> > > &parent_child_connect,
+ std::vector<unsigned int> &n_owned_level_cells,
+ std::vector<std::vector<std::vector<unsigned short> > > &dirichlet_indices,
+ std::vector<AlignedVector<VectorizedArray<Number> > > &weights_on_refined,
+ std::vector<std::vector<std::pair<unsigned int, unsigned int> > > ©_indices_global_mine,
+ MGLevelObject<LinearAlgebra::distributed::Vector<Number> > &ghosted_level_vector);
+
}
}
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_tools.h>
#include <deal.II/multigrid/mg_transfer_internal.h>
DEAL_II_NAMESPACE_OPEN
{
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;
{}
};
+ // Internal function for filling the copy indices from global to level
+ // indices
template <int dim, int spacedim>
void fill_copy_indices(const dealii::DoFHandler<dim,spacedim> &mg_dof,
const MGConstrainedDoFs *mg_constrained_dofs,
for (unsigned int level=0; level<copy_indices_global_mine.size(); ++level)
std::sort(copy_indices_global_mine[level].begin(), copy_indices_global_mine[level].end(), compare);
}
+
+
+
+ // initialize the vectors needed for the transfer (and merge with the
+ // content in copy_indices_global_mine)
+ template <typename Number>
+ void
+ reinit_ghosted_vector(const IndexSet &locally_owned,
+ std::vector<types::global_dof_index> &ghosted_level_dofs,
+ const MPI_Comm &communicator,
+ LinearAlgebra::distributed::Vector<Number> &ghosted_level_vector,
+ std::vector<std::pair<unsigned int,unsigned int> > ©_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<const Utilities::MPI::Partitioner> part
+ = ghosted_level_vector.get_partitioner();
+ ghosted_dofs.add_indices(part->ghost_indices());
+ for (unsigned int i=0; i<copy_indices_global_mine.size(); ++i)
+ copy_indices_global_mine[i].second =
+ locally_owned.n_elements() +
+ ghosted_dofs.index_within_set(part->local_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<types::global_dof_index> &mine,
+ const std::vector<types::global_dof_index> &remote,
+ std::vector<unsigned int> &localized_indices)
+ {
+ localized_indices.resize(mine.size()+remote.size(),
+ numbers::invalid_unsigned_int);
+ for (unsigned int i=0; i<mine.size(); ++i)
+ if (mine[i] != numbers::invalid_dof_index)
+ localized_indices[i] = part.global_to_local(mine[i]);
+
+ for (unsigned int i=0; i<remote.size(); ++i)
+ if (remote[i] != numbers::invalid_dof_index)
+ localized_indices[i+mine.size()] = part.global_to_local(remote[i]);
+ }
+
+ // given the collection of child cells in lexicographic ordering as seen
+ // from the parent, compute the first index of the given child
+ template <int dim>
+ 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<dim; ++d)
+ {
+ c_tensor_index[d] = tmp % 2;
+ tmp /= 2;
+ }
+ const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
+ unsigned int factor = 1;
+ unsigned int shift = fe_shift_1d * c_tensor_index[0];
+ for (unsigned int d=1; d<dim; ++d)
+ {
+ factor *= n_child_dofs_1d;
+ shift = shift + factor * fe_shift_1d * c_tensor_index[d];
+ }
+ return shift;
+ }
+
+ // puts the indices on the given child cell in lexicographic ordering with
+ // respect to the collection of all child cells as seen from the parent
+ template <int dim>
+ void add_child_indices(const unsigned int child,
+ const unsigned int fe_shift_1d,
+ const unsigned int fe_degree,
+ const std::vector<unsigned int> &lexicographic_numbering,
+ const std::vector<types::global_dof_index> &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<dim>(child, fe_shift_1d, fe_degree);
+ const unsigned int n_components =
+ local_dof_indices.size()/Utilities::fixed_power<dim>(fe_degree+1);
+ types::global_dof_index *indices = target_indices + shift;
+ const unsigned int n_scalar_cell_dofs = Utilities::fixed_power<dim>(n_child_dofs_1d);
+ for (unsigned int c=0, m=0; c<n_components; ++c)
+ for (unsigned int k=0; k<(dim>2 ? (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 <int dim, typename Number>
+ void setup_element_info(ElementInfo<Number> &elem_info,const FiniteElement<1> &fe,
+ const dealii::DoFHandler<dim> &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<dim>(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<unsigned int> renumbering(fe.dofs_per_cell);
+ {
+ AssertIndexRange(fe.dofs_per_vertex, 2);
+ renumbering[0] = 0;
+ for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+ renumbering[i+fe.dofs_per_vertex] =
+ GeometryInfo<1>::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<Point<1> > basic_support_points = fe.get_unit_support_points();
+ Assert(fe.dofs_per_vertex == 0 || fe.dofs_per_vertex == 1,
+ ExcNotImplemented());
+ std::vector<Point<1> > 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<GeometryInfo<1>::max_children_per_cell; ++c)
+ for (unsigned int j=0; j<basic_support_points.size(); ++j)
+ points_refined[shift*c+j][0] =
+ c*0.5 + 0.5 * basic_support_points[renumbering[j]][0];
+
+ unsigned int n_child_dofs_1d = points_refined.size();
+ elem_info.n_child_cell_dofs = elem_info.n_components*Utilities::fixed_power<dim>(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<GeometryInfo<1>::max_children_per_cell; ++c)
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ Assert(std::abs(elem_info.shape_info.shape_values[i*n_child_dofs_1d+j+c*shift][0] -
+ fe.get_prolongation_matrix(c)(renumbering[j],renumbering[i]))
+ < std::max(2.*(double)std::numeric_limits<Number>::epsilon(),1e-12),
+ ExcInternalError());
+ }
+
+
+ // Sets up most of the internal data structures of the MGTransferMatrixFree
+ // class
+ template <int dim, typename Number>
+ void setup_transfer(const dealii::DoFHandler<dim> &mg_dof,
+ const MGConstrainedDoFs *mg_constrained_dofs,
+ ElementInfo<Number> &elem_info,
+ std::vector<std::vector<unsigned int> > &level_dof_indices,
+ std::vector<std::vector<std::pair<unsigned int,unsigned int> > > &parent_child_connect,
+ std::vector<unsigned int> &n_owned_level_cells,
+ std::vector<std::vector<std::vector<unsigned short> > > &dirichlet_indices,
+ std::vector<AlignedVector<VectorizedArray<Number> > > &weights_on_refined,
+ std::vector<std::vector<std::pair<unsigned int, unsigned int> > > ©_indices_global_mine,
+ MGLevelObject<LinearAlgebra::distributed::Vector<Number> > &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<dim> &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<FiniteElement<1> > 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<std::vector<unsigned int> > coarse_level_indices(n_levels-1);
+ for (unsigned int level=0; level<std::min(tria.n_levels(),n_levels-1); ++level)
+ coarse_level_indices[level].resize(tria.n_raw_cells(level),
+ numbers::invalid_unsigned_int);
+ std::vector<types::global_dof_index> 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<types::global_dof_index> global_level_dof_indices;
+ std::vector<types::global_dof_index> global_level_dof_indices_remote;
+ std::vector<types::global_dof_index> ghosted_level_dofs;
+ std::vector<types::global_dof_index> global_level_dof_indices_l0;
+ std::vector<types::global_dof_index> ghosted_level_dofs_l0;
+
+ // step 2.1: loop over the cells on the coarse side
+ typename dealii::DoFHandler<dim>::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<GeometryInfo<dim>::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<dim>::max_children_per_cell);
+ std::vector<types::global_dof_index> &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<GeometryInfo<dim>::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<local_dof_indices.size(); ++i)
+ if (!owned_level_dofs.is_element(local_dof_indices[i]))
+ ghosted_level_dofs.push_back(local_dof_indices[i]);
+
+ add_child_indices<dim>(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<unsigned short>(-1));
+
+ // set Dirichlet boundary conditions (as a list of
+ // constrained DoFs) for the child
+ if (mg_constrained_dofs != 0)
+ for (unsigned int i=0; i<mg_dof.get_fe().dofs_per_cell; ++i)
+ if (mg_constrained_dofs->is_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<unsigned int>(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<local_dof_indices.size(); ++i)
+ if (!owned_level_dofs_l0.is_element(local_dof_indices[i]))
+ ghosted_level_dofs_l0.push_back(local_dof_indices[i]);
+
+ const std::size_t start_index = global_level_dof_indices_l0.size();
+ global_level_dof_indices_l0.resize(start_index+elem_info.n_child_cell_dofs,
+ numbers::invalid_dof_index);
+ add_child_indices<dim>(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<unsigned short>());
+ if (mg_constrained_dofs != 0)
+ for (unsigned int i=0; i<mg_dof.get_fe().dofs_per_cell; ++i)
+ if (mg_constrained_dofs->is_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<std::pair<unsigned int,unsigned int> >::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<dim,dim> *ptria =
+ (dynamic_cast<const parallel::Triangulation<dim,dim>*> (&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<parent_child_connect[0].size(); ++i)
+ parent_child_connect[0][i] = std::make_pair(i, 0U);
+
+ reinit_ghosted_vector (mg_dof.locally_owned_mg_dofs(0),
+ ghosted_level_dofs_l0, communicator,
+ ghosted_level_vector[0],
+ copy_indices_global_mine[0]);
+
+ copy_indices_to_mpi_local_numbers(*ghosted_level_vector[0].get_partitioner(),
+ global_level_dof_indices_l0,
+ std::vector<types::global_dof_index>(),
+ 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_levels; ++level)
+ {
+ ghosted_level_vector[level] = 0;
+ for (unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
+ for (unsigned int j=0; j<elem_info.n_child_cell_dofs; ++j)
+ ghosted_level_vector[level].local_element(level_dof_indices[level][elem_info.n_child_cell_dofs*c+j]) += Number(1.);
+ ghosted_level_vector[level].compress(VectorOperation::add);
+ ghosted_level_vector[level].update_ghost_values();
+
+ const unsigned int vec_size = VectorizedArray<Number>::n_array_elements;
+ std::vector<unsigned int> degree_to_3 (n_child_dofs_1d);
+ degree_to_3[0] = 0;
+ for (unsigned int i=1; i<n_child_dofs_1d-1; ++i)
+ degree_to_3[i] = 1;
+ degree_to_3.back() = 2;
+
+ // we only store 3^dim weights because all dofs on a line have the
+ // same valence, and all dofs on a quad have the same valence.
+ weights_on_refined[level-1].resize(((n_owned_level_cells[level-1]+vec_size-1)/vec_size)*Utilities::fixed_power<dim>(3));
+ for (unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
+ {
+ const unsigned int comp = c/vec_size;
+ const unsigned int v = c%vec_size;
+
+ for (unsigned int k=0, m=0; k<(dim>2 ? 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<n_child_dofs_1d; ++i, ++m)
+ weights_on_refined[level-1][comp*Utilities::fixed_power<dim>(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]);
+ }
+ }
+ }
+
+ }
+
}
}
\}
\}
}
+
+
+for (deal_II_dimension : DIMENSIONS)
+{
+ namespace internal
+ \{
+ namespace MGTransfer
+ \{
+
+ template
+ unsigned int
+ compute_shift_within_children<deal_II_dimension>(const unsigned int,
+ const unsigned int,
+ const unsigned int);
+ \}
+ \}
+}
+
+for (S : REAL_SCALARS)
+{
+ namespace internal
+ \{
+ namespace MGTransfer
+ \{
+
+ template
+ struct ElementInfo<S>;
+ \}
+ \}
+}
+
+for (deal_II_dimension : DIMENSIONS; S : REAL_SCALARS)
+{
+ namespace internal
+ \{
+ namespace MGTransfer
+ \{
+
+ template
+ void setup_transfer<deal_II_dimension>(const dealii::DoFHandler<deal_II_dimension>&,
+ const MGConstrainedDoFs*,
+ ElementInfo<S>&,
+ std::vector<std::vector<unsigned int> >&,
+ std::vector<std::vector<std::pair<unsigned int,unsigned int> > >&,
+ std::vector<unsigned int>&,
+ std::vector<std::vector<std::vector<unsigned short> > >&,
+ std::vector<AlignedVector<VectorizedArray<S> > >&,
+ std::vector<std::vector<std::pair<unsigned int, unsigned int> > >&,
+ MGLevelObject<LinearAlgebra::distributed::Vector<S> >&);
+ \}
+ \}
+}
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe.h>
-#include <deal.II/fe/fe_tools.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/multigrid/mg_tools.h>
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+#include <deal.II/multigrid/mg_transfer_internal.h>
#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/matrix_free/fe_evaluation.h>
}
-
-namespace
-{
- // given the collection of child cells in lexicographic ordering as seen
- // from the parent, compute the first index of the given child
- template <int dim>
- 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<dim; ++d)
- {
- c_tensor_index[d] = tmp % 2;
- tmp /= 2;
- }
- const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
- unsigned int factor = 1;
- unsigned int shift = fe_shift_1d * c_tensor_index[0];
- for (unsigned int d=1; d<dim; ++d)
- {
- factor *= n_child_dofs_1d;
- shift = shift + factor * fe_shift_1d * c_tensor_index[d];
- }
- return shift;
- }
-
-
-
- // puts the indices on the given child cell in lexicographic ordering with
- // respect to the collection of all child cells as seen from the parent
- template <int dim>
- void add_child_indices(const unsigned int child,
- const unsigned int fe_shift_1d,
- const unsigned int fe_degree,
- const std::vector<unsigned int> &lexicographic_numbering,
- const std::vector<types::global_dof_index> &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<dim>(child, fe_shift_1d, fe_degree);
- const unsigned int n_components =
- local_dof_indices.size()/Utilities::fixed_power<dim>(fe_degree+1);
- types::global_dof_index *indices = target_indices + shift;
- const unsigned int n_scalar_cell_dofs = Utilities::fixed_power<dim>(n_child_dofs_1d);
- for (unsigned int c=0, m=0; c<n_components; ++c)
- for (unsigned int k=0; k<(dim>2 ? (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 <typename Number>
- void
- reinit_ghosted_vector(const IndexSet &locally_owned,
- std::vector<types::global_dof_index> &ghosted_level_dofs,
- const MPI_Comm &communicator,
- LinearAlgebra::distributed::Vector<Number> &ghosted_level_vector,
- std::vector<std::pair<unsigned int,unsigned int> > ©_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<const Utilities::MPI::Partitioner> part
- = ghosted_level_vector.get_partitioner();
- ghosted_dofs.add_indices(part->ghost_indices());
- for (unsigned int i=0; i<copy_indices_global_mine.size(); ++i)
- copy_indices_global_mine[i].second =
- locally_owned.n_elements() +
- ghosted_dofs.index_within_set(part->local_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<types::global_dof_index> &mine,
- const std::vector<types::global_dof_index> &remote,
- std::vector<unsigned int> &localized_indices)
- {
- localized_indices.resize(mine.size()+remote.size(),
- numbers::invalid_unsigned_int);
- for (unsigned int i=0; i<mine.size(); ++i)
- if (mine[i] != numbers::invalid_dof_index)
- localized_indices[i] = part.global_to_local(mine[i]);
-
- for (unsigned int i=0; i<remote.size(); ++i)
- if (remote[i] != numbers::invalid_dof_index)
- localized_indices[i+mine.size()] = part.global_to_local(remote[i]);
- }
-}
-
-
-
template <int dim, typename Number>
void MGTransferMatrixFree<dim,Number>::build
(const DoFHandler<dim,dim> &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<dim> &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<FiniteElement<1> > 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<dim>(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<unsigned int> renumbering(fe.dofs_per_cell);
- {
- AssertIndexRange(fe.dofs_per_vertex, 2);
- renumbering[0] = 0;
- for (unsigned int i=0; i<fe.dofs_per_line; ++i)
- renumbering[i+fe.dofs_per_vertex] =
- GeometryInfo<1>::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<Point<1> > basic_support_points = fe.get_unit_support_points();
- Assert(fe.dofs_per_vertex == 0 || fe.dofs_per_vertex == 1,
- ExcNotImplemented());
- std::vector<Point<1> > 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<GeometryInfo<1>::max_children_per_cell; ++c)
- for (unsigned int j=0; j<basic_support_points.size(); ++j)
- points_refined[shift*c+j][0] =
- c*0.5 + 0.5 * basic_support_points[renumbering[j]][0];
-
- n_child_dofs_1d = points_refined.size();
- n_child_cell_dofs = n_components*Utilities::fixed_power<dim>(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<GeometryInfo<1>::max_children_per_cell; ++c)
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
- Assert(std::abs(shape_info.shape_values[i*n_child_dofs_1d+j+c*shift][0] -
- fe.get_prolongation_matrix(c)(renumbering[j],renumbering[i]))
- < std::max(2.*(double)std::numeric_limits<Number>::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<std::vector<unsigned int> > coarse_level_indices(n_levels-1);
- for (unsigned int level=0; level<std::min(tria.n_levels(),n_levels-1); ++level)
- coarse_level_indices[level].resize(tria.n_raw_cells(level),
- numbers::invalid_unsigned_int);
- std::vector<types::global_dof_index> 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<types::global_dof_index> global_level_dof_indices;
- std::vector<types::global_dof_index> global_level_dof_indices_remote;
- std::vector<types::global_dof_index> ghosted_level_dofs;
- std::vector<types::global_dof_index> global_level_dof_indices_l0;
- std::vector<types::global_dof_index> ghosted_level_dofs_l0;
-
- // step 2.1: loop over the cells on the coarse side
- typename DoFHandler<dim>::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<GeometryInfo<dim>::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<dim>::max_children_per_cell);
- std::vector<types::global_dof_index> &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<GeometryInfo<dim>::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<local_dof_indices.size(); ++i)
- if (!owned_level_dofs.is_element(local_dof_indices[i]))
- ghosted_level_dofs.push_back(local_dof_indices[i]);
-
- add_child_indices<dim>(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<unsigned short>(-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; i<mg_dof.get_fe().dofs_per_cell; ++i)
- if (this->mg_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<unsigned int>(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<local_dof_indices.size(); ++i)
- if (!owned_level_dofs_l0.is_element(local_dof_indices[i]))
- ghosted_level_dofs_l0.push_back(local_dof_indices[i]);
-
- const std::size_t start_index = global_level_dof_indices_l0.size();
- global_level_dof_indices_l0.resize(start_index+n_child_cell_dofs,
- numbers::invalid_dof_index);
- add_child_indices<dim>(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<unsigned short>());
- if (this->mg_constrained_dofs != 0)
- for (unsigned int i=0; i<mg_dof.get_fe().dofs_per_cell; ++i)
- if (this->mg_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<std::pair<unsigned int,unsigned int> >::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<dim,dim> *ptria =
- (dynamic_cast<const parallel::Triangulation<dim,dim>*> (&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; i<parent_child_connect[0].size(); ++i)
- parent_child_connect[0][i] = std::make_pair(i, 0U);
-
- reinit_ghosted_vector(mg_dof.locally_owned_mg_dofs(0),
- ghosted_level_dofs_l0, communicator,
- this->ghosted_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<types::global_dof_index>(),
- 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_levels; ++level)
- {
- this->ghosted_level_vector[level] = 0;
- for (unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
- for (unsigned int j=0; j<n_child_cell_dofs; ++j)
- this->ghosted_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<Number>::n_array_elements;
- std::vector<unsigned int> degree_to_3 (n_child_dofs_1d);
- degree_to_3[0] = 0;
- for (unsigned int i=1; i<n_child_dofs_1d-1; ++i)
- degree_to_3[i] = 1;
- degree_to_3.back() = 2;
-
- // we only store 3^dim weights because all dofs on a line have the same
- // valence, and all dofs on a quad have the same valence.
- weights_on_refined[level-1].resize(((n_owned_level_cells[level-1]+vec_size-1)/vec_size)*Utilities::fixed_power<dim>(3));
- for (unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
- {
- const unsigned int comp = c/vec_size;
- const unsigned int v = c%vec_size;
-
- for (unsigned int k=0, m=0; k<(dim>2 ? 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<n_child_dofs_1d; ++i, ++m)
- weights_on_refined[level-1][comp*Utilities::fixed_power<dim>(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<Number> elem_info;
+ internal::MGTransfer::setup_transfer<dim,Number>(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);
}
// read from source vector
for (unsigned int v=0; v<n_chunks; ++v)
{
- const unsigned int shift = compute_shift_within_children<dim>
+ const unsigned int shift = internal::MGTransfer::compute_shift_within_children<dim>
(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];
// write into dst vector
for (unsigned int v=0; v<n_chunks; ++v)
{
- const unsigned int shift = compute_shift_within_children<dim>
+ const unsigned int shift = internal::MGTransfer::compute_shift_within_children<dim>
(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*
projects / dealii.git / commitdiff