const AffineConstraints<double> &
get_user_constraint_matrix(const unsigned int level) const;
+ /**
+ * Merge selected constraints of a specifiedlevel into a given single
+ * AffineConstraints object.
+ *
+ * @param constraints AffineConstraints object to be filled.
+ * @param level Refinement to be considered.
+ * @param add_boundary_indices Add boundary indices.
+ * @param add_refinement_edge_indices Add refinement-edge indices.
+ * @param add_level_constraints Add level constraints including the one passed
+ * during initialize() and periodicy constraints.
+ * @param add_user_constraints Add user constraints.
+ */
+ template <typename Number>
+ void
+ merge_constraints(AffineConstraints<Number> &constraints,
+ const unsigned int level,
+ const bool add_boundary_indices,
+ const bool add_refinement_edge_indices,
+ const bool add_level_constraints,
+ const bool add_user_constraints) const;
+
private:
/**
* The indices of boundary dofs for each level.
}
+template <typename Number>
+inline void
+MGConstrainedDoFs::merge_constraints(AffineConstraints<Number> &constraints,
+ const unsigned int level,
+ const bool add_boundary_indices,
+ const bool add_refinement_edge_indices,
+ const bool add_level_constraints,
+ const bool add_user_constraints) const
+{
+ constraints.clear();
+
+ // determine local lines
+ IndexSet index_set(this->get_refinement_edge_indices(level).size());
+
+ if (add_boundary_indices && this->have_boundary_indices())
+ index_set.add_indices(this->get_boundary_indices(level));
+
+ if (add_refinement_edge_indices)
+ index_set.add_indices(this->get_refinement_edge_indices(level));
+
+ if (add_level_constraints)
+ index_set.add_indices(this->get_level_constraints(level).get_local_lines());
+
+ if (add_user_constraints)
+ index_set.add_indices(
+ this->get_user_constraint_matrix(level).get_local_lines());
+
+ constraints.reinit(index_set);
+
+ // merge constraints
+ if (add_boundary_indices && this->have_boundary_indices())
+ constraints.add_lines(this->get_boundary_indices(level));
+
+ if (add_refinement_edge_indices)
+ constraints.add_lines(this->get_refinement_edge_indices(level));
+
+ if (add_level_constraints)
+ constraints.merge(this->get_level_constraints(level),
+ AffineConstraints<Number>::left_object_wins,
+ true);
+
+ if (add_user_constraints)
+ constraints.merge(this->get_user_constraint_matrix(level),
+ AffineConstraints<Number>::left_object_wins,
+ true);
+
+ // finalize setup
+ constraints.close();
+}
+
+
DEAL_II_NAMESPACE_CLOSE
/**
* Implementation of the MGTransferBase. In contrast to
* other multigrid transfer operators, the user can provide separate
- * transfer operators of type MGTwoLevelTransfer between each level.
+ * transfer operators of type MGTwoLevelTransfer between each level. The
+ * sequence of functions calls for setup is:
+ * @code
+ * MGTransferGlobalCoarsening mg_transfer;
+ * mg_transfer.intitialize_two_level_transfers(two_level_transfers);
+ * mg_transfer.build(partitioners);
+ * @endcode
+ *
+ * Alternatively, this class can also be set up as in the case of
+ * MGTransferMatrixFree:
+ * @code
+ * MGTransferGlobalCoarsening mg_transfer;
+ * mg_transfer.initialize_constraints(mg_constrained_dofs);
+ * mg_transfer.build(dof_handler, partitioners);
+ * @endcode
+ * However, this is way to set up is currently only working for globally
+ * refined meshes.
*
* This class currently only works for the tensor-product finite elements
* FE_Q and FE_DGQ and simplex elements FE_SimplexP and FE_SimplexDGP as well as
*/
using Number = typename VectorType::value_type;
+ /**
+ * Default constructor.
+ */
+ MGTransferGlobalCoarsening() = default;
+
+ /**
+ * @name Global coarsening.
+ */
+ /** @{ */
+
/**
* Constructor taking a collection of transfer operators (with the coarsest
* level kept empty in @p transfer) and an optional function that initializes the
const std::function<void(const unsigned int, VectorType &)>
&initialize_dof_vector = {});
+ /**
+ * Set two-level transfers.
+ */
+ template <typename MGTwoLevelTransferObject>
+ void
+ intitialize_two_level_transfers(
+ const MGLevelObject<MGTwoLevelTransferObject> &transfer);
+
/**
* Similar function to MGTransferMatrixFree::build() with the difference that
* the information for the prolongation for each level has already been built.
build(const std::function<void(const unsigned int, VectorType &)>
&initialize_dof_vector);
+ /** @} */
+
+ /**
+ * @name Local smoothing.
+ */
+ /** @{ */
+
+ /**
+ * Constructor with constraints. Equivalent to the default constructor
+ * followed by initialize_constraints().
+ *
+ * @note See also MGTransferMatrixFree.
+ */
+ MGTransferGlobalCoarsening(const MGConstrainedDoFs &mg_constrained_dofs);
+
+ /**
+ * Initialize the constraints to be used in build().
+ *
+ * @note See also MGTransferMatrixFree.
+ */
+ void
+ initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs);
+
+ /**
+ * Actually build the information for the prolongation for each level.
+ *
+ * @note See also MGTransferMatrixFree.
+ */
+ void
+ build(const DoFHandler<dim> &dof_handler,
+ const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
+ &external_partitioners = {});
+
+ /**
+ * Same as above but taking a lambda for initializing vector instead of
+ * partitioners.
+ *
+ * @note See also MGTransferMatrixFree.
+ */
+ void
+ build(const DoFHandler<dim> &dof_handler,
+ const std::function<void(const unsigned int, VectorType &)>
+ &initialize_dof_vector);
+
+ /** @} */
+
+ /**
+ * @name Tranfer functions.
+ */
+ /** @{ */
+
/**
* Perform prolongation.
*/
MGLevelObject<VectorType> &dst,
const InVector & src) const;
+ /** @} */
+
+ /**
+ * @name Utility functions.
+ */
+ /** @{ */
+
/**
* Return the memory consumption of the allocated memory in this class.
*
unsigned int
max_level() const;
+ /** @} */
+
private:
+ /**
+ * Initial internal transfer operator.
+ *
+ * @note See also MGTransferMatrixFree.
+ */
+ void
+ intitialize_internal_transfer(
+ const DoFHandler<dim> & dof_handler,
+ const SmartPointer<const MGConstrainedDoFs> &mg_constrained_dofs);
+
+ /**
+ * Set references to two-level transfer operators to be used.
+ */
+ template <typename MGTwoLevelTransferObject>
+ void
+ intitialize_transfer_references(
+ const MGLevelObject<MGTwoLevelTransferObject> &transfer);
+
/**
* Function to initialize internal level vectors.
*/
void
initialize_dof_vector(const unsigned int level, VectorType &vector) const;
+ /**
+ * MGConstrainedDoFs passed during build().
+ *
+ * @note See also MGTransferMatrixFree.
+ */
+ SmartPointer<const MGConstrainedDoFs> mg_constrained_dofs;
+
+ /**
+ * Internal transfer operator.
+ *
+ * @note See also MGTransferMatrixFree.
+ */
+ MGLevelObject<MGTwoLevelTransfer<dim, VectorType>> internal_transfer;
+
/**
* Collection of the two-level transfer operators.
*/
const MGLevelObject<MGTwoLevelTransferObject> &transfer,
const std::function<void(const unsigned int, VectorType &)>
&initialize_dof_vector)
+{
+ this->intitialize_transfer_references(transfer);
+ this->build(initialize_dof_vector);
+}
+
+
+
+template <int dim, typename VectorType>
+MGTransferGlobalCoarsening<dim, VectorType>::MGTransferGlobalCoarsening(
+ const MGConstrainedDoFs &mg_constrained_dofs)
+{
+ this->initialize_constraints(mg_constrained_dofs);
+}
+
+
+
+template <int dim, typename VectorType>
+template <typename MGTwoLevelTransferObject>
+void
+MGTransferGlobalCoarsening<dim, VectorType>::intitialize_two_level_transfers(
+ const MGLevelObject<MGTwoLevelTransferObject> &transfer)
+{
+ this->intitialize_transfer_references(transfer);
+}
+
+
+
+template <int dim, typename VectorType>
+void
+MGTransferGlobalCoarsening<dim, VectorType>::initialize_constraints(
+ const MGConstrainedDoFs &mg_constrained_dofs)
+{
+ this->mg_constrained_dofs = &mg_constrained_dofs;
+}
+
+
+
+template <int dim, typename VectorType>
+void
+MGTransferGlobalCoarsening<dim, VectorType>::intitialize_internal_transfer(
+ const DoFHandler<dim> & dof_handler,
+ const SmartPointer<const MGConstrainedDoFs> &mg_constrained_dofs)
+{
+ const unsigned int min_level = 0;
+ const unsigned int max_level =
+ dof_handler.get_triangulation().n_global_levels() - 1;
+
+ MGLevelObject<AffineConstraints<typename VectorType::value_type>> constraints(
+ min_level, max_level);
+
+ if (mg_constrained_dofs)
+ for (unsigned int l = min_level; l <= max_level; ++l)
+ mg_constrained_dofs->merge_constraints(
+ constraints[l],
+ l,
+ /*add_boundary_indices*/ true,
+ /*add_refinement_edge_indices*/ false,
+ /*add_level_constraints*/ true,
+ /*add_user_constraints*/ true);
+
+ this->internal_transfer.resize(min_level, max_level);
+
+ for (unsigned int l = min_level; l < max_level; ++l)
+ internal_transfer[l + 1].reinit_geometric_transfer(
+ dof_handler, dof_handler, constraints[l + 1], constraints[l], l + 1, l);
+}
+
+
+
+template <int dim, typename VectorType>
+template <typename MGTwoLevelTransferObject>
+void
+MGTransferGlobalCoarsening<dim, VectorType>::intitialize_transfer_references(
+ const MGLevelObject<MGTwoLevelTransferObject> &transfer)
{
const unsigned int min_level = transfer.min_level();
const unsigned int max_level = transfer.max_level();
this->transfer[l] = &const_cast<MGTwoLevelTransferBase<VectorType> &>(
static_cast<const MGTwoLevelTransferBase<VectorType> &>(
Utilities::get_underlying_value(transfer[l])));
-
- this->build(initialize_dof_vector);
}
+template <int dim, typename VectorType>
+void
+MGTransferGlobalCoarsening<dim, VectorType>::build(
+ const DoFHandler<dim> &dof_handler,
+ const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
+ &external_partitioners)
+{
+ this->intitialize_internal_transfer(dof_handler, mg_constrained_dofs);
+ this->intitialize_transfer_references(internal_transfer);
+ this->build(external_partitioners);
+}
+
+
+
+template <int dim, typename VectorType>
+void
+MGTransferGlobalCoarsening<dim, VectorType>::build(
+ const DoFHandler<dim> &dof_handler,
+ const std::function<void(const unsigned int, VectorType &)>
+ &initialize_dof_vector)
+{
+ this->intitialize_internal_transfer(dof_handler, mg_constrained_dofs);
+ this->intitialize_transfer_references(internal_transfer);
+ this->build(initialize_dof_vector);
+}
+
+
+
template <int dim, typename VectorType>
void
MGTransferGlobalCoarsening<dim, VectorType>::initialize_dof_vector(
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+/**
+ * Test global-coarsening multigrid so that it works also on refinement levels.
+ */
+
+#include "multigrid_util.h"
+
+template <int dim, typename Number = double>
+void
+test(const unsigned int n_refinements, const unsigned int fe_degree_fine)
+{
+ using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+ parallel::distributed::Triangulation<dim> tria(
+ MPI_COMM_WORLD,
+ Triangulation<dim>::MeshSmoothing::none,
+ parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+ GridGenerator::subdivided_hyper_cube(tria, 2);
+ tria.refine_global(n_refinements);
+
+ const unsigned int min_level = 0;
+ const unsigned int max_level = n_refinements;
+
+ MGLevelObject<AffineConstraints<Number>> constraints(min_level, max_level);
+ MGLevelObject<Operator<dim, Number>> operators(min_level, max_level);
+
+ std::unique_ptr<FiniteElement<dim>> fe;
+ std::unique_ptr<Quadrature<dim>> quad;
+ std::unique_ptr<Mapping<dim>> mapping;
+
+ fe = std::make_unique<FE_Q<dim>>(fe_degree_fine);
+ quad = std::make_unique<QGauss<dim>>(fe_degree_fine + 1);
+ mapping = std::make_unique<MappingFE<dim>>(FE_Q<dim>(1));
+
+ // set up dofhandler
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(*fe);
+ dof_handler.distribute_mg_dofs();
+
+ const std::set<types::boundary_id> dirichlet_boundary = {0};
+ MGConstrainedDoFs mg_constrained_dofs;
+ mg_constrained_dofs.initialize(dof_handler);
+ mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,
+ dirichlet_boundary);
+
+ // set up levels
+ for (auto l = min_level; l <= max_level; ++l)
+ {
+ auto &constraint = constraints[l];
+ auto &op = operators[l];
+
+ // set up constraints
+ IndexSet relevant_dofs;
+ DoFTools::extract_locally_relevant_level_dofs(dof_handler,
+ l,
+ relevant_dofs);
+ constraint.reinit(relevant_dofs);
+ constraint.add_lines(mg_constrained_dofs.get_boundary_indices(l));
+ constraint.close();
+
+ constraint.print(std::cout);
+
+ // set up operator
+ op.reinit(*mapping, dof_handler, *quad, constraint, l);
+ }
+
+ // set up transfer operator
+ MGTransferGlobalCoarsening<dim, VectorType> transfer(mg_constrained_dofs);
+
+ transfer.build(dof_handler, [&](const auto l, auto &vec) {
+ operators[l].initialize_dof_vector(vec);
+ });
+
+ GMGParameters mg_data; // TODO
+
+ VectorType dst, src;
+ operators[max_level].initialize_dof_vector(dst);
+ operators[max_level].initialize_dof_vector(src);
+
+ operators[max_level].rhs(src);
+
+ ReductionControl solver_control(
+ mg_data.maxiter, mg_data.abstol, mg_data.reltol, false, false);
+
+ mg_solve(solver_control,
+ dst,
+ src,
+ mg_data,
+ dof_handler,
+ operators[max_level],
+ operators,
+ transfer);
+
+ constraints[max_level].distribute(dst);
+
+ deallog << dim << ' ' << fe_degree_fine << ' ' << n_refinements << ' '
+ << solver_control.last_step() << std::endl;
+
+ return;
+
+ static unsigned int counter = 0;
+
+ MGLevelObject<VectorType> results(min_level, max_level);
+
+ transfer.interpolate_to_mg(dof_handler, results, dst);
+
+ for (unsigned int l = min_level; l <= max_level; ++l)
+ {
+ DataOut<dim> data_out;
+
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(
+ results[l],
+ "solution",
+ DataOut_DoFData<dim, dim>::DataVectorType::type_dof_data);
+ data_out.build_patches(*mapping, 2);
+
+ std::ofstream output("test." + std::to_string(dim) + "." +
+ std::to_string(counter) + "." + std::to_string(l) +
+ ".vtk");
+ data_out.write_vtk(output);
+ }
+
+ counter++;
+}
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+ MPILogInitAll all;
+
+ deallog.precision(8);
+
+ for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
+ for (unsigned int degree = 2; degree <= 4; ++degree)
+ test<2>(n_refinements, degree);
+}
projects / dealii.git / commitdiff