--- /dev/null
+Improved: MGTransferMF can now renumerate DoFs during
+MGTransferMF::copy_to_mg(), MGTransferMF::copy_to_mg(),
+and MGTransferMF::interpolate_to_mg() if the outer
+solver and the multigrid preconditioner have been set up
+with different DoFHandler instances.
+<br>
+(Peter Munch, Laura Prieto Saavedra, 2023/10/03)
template <int dim, int spacedim>
class Base;
}
+
+template <int dim, typename Number>
+class MGTransferMF;
#endif
*/
unsigned int n_components;
+ /**
+ * Pointer to the DoFHandler object used during initialization.
+ */
+ SmartPointer<const DoFHandler<dim>> dof_handler_fine;
+
+ /**
+ * Muligird level used during initialization.
+ */
+ unsigned int mg_level_fine;
+
friend class internal::MGTwoLevelTransferImplementation;
+
+ friend class MGTransferMF<dim, Number>;
};
LinearAlgebra::distributed::Vector<Number> &dst,
const LinearAlgebra::distributed::Vector<Number> &src) const;
+ /**
+ * Pointer to the DoFHandler object used during initialization.
+ */
+ SmartPointer<const DoFHandler<dim>> dof_handler_fine;
+
+ /**
+ * Multigrid level used during initialization.
+ */
+ unsigned int mg_level_fine;
+
/**
* Object to evaluate shape functions on one mesh on visited support points of
* the other mesh.
* point.
*/
std::vector<unsigned int> level_dof_indices_fine_ptrs;
+
+ friend class MGTransferMF<dim, Number>;
};
*
* @note See also MGTransferMatrixFree.
*/
- MGTransferMF() = default;
+ MGTransferMF();
/**
* @name Global coarsening.
void
initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs);
+ /** @} */
+
+ /**
+ * @name Global coarsening and local smoothing.
+ */
+ /** @{ */
+
/**
* Actually build the information for the prolongation for each level.
*
- * @note See also MGTransferMatrixFree.
+ * @note In the case of global coarsening, you can pass into this function
+ * a @p dof_handler with different DoF numbering as the one used within the
+ * provided two-level transfer objects. In this case, vector entries are
+ * permuted during copy_to_mg(), copy_from_mg(), and interpolate_to_mg().
*/
void
build(const DoFHandler<dim> &dof_handler,
* Same as above but taking a lambda for initializing vector instead of
* partitioners.
*
- * @note See also MGTransferMatrixFree.
+ * @note In the case of global coarsening, you can pass into this function
+ * a @p dof_handler with different DoF numbering as the one used within the
+ * provided two-level transfer objects. In this case, vector entries are
+ * permuted during copy_to_mg(), copy_from_mg(), and interpolate_to_mg().
*/
void
build(const DoFHandler<dim> &dof_handler,
*/
template <class InVector>
void
- interpolate_to_mg(MGLevelObject<VectorType> &dst, const InVector &src) const;
+ interpolate_to_mg(const DoFHandler<dim> &dof_handler,
+ MGLevelObject<VectorType> &dst,
+ const InVector &src) const;
/**
- * Like the above function but with a user-provided DoFHandler as
- * additional argument. However, this DoFHandler is not used internally, but
- * is required to be able to use MGTransferMF and
- * MGTransferMatrixFree as template argument.
+ * Interpolate fine-mesh field @p src to each multigrid level in
+ * @p dof_handler and store the result in @p dst.
*/
template <class InVector>
void
- interpolate_to_mg(const DoFHandler<dim> &dof_handler,
- MGLevelObject<VectorType> &dst,
- const InVector &src) const;
+ interpolate_to_mg(MGLevelObject<VectorType> &dst, const InVector &src) const;
/** @} */
const DoFHandler<dim> &dof_handler,
const SmartPointer<const MGConstrainedDoFs> &mg_constrained_dofs);
+ /**
+ * Retreave finest DoFHandler from two-level transfer objects.
+ */
+ std::pair<const DoFHandler<dim> *, unsigned int>
+ get_dof_handler_fine() const;
+
+ /**
+ * Initialize copy indices for MGTransferMF::copy_to_mg(),
+ * MGTransferMF::copy_to_mg(), and MGTransferMF::interpolate_to_mg()
+ * in the case of global coarsening.
+ */
+ void
+ fill_and_communicate_copy_indices_global_coarsening(
+ const DoFHandler<dim> &dof_handler);
+
/**
* Set references to two-level transfer operators to be used.
*/
const InVector &vector_reference,
const bool omit_zeroing_entries) const;
+ /**
+ * Check that the internal DoFHandler is compatible with the external one
+ * provided by copy_to_mg(), copy_from_mg() and interpolate_to_mg()
+ * used, e.g., by PreconditionMG.
+ */
+ void
+ assert_dof_handler(const DoFHandler<dim> &dof_handler_out) const;
+
/**
* Internal transfer operator.
*
const std::function<void(const unsigned int, VectorType &)>
&initialize_dof_vector)
{
+ this->transfer.clear();
+ this->internal_transfer.clear();
+
this->initialize_transfer_references(transfer);
this->build(initialize_dof_vector);
}
MGLevelObject<VectorType> &dst,
const InVector &src) const
{
- (void)dof_handler;
+ assert_dof_handler(dof_handler);
for (unsigned int level = dst.min_level(); level <= dst.max_level(); ++level)
{
OutVector &dst,
const MGLevelObject<VectorType> &src) const
{
- (void)dof_handler;
+ assert_dof_handler(dof_handler);
if (this->perform_plain_copy)
{
MGTransferMF<dim, Number>::interpolate_to_mg(MGLevelObject<VectorType> &dst,
const InVector &src) const
{
+ DoFHandler<dim> dof_handler_dummy;
+
+ this->interpolate_to_mg(dof_handler_dummy, dst, src);
+}
+
+
+
+template <int dim, typename Number>
+template <class InVector>
+void
+MGTransferMF<dim, Number>::interpolate_to_mg(const DoFHandler<dim> &dof_handler,
+ MGLevelObject<VectorType> &dst,
+ const InVector &src) const
+{
+ assert_dof_handler(dof_handler);
+
const unsigned int min_level = transfer.min_level();
const unsigned int max_level = transfer.max_level();
}
}
-
-
-template <int dim, typename Number>
-template <class InVector>
-void
-MGTransferMF<dim, Number>::interpolate_to_mg(const DoFHandler<dim> &dof_handler,
- MGLevelObject<VectorType> &dst,
- const InVector &src) const
-{
- (void)dof_handler;
-
- this->interpolate_to_mg(dst, src);
-}
-
#endif
DEAL_II_NAMESPACE_CLOSE
(mg_level_coarse + 1 == mg_level_fine),
ExcNotImplemented());
- // if (mg_level_fine != numbers::invalid_unsigned_int)
- // AssertIndexRange(mg_level_fine,
- // MGTools::max_level_for_coarse_mesh(
- // dof_handler_fine.get_triangulation()) +
- // 1);
- //
- // if (mg_level_coarse != numbers::invalid_unsigned_int)
- // AssertIndexRange(mg_level_coarse,
- // MGTools::max_level_for_coarse_mesh(
- // dof_handler_coarse.get_triangulation()) +
- // 1);
+ transfer.dof_handler_fine = &dof_handler_fine;
+ transfer.mg_level_fine = mg_level_fine;
std::unique_ptr<FineDoFHandlerViewBase<dim>> dof_handler_fine_view;
"(numbers::invalid_unsigned_int) or on refinement levels without "
"hanging nodes."));
+ transfer.dof_handler_fine = &dof_handler_fine;
+ transfer.mg_level_fine = mg_level_fine;
+
std::unique_ptr<FineDoFHandlerViewBase<dim>> dof_handler_fine_view;
if (internal::p_transfer_involves_repartitioning(dof_handler_fine,
+template <int dim, typename Number>
+MGTransferMF<dim, Number>::MGTransferMF()
+{
+ this->transfer.clear();
+ this->internal_transfer.clear();
+}
+
+
+
template <int dim, typename Number>
MGTransferMF<dim, Number>::MGTransferMF(
const MGConstrainedDoFs &mg_constrained_dofs)
{
+ this->transfer.clear();
+ this->internal_transfer.clear();
this->initialize_constraints(mg_constrained_dofs);
}
+template <int dim, typename Number>
+std::pair<const DoFHandler<dim> *, unsigned int>
+MGTransferMF<dim, Number>::get_dof_handler_fine() const
+{
+ if (this->transfer.n_levels() <=
+ 1) // single level: the information cannot be retreaved
+ return {nullptr, numbers::invalid_unsigned_int};
+
+ if (const auto t = dynamic_cast<
+ const MGTwoLevelTransfer<dim,
+ LinearAlgebra::distributed::Vector<Number>> *>(
+ this->transfer[this->transfer.max_level()].get()))
+ {
+ return {t->dof_handler_fine, t->mg_level_fine};
+ }
+ else if (const auto t = dynamic_cast<const MGTwoLevelTransferNonNested<
+ dim,
+ LinearAlgebra::distributed::Vector<Number>> *>(
+ this->transfer[this->transfer.max_level()].get()))
+ {
+ return {t->dof_handler_fine, t->mg_level_fine};
+ }
+ else
+ {
+ Assert(false, ExcNotImplemented());
+ return {nullptr, numbers::invalid_unsigned_int};
+ }
+}
+
+
+
+template <int dim, typename Number>
+void
+MGTransferMF<dim, Number>::fill_and_communicate_copy_indices_global_coarsening(
+ const DoFHandler<dim> &dof_handler_out)
+{
+ const auto [dof_handler_in, level_in] = get_dof_handler_fine();
+
+ if ((dof_handler_in == nullptr) || (dof_handler_in == &dof_handler_out))
+ return; // nothing to do
+
+ this->copy_indices.resize(1);
+ this->copy_indices[0].reinit(2, dof_handler_out.n_locally_owned_dofs());
+
+ std::vector<types::global_dof_index> dof_indices_in;
+ std::vector<types::global_dof_index> dof_indices_out;
+
+ this->perform_plain_copy = true;
+
+ const auto &is_out = (level_in == numbers::invalid_unsigned_int) ?
+ dof_handler_out.locally_owned_dofs() :
+ dof_handler_out.locally_owned_mg_dofs(level_in);
+
+ const auto &is_in = (level_in == numbers::invalid_unsigned_int) ?
+ dof_handler_in->locally_owned_dofs() :
+ dof_handler_in->locally_owned_mg_dofs(level_in);
+
+ internal::loop_over_active_or_level_cells(
+ dof_handler_in->get_triangulation(), level_in, [&](const auto &cell) {
+ const auto cell_id = cell->id();
+
+ Assert(
+ dof_handler_out.get_triangulation().contains_cell(cell_id),
+ ExcMessage(
+ "DoFHandler instances used for set up of MGTransferMF and copy_to_mg(), "
+ "copy_from_mg(), or interpolate_to_mg() are not compatible."));
+
+ if (level_in == numbers::invalid_unsigned_int)
+ {
+ const auto cell_in = cell->as_dof_handler_iterator(*dof_handler_in);
+ const auto cell_out = dof_handler_out.get_triangulation()
+ .create_cell_iterator(cell_id)
+ ->as_dof_handler_iterator(dof_handler_out);
+
+ AssertDimension(cell_in->get_fe().n_dofs_per_cell(),
+ cell_out->get_fe().n_dofs_per_cell());
+
+ dof_indices_in.resize(cell_in->get_fe().n_dofs_per_cell());
+ dof_indices_out.resize(cell_out->get_fe().n_dofs_per_cell());
+
+ cell_in->get_dof_indices(dof_indices_in);
+ cell_out->get_dof_indices(dof_indices_out);
+ }
+ else
+ {
+ const auto cell_in =
+ cell->as_dof_handler_level_iterator(*dof_handler_in);
+ const auto cell_out =
+ dof_handler_out.get_triangulation()
+ .create_cell_iterator(cell_id)
+ ->as_dof_handler_level_iterator(dof_handler_out);
+
+ AssertDimension(cell_in->get_fe().n_dofs_per_cell(),
+ cell_out->get_fe().n_dofs_per_cell());
+
+ dof_indices_in.resize(cell_in->get_fe().n_dofs_per_cell());
+ dof_indices_out.resize(cell_out->get_fe().n_dofs_per_cell());
+
+ cell_in->get_mg_dof_indices(dof_indices_in);
+ cell_out->get_mg_dof_indices(dof_indices_out);
+ }
+
+ this->perform_plain_copy &= (dof_indices_in == dof_indices_out);
+
+ for (unsigned int i = 0; i < dof_indices_in.size(); ++i)
+ if (is_out.is_element(dof_indices_out[i]))
+ this->copy_indices[0](1,
+ is_out.index_within_set(dof_indices_out[i])) =
+ is_in.index_within_set(dof_indices_in[i]);
+ });
+
+
+ this->perform_plain_copy =
+ Utilities::MPI::max(this->perform_plain_copy ? 1 : 0,
+ dof_handler_out.get_communicator()) != 0;
+
+ if (this->perform_plain_copy)
+ {
+ this->copy_indices.clear();
+ }
+ else
+ {
+ this->perform_renumbered_plain_copy = true;
+ this->solution_copy_indices = this->copy_indices;
+ }
+}
+
+
+
template <int dim, typename Number>
void
MGTransferMF<dim, Number>::build(
const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
&external_partitioners)
{
- this->initialize_internal_transfer(dof_handler, this->mg_constrained_dofs);
- this->initialize_transfer_references(internal_transfer);
+ const bool use_local_smoothing =
+ this->transfer.n_levels() == 0 || this->internal_transfer.n_levels() > 0;
+
+ if (use_local_smoothing)
+ {
+ this->initialize_internal_transfer(dof_handler,
+ this->mg_constrained_dofs);
+ this->initialize_transfer_references(internal_transfer);
+ }
+
this->build(external_partitioners);
- this->fill_and_communicate_copy_indices(dof_handler);
+
+ if (use_local_smoothing)
+ this->fill_and_communicate_copy_indices(dof_handler);
+ else
+ this->fill_and_communicate_copy_indices_global_coarsening(dof_handler);
}
const std::function<void(const unsigned int, VectorType &)>
&initialize_dof_vector)
{
- this->initialize_internal_transfer(dof_handler, this->mg_constrained_dofs);
- this->initialize_transfer_references(internal_transfer);
+ const bool use_local_smoothing =
+ this->transfer.n_levels() == 0 || this->internal_transfer.n_levels() > 0;
+
+ if (use_local_smoothing)
+ {
+ this->initialize_internal_transfer(dof_handler,
+ this->mg_constrained_dofs);
+ this->initialize_transfer_references(internal_transfer);
+ }
+
this->build(initialize_dof_vector);
- this->fill_and_communicate_copy_indices(dof_handler);
+
+ if (use_local_smoothing)
+ this->fill_and_communicate_copy_indices(dof_handler);
+ else
+ this->fill_and_communicate_copy_indices_global_coarsening(dof_handler);
}
+template <int dim, typename Number>
+void
+MGTransferMF<dim, Number>::assert_dof_handler(
+ const DoFHandler<dim> &dof_handler_out) const
+{
+#ifndef DEBUG
+ (void)dof_handler_out;
+#else
+
+ const auto [dof_handler_in, level_in] = get_dof_handler_fine();
+
+ if ((dof_handler_out.n_dofs() == 0) || // dummy DoFHandler
+ (dof_handler_in == nullptr) || // single level
+ (dof_handler_in == &dof_handler_out // same DoFHandler
+ ))
+ return; // nothing to do
+
+ if (this->perform_plain_copy)
+ {
+ // global-coarsening path: compare indices of cells
+
+ std::vector<types::global_dof_index> dof_indices_in;
+ std::vector<types::global_dof_index> dof_indices_out;
+
+ internal::loop_over_active_or_level_cells(
+ dof_handler_in->get_triangulation(), level_in, [&](const auto &cell) {
+ const auto cell_id = cell->id();
+
+ Assert(
+ dof_handler_out.get_triangulation().contains_cell(cell_id),
+ ExcMessage(
+ "DoFHandler instances used for set up of MGTransferMF and copy_to_mg(), "
+ "copy_from_mg(), or interpolate_to_mg() are not compatible."));
+
+ if (level_in == numbers::invalid_unsigned_int)
+ {
+ const auto cell_in =
+ cell->as_dof_handler_iterator(*dof_handler_in);
+ const auto cell_out =
+ dof_handler_out.get_triangulation()
+ .create_cell_iterator(cell_id)
+ ->as_dof_handler_iterator(dof_handler_out);
+
+ AssertDimension(cell_in->get_fe().n_dofs_per_cell(),
+ cell_out->get_fe().n_dofs_per_cell());
+
+ dof_indices_in.resize(cell_in->get_fe().n_dofs_per_cell());
+ dof_indices_out.resize(cell_out->get_fe().n_dofs_per_cell());
+
+ cell_in->get_dof_indices(dof_indices_in);
+ cell_out->get_dof_indices(dof_indices_out);
+ }
+ else
+ {
+ const auto cell_in =
+ cell->as_dof_handler_level_iterator(*dof_handler_in);
+ const auto cell_out =
+ dof_handler_out.get_triangulation()
+ .create_cell_iterator(cell_id)
+ ->as_dof_handler_level_iterator(dof_handler_out);
+
+ AssertDimension(cell_in->get_fe().n_dofs_per_cell(),
+ cell_out->get_fe().n_dofs_per_cell());
+
+ dof_indices_in.resize(cell_in->get_fe().n_dofs_per_cell());
+ dof_indices_out.resize(cell_out->get_fe().n_dofs_per_cell());
+
+ cell_in->get_mg_dof_indices(dof_indices_in);
+ cell_out->get_mg_dof_indices(dof_indices_out);
+ }
+
+ Assert(
+ dof_indices_in == dof_indices_out,
+ ExcMessage(
+ "DoFHandler instances used for set up of MGTransferMF and copy_to_mg(), "
+ "copy_from_mg(), or interpolate_to_mg() are not compatible."));
+ });
+ }
+ else if (this->perform_renumbered_plain_copy)
+ {
+ // nothing to do
+ }
+#endif
+}
+
+
+
template <int dim, typename Number>
std::size_t
MGTransferMF<dim, Number>::memory_consumption() const
dof_handler_fine.get_fe().n_components() > 0,
ExcNotImplemented());
+ this->dof_handler_fine = &dof_handler_fine;
+ this->mg_level_fine = numbers::invalid_unsigned_int;
+
this->fine_element_is_continuous =
dof_handler_fine.get_fe().n_dofs_per_vertex() > 0;
--- /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.
+//
+// ---------------------------------------------------------------------
+
+
+/**
+ * Similar to test multigrid_01 but the DoFs are renumbered in the DoFHandler
+ * used by the outer solver.
+ */
+
+#include <deal.II/dofs/dof_renumbering.h>
+
+#include "multigrid_util.h"
+
+template <int dim, typename Number = double>
+void
+test(const unsigned int n_refinements,
+ const unsigned int fe_degree_fine,
+ const bool do_simplex_mesh,
+ const unsigned int mesh_type)
+{
+ using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+ Triangulation<dim> tria;
+ if (do_simplex_mesh)
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, 2);
+ else
+ GridGenerator::subdivided_hyper_cube(tria, 2);
+
+ if (mesh_type == 0)
+ {
+ tria.refine_global(n_refinements);
+ }
+ else if (mesh_type == 1)
+ {
+ for (unsigned int i = 1; i < n_refinements; ++i)
+ {
+ for (auto cell : tria.active_cell_iterators())
+ if (cell->is_locally_owned())
+ {
+ bool flag = true;
+ for (int d = 0; d < dim; ++d)
+ if (cell->center()[d] > 0.5)
+ flag = false;
+ if (flag)
+ cell->set_refine_flag();
+ }
+ tria.execute_coarsening_and_refinement();
+ }
+ }
+ else
+ AssertThrow(false, ExcNotImplemented());
+
+ ////////////////////////////////////////////////////////////////////////
+
+ std::unique_ptr<FiniteElement<dim>> fe;
+ std::unique_ptr<Quadrature<dim>> quad;
+ std::unique_ptr<Mapping<dim>> mapping;
+
+ if (do_simplex_mesh)
+ {
+ fe = std::make_unique<FE_SimplexP<dim>>(fe_degree_fine);
+ quad = std::make_unique<QGaussSimplex<dim>>(fe_degree_fine + 1);
+ mapping = std::make_unique<MappingFE<dim>>(FE_SimplexP<dim>(1));
+ }
+ else
+ {
+ 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));
+ }
+
+ DoFHandler<dim> fine_dof_handler(tria);
+ fine_dof_handler.distribute_dofs(*fe);
+ DoFRenumbering::Cuthill_McKee(fine_dof_handler);
+
+ AffineConstraints<Number> fine_constraints;
+ const IndexSet locally_relevant_dofs =
+ DoFTools::extract_locally_relevant_dofs(fine_dof_handler);
+ fine_constraints.reinit(locally_relevant_dofs);
+ VectorTools::interpolate_boundary_values(*mapping,
+ fine_dof_handler,
+ 0,
+ Functions::ZeroFunction<dim>(),
+ fine_constraints);
+ DoFTools::make_hanging_node_constraints(fine_dof_handler, fine_constraints);
+ fine_constraints.close();
+
+ // set up operator
+ Operator<dim, Number> fine_operator;
+ fine_operator.reinit(*mapping, fine_dof_handler, *quad, fine_constraints);
+
+ ////////////////////////////////////////////////////////////////////////
+
+ const auto level_degrees =
+ MGTransferGlobalCoarseningTools::create_polynomial_coarsening_sequence(
+ fe_degree_fine,
+ MGTransferGlobalCoarseningTools::PolynomialCoarseningSequenceType::
+ bisect);
+
+ const unsigned int min_level = 0;
+ const unsigned int max_level = level_degrees.size() - 1;
+
+ MGLevelObject<DoFHandler<dim>> dof_handlers(min_level, max_level, tria);
+ MGLevelObject<AffineConstraints<Number>> constraints(min_level, max_level);
+ MGLevelObject<MGTwoLevelTransfer<dim, VectorType>> transfers(min_level,
+ max_level);
+ MGLevelObject<Operator<dim, Number>> operators(min_level, max_level);
+
+ std::unique_ptr<Mapping<dim>> mapping_;
+
+ // set up levels
+ for (auto l = min_level; l <= max_level; ++l)
+ {
+ auto &dof_handler = dof_handlers[l];
+ auto &constraint = constraints[l];
+ auto &op = operators[l];
+
+ std::unique_ptr<FiniteElement<dim>> fe;
+ std::unique_ptr<Quadrature<dim>> quad;
+ std::unique_ptr<Mapping<dim>> mapping;
+
+ if (do_simplex_mesh)
+ {
+ fe = std::make_unique<FE_SimplexP<dim>>(level_degrees[l]);
+ quad = std::make_unique<QGaussSimplex<dim>>(level_degrees[l] + 1);
+ mapping = std::make_unique<MappingFE<dim>>(FE_SimplexP<dim>(1));
+ }
+ else
+ {
+ fe = std::make_unique<FE_Q<dim>>(level_degrees[l]);
+ quad = std::make_unique<QGauss<dim>>(level_degrees[l] + 1);
+ mapping = std::make_unique<MappingFE<dim>>(FE_Q<dim>(1));
+ }
+
+ if (l == max_level)
+ mapping_ = mapping->clone();
+
+ // set up dofhandler
+ dof_handler.distribute_dofs(*fe);
+
+ // set up constraints
+ const IndexSet locally_relevant_dofs =
+ DoFTools::extract_locally_relevant_dofs(dof_handler);
+ constraint.reinit(locally_relevant_dofs);
+ VectorTools::interpolate_boundary_values(
+ *mapping, dof_handler, 0, Functions::ZeroFunction<dim>(), constraint);
+ DoFTools::make_hanging_node_constraints(dof_handler, constraint);
+ constraint.close();
+
+ // set up operator
+ op.reinit(*mapping, dof_handler, *quad, constraint);
+ }
+
+ // set up transfer operator
+ for (unsigned int l = min_level; l < max_level; ++l)
+ transfers[l + 1].reinit(dof_handlers[l + 1],
+ dof_handlers[l],
+ constraints[l + 1],
+ constraints[l]);
+
+ MGTransferGlobalCoarsening<dim, VectorType> transfer;
+ transfer.initialize_two_level_transfers(transfers);
+ transfer.build(fine_dof_handler, [&](const auto l, auto &vec) {
+ operators[l].initialize_dof_vector(vec);
+ });
+
+ GMGParameters mg_data; // TODO
+
+ VectorType dst, src;
+ fine_operator.initialize_dof_vector(dst);
+ fine_operator.initialize_dof_vector(src);
+ fine_operator.rhs(src);
+
+ ReductionControl solver_control(
+ mg_data.maxiter, mg_data.abstol, mg_data.reltol, false, false);
+
+ mg_solve(solver_control,
+ dst,
+ src,
+ mg_data,
+ fine_dof_handler,
+ fine_operator,
+ operators,
+ transfer);
+
+ fine_constraints.distribute(dst);
+
+ deallog << dim << ' ' << fe_degree_fine << ' ' << n_refinements << ' '
+ << (do_simplex_mesh ? "tri " : "quad") << ' '
+ << solver_control.last_step() << std::endl;
+
+ static unsigned int counter = 0;
+
+ MGLevelObject<VectorType> results(min_level, max_level);
+
+ transfer.interpolate_to_mg(dof_handlers[max_level], results, dst);
+
+ for (unsigned int l = min_level; l <= max_level; ++l)
+ {
+ deallog << "Norm interpolated solution on level " << l << ": "
+ << results[l].l2_norm() << std::endl;
+ DataOut<dim> data_out;
+
+ data_out.attach_dof_handler(dof_handlers[l]);
+ 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 mesh_type = 0; mesh_type < 2; ++mesh_type)
+ {
+ for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
+ for (unsigned int degree = 2; degree <= 4; ++degree)
+ test<2>(n_refinements, degree, false /*quadrilateral*/, mesh_type);
+
+ for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
+ for (unsigned int degree = 2; degree <= 2; ++degree)
+ test<2>(n_refinements, degree, true /*triangle*/, mesh_type);
+ }
+}
--- /dev/null
+
+DEAL:0::2 2 2 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.33877875
+DEAL:0::Norm interpolated solution on level 1: 0.66014529
+DEAL:0::2 3 2 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.33872299
+DEAL:0::Norm interpolated solution on level 1: 0.99016130
+DEAL:0::2 4 2 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.33877890
+DEAL:0::Norm interpolated solution on level 1: 0.66010393
+DEAL:0::Norm interpolated solution on level 2: 1.3201957
+DEAL:0::2 2 3 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.66459543
+DEAL:0::Norm interpolated solution on level 1: 1.3203627
+DEAL:0::2 3 3 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.66458798
+DEAL:0::Norm interpolated solution on level 1: 1.9805368
+DEAL:0::2 4 3 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.66459725
+DEAL:0::Norm interpolated solution on level 1: 1.3203574
+DEAL:0::Norm interpolated solution on level 2: 2.6407133
+DEAL:0::2 2 4 quad 3
+DEAL:0::Norm interpolated solution on level 0: 1.3225826
+DEAL:0::Norm interpolated solution on level 1: 2.6407347
+DEAL:0::2 3 4 quad 3
+DEAL:0::Norm interpolated solution on level 0: 1.3225816
+DEAL:0::Norm interpolated solution on level 1: 3.9611011
+DEAL:0::2 4 4 quad 3
+DEAL:0::Norm interpolated solution on level 0: 1.3225830
+DEAL:0::Norm interpolated solution on level 1: 2.6407340
+DEAL:0::Norm interpolated solution on level 2: 5.2814679
+DEAL:0::2 2 2 tri 3
+DEAL:0::Norm interpolated solution on level 0: 0.33806168
+DEAL:0::Norm interpolated solution on level 1: 0.66011744
+DEAL:0::2 2 3 tri 3
+DEAL:0::Norm interpolated solution on level 0: 0.66445481
+DEAL:0::Norm interpolated solution on level 1: 1.3203595
+DEAL:0::2 2 4 tri 3
+DEAL:0::Norm interpolated solution on level 0: 1.3225585
+DEAL:0::Norm interpolated solution on level 1: 2.6407344
+DEAL:0::2 2 2 quad 2
+DEAL:0::Norm interpolated solution on level 0: 0.13157424
+DEAL:0::Norm interpolated solution on level 1: 0.24449632
+DEAL:0::2 3 2 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.13144458
+DEAL:0::Norm interpolated solution on level 1: 0.35587970
+DEAL:0::2 4 2 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.13139421
+DEAL:0::Norm interpolated solution on level 1: 0.24404614
+DEAL:0::Norm interpolated solution on level 2: 0.46558561
+DEAL:0::2 2 3 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.23622379
+DEAL:0::Norm interpolated solution on level 1: 0.45479680
+DEAL:0::2 3 3 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.23657515
+DEAL:0::Norm interpolated solution on level 1: 0.66773399
+DEAL:0::2 4 3 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.23652211
+DEAL:0::Norm interpolated solution on level 1: 0.45456995
+DEAL:0::Norm interpolated solution on level 2: 0.87733284
+DEAL:0::2 2 4 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.43119420
+DEAL:0::Norm interpolated solution on level 1: 0.85199621
+DEAL:0::2 3 4 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.43119824
+DEAL:0::Norm interpolated solution on level 1: 1.2615628
+DEAL:0::2 4 4 quad 3
+DEAL:0::Norm interpolated solution on level 0: 0.43118187
+DEAL:0::Norm interpolated solution on level 1: 0.85199072
+DEAL:0::Norm interpolated solution on level 2: 1.6669458
+DEAL:0::2 2 2 tri 3
+DEAL:0::Norm interpolated solution on level 0: 0.13575160
+DEAL:0::Norm interpolated solution on level 1: 0.24280231
+DEAL:0::2 2 3 tri 3
+DEAL:0::Norm interpolated solution on level 0: 0.21915908
+DEAL:0::Norm interpolated solution on level 1: 0.41924192
+DEAL:0::2 2 4 tri 4
+DEAL:0::Norm interpolated solution on level 0: 0.40718520
+DEAL:0::Norm interpolated solution on level 1: 0.80533381
projects / dealii.git / commitdiff