Skip to content

Commit

Permalink
MGTwoLevelTransferNonNested: simplices with multiple components
Browse files Browse the repository at this point in the history
  • Loading branch information
peterrum committed Sep 27, 2023
1 parent 368e863 commit c66425a
Show file tree
Hide file tree
Showing 16 changed files with 269 additions and 71 deletions.
54 changes: 33 additions & 21 deletions include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
Original file line number Diff line number Diff line change
Expand Up @@ -4472,18 +4472,20 @@ namespace internal
ExcMessage("Function expects FE_DGQ, FE_Q, FE_SimplexP, or "
"FE_SimplexDGP in dof_handler."));

Assert((dynamic_cast<const FE_Q<dim, spacedim> *>(
&dof_handler_support_points.get_fe()) != nullptr ||
dynamic_cast<const FE_SimplexP<dim, spacedim> *>(
&dof_handler_support_points.get_fe()) != nullptr) ||
((dynamic_cast<const FE_DGQ<dim, spacedim> *>(
&dof_handler_support_points.get_fe()) != nullptr ||
dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(
&dof_handler_support_points.get_fe()) != nullptr) &&
dof_handler_support_points.get_fe().degree == 0),
ExcMessage(
"Function expects (FE_DGQ||FE_SimplexDGP)&&degree==0 or "
"(FE_Q||FE_SimplexP) in dof_handler_support_points."));
Assert(
(dynamic_cast<const FE_Q<dim, spacedim> *>(
&dof_handler_support_points.get_fe().base_element(0)) != nullptr ||
dynamic_cast<const FE_SimplexP<dim, spacedim> *>(
&dof_handler_support_points.get_fe().base_element(0)) != nullptr) ||
((dynamic_cast<const FE_DGQ<dim, spacedim> *>(
&dof_handler_support_points.get_fe().base_element(0)) !=
nullptr ||
dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(
&dof_handler_support_points.get_fe().base_element(0)) !=
nullptr) &&
dof_handler_support_points.get_fe().degree == 0),
ExcMessage("Function expects (FE_DGQ||FE_SimplexDGP)&&degree==0 or "
"(FE_Q||FE_SimplexP) in dof_handler_support_points."));

Assert(
dof_handler_support_points.get_fe().n_components() == 1,
Expand Down Expand Up @@ -4520,7 +4522,8 @@ namespace internal
const bool needs_conversion =
dof_handler.get_fe().conforming_space ==
FiniteElementData<dim>::Conformity::L2 &&
(dof_handler.get_fe().degree > 0);
(dof_handler.get_fe().degree > 0) &&
dof_handler.get_fe().reference_cell().is_hyper_cube();
std::vector<unsigned int> lexicographic_to_hierarchic;
if (needs_conversion)
lexicographic_to_hierarchic =
Expand Down Expand Up @@ -4658,8 +4661,10 @@ namespace internal
const auto n_components = fe.n_components();

if (n_components == 1 &&
(fe.conforming_space == FiniteElementData<dim>::Conformity::H1 ||
degree == 0))
((fe.reference_cell().is_hyper_cube() ||
fe.reference_cell().is_simplex()) &&
(fe.conforming_space == FiniteElementData<dim>::Conformity::H1 ||
degree == 0)))
{
// in case a DG space of order 0 is provided, DoFs indices are always
// uniquely assigned to support points (they are always defined in the
Expand All @@ -4678,7 +4683,13 @@ namespace internal
auto dof_handler_support_points =
std::make_shared<DoFHandler<dim, spacedim>>(tria);

if (degree == 0)
if (fe.reference_cell().is_simplex() && (degree == 0))
dof_handler_support_points->distribute_dofs(
FE_SimplexDGP<dim, spacedim>(degree));
else if (fe.reference_cell().is_simplex())
dof_handler_support_points->distribute_dofs(
FE_SimplexP<dim, spacedim>(degree));
else if (degree == 0)
dof_handler_support_points->distribute_dofs(
FE_DGQ<dim, spacedim>(degree));
else
Expand Down Expand Up @@ -4893,12 +4904,13 @@ MGTwoLevelTransferNonNested<dim, LinearAlgebra::distributed::Vector<Number>>::
if (const auto fe = dynamic_cast<const FE_Q<dim> *>(&fe_base))
fe_coarse = std::make_unique<FESystem<dim>>(FE_DGQ<dim>(fe->get_degree()),
n_components);
else if (const auto fe = dynamic_cast<const FE_DGQ<dim> *>(&fe_base))
fe_coarse = fe->clone();
else if (const auto fe = dynamic_cast<const FE_SimplexP<dim> *>(&fe_base))
fe_coarse = fe->clone();
else if (const auto fe = dynamic_cast<const FE_SimplexDGP<dim> *>(&fe_base))
fe_coarse = fe->clone();
fe_coarse =
std::make_unique<FESystem<dim>>(FE_SimplexDGP<dim>(fe->get_degree()),
n_components);
else if (dynamic_cast<const FE_DGQ<dim> *>(&fe_base) ||
dynamic_cast<const FE_SimplexP<dim> *>(&fe_base))
fe_coarse = dof_handler_coarse.get_fe().clone();
else
AssertThrow(false, ExcMessage(dof_handler_coarse.get_fe().get_name()));
}
Expand Down
2 changes: 1 addition & 1 deletion tests/multigrid-global-coarsening/multigrid_a_01.cc
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ test(const unsigned int n_refinements,
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);
MGLevelObject<Operator<dim, 1, Number>> operators(min_level, max_level);

std::unique_ptr<Mapping<dim>> mapping_;

Expand Down
2 changes: 1 addition & 1 deletion tests/multigrid-global-coarsening/multigrid_a_02.cc
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ test(const unsigned int n_refinements, const unsigned int fe_degree_fine)
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);
MGLevelObject<Operator<dim, 1, Number>> operators(min_level, max_level);

std::unique_ptr<FiniteElement<dim>> fe;
std::unique_ptr<Quadrature<dim>> quad;
Expand Down
2 changes: 1 addition & 1 deletion tests/multigrid-global-coarsening/multigrid_a_03.cc
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ test(const unsigned int n_refinements, const unsigned int fe_degree_fine)
const unsigned int max_level = n_refinements;

MGLevelObject<AffineConstraints<Number>> constraints(min_level, max_level);
MGLevelObject<Operator<dim, Number>> operators(min_level, max_level);
MGLevelObject<Operator<dim, 1, Number>> operators(min_level, max_level);

std::unique_ptr<FiniteElement<dim>> fe;
std::unique_ptr<Quadrature<dim>> quad;
Expand Down
2 changes: 1 addition & 1 deletion tests/multigrid-global-coarsening/multigrid_p_01.cc
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ test(const unsigned int n_refinements,
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);
MGLevelObject<Operator<dim, 1, Number>> operators(min_level, max_level);

std::unique_ptr<Mapping<dim>> mapping_;

Expand Down
2 changes: 1 addition & 1 deletion tests/multigrid-global-coarsening/multigrid_p_02.cc
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ test(const unsigned int n_refinements, const unsigned int fe_degree_fine)
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);
MGLevelObject<Operator<dim, 1, Number>> operators(min_level, max_level);

std::unique_ptr<Mapping<dim>> mapping_;

Expand Down
2 changes: 1 addition & 1 deletion tests/multigrid-global-coarsening/multigrid_p_03.cc
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@ test(const unsigned int n_refinements,
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);
MGLevelObject<Operator<dim, 1, Number>> operators(min_level, max_level);

std::unique_ptr<Mapping<dim>> mapping_;

Expand Down
23 changes: 20 additions & 3 deletions tests/multigrid-global-coarsening/multigrid_util.h
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@

#include "../tests.h"

template <int dim_, typename Number>
template <int dim_, int n_components = dim_, typename Number = double>
class Operator : public Subscriptor
{
public:
Expand All @@ -68,7 +68,7 @@ class Operator : public Subscriptor

static const int dim = dim_;

using FECellIntegrator = FEEvaluation<dim, -1, 0, 1, number>;
using FECellIntegrator = FEEvaluation<dim, -1, 0, n_components, Number>;

void
reinit(const Mapping<dim> &mapping,
Expand Down Expand Up @@ -195,7 +195,24 @@ class Operator : public Subscriptor
{
phi.reinit(cell);
for (unsigned int q = 0; q < phi.n_q_points; ++q)
phi.submit_value(1.0, q);
{
if constexpr (n_components == 1)
{
phi.submit_value(1.0, q);
}
else
{
Tensor<1, n_components, VectorizedArray<Number>> temp;
for (unsigned int v = 0;
v < VectorizedArray<Number>::size();
++v)
{
for (unsigned int i = 0; i < n_components; i++)
temp[i][v] = 1.;
}
phi.submit_value(temp, q);
}
}

phi.integrate_scatter(EvaluationFlags::values, dst);
}
Expand Down
4 changes: 2 additions & 2 deletions tests/multigrid-global-coarsening/non_nested_multigrid_01.cc
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,8 @@ test(const unsigned int n_refinements,
MGLevelObject<AffineConstraints<Number>> constraints(min_level, max_level);
MGLevelObject<std::unique_ptr<Mapping<dim>>> mappings(min_level, max_level);
MGLevelObject<std::shared_ptr<MGTwoLevelTransferNonNested<dim, VectorType>>>
transfers(min_level, max_level);
MGLevelObject<Operator<dim, Number>> operators(min_level, max_level);
transfers(min_level, max_level);
MGLevelObject<Operator<dim, 1, Number>> operators(min_level, max_level);


// set up levels
Expand Down
4 changes: 2 additions & 2 deletions tests/multigrid-global-coarsening/non_nested_multigrid_02.cc
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,8 @@ test(const unsigned int n_refinements,
MGLevelObject<AffineConstraints<Number>> constraints(min_level, max_level);
MGLevelObject<MappingQ1<dim>> mappings(min_level, max_level);
MGLevelObject<std::shared_ptr<MGTwoLevelTransferNonNested<dim, VectorType>>>
transfers(min_level, max_level);
MGLevelObject<Operator<dim, Number>> operators(min_level, max_level);
transfers(min_level, max_level);
MGLevelObject<Operator<dim, 1, Number>> operators(min_level, max_level);


// set up levels
Expand Down
4 changes: 2 additions & 2 deletions tests/multigrid-global-coarsening/non_nested_multigrid_03.cc
Original file line number Diff line number Diff line change
Expand Up @@ -37,8 +37,8 @@ test(const unsigned int n_refinements, const unsigned int fe_degree_fine)
MGLevelObject<AffineConstraints<Number>> constraints(min_level, max_level);
MGLevelObject<MappingQ1<dim>> mappings(min_level, max_level);
MGLevelObject<std::shared_ptr<MGTwoLevelTransferNonNested<dim, VectorType>>>
transfers(min_level, max_level);
MGLevelObject<Operator<dim, Number>> operators(min_level, max_level);
transfers(min_level, max_level);
MGLevelObject<Operator<dim, 1, Number>> operators(min_level, max_level);


// set up levels
Expand Down
4 changes: 2 additions & 2 deletions tests/multigrid-global-coarsening/non_nested_multigrid_04.cc
Original file line number Diff line number Diff line change
Expand Up @@ -97,8 +97,8 @@ test(const unsigned int n_refinements,
MGLevelObject<AffineConstraints<Number>> constraints(min_level, max_level);
MGLevelObject<std::unique_ptr<Mapping<dim>>> mappings(min_level, max_level);
MGLevelObject<std::shared_ptr<MGTwoLevelTransferNonNested<dim, VectorType>>>
transfers(min_level, max_level);
MGLevelObject<Operator<dim, Number>> operators(min_level, max_level);
transfers(min_level, max_level);
MGLevelObject<Operator<dim, 1, Number>> operators(min_level, max_level);


// set up levels
Expand Down
160 changes: 160 additions & 0 deletions tests/multigrid-global-coarsening/non_nested_multigrid_05.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,160 @@
// // ---------------------------------------------------------------------
// //
// // Copyright (C) 2020 - 2023 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 for a uniformly refined mesh both for
* simplicial elements and multiple components.
*/

#include <deal.II/fe/fe_system.h>

#include <deal.II/grid/grid_tools.h>

#include "multigrid_util.h"

template <int dim, int n_components = 1, typename Number = double>
void
test(const unsigned int n_refinements,
const unsigned int fe_degree_fine,
const double factor)
{
using VectorType = LinearAlgebra::distributed::Vector<Number>;

const unsigned int min_level = 0;
const unsigned int max_level = n_refinements;

MGLevelObject<Triangulation<dim>> triangulations(min_level, max_level);
MGLevelObject<DoFHandler<dim>> dof_handlers(min_level, max_level);
MGLevelObject<AffineConstraints<Number>> constraints(min_level, max_level);
MGLevelObject<std::unique_ptr<Mapping<dim>>> mappings(min_level, max_level);
MGLevelObject<std::shared_ptr<MGTwoLevelTransferNonNested<dim, VectorType>>>
transfers(min_level, max_level);
MGLevelObject<Operator<dim, n_components, Number>> operators(min_level,
max_level);


// set up levels
for (auto l = min_level; l <= max_level; ++l)
{
auto &tria = triangulations[l];
auto &dof_handler = dof_handlers[l];
auto &constraint = constraints[l];
auto &op = operators[l];

std::unique_ptr<FESystem<dim>> fe =
std::make_unique<FESystem<dim>>(FE_SimplexP<dim>(fe_degree_fine),
n_components);
std::unique_ptr<Quadrature<dim>> quad =
std::make_unique<QGaussSimplex<dim>>(fe_degree_fine + 1);
std::unique_ptr<Mapping<dim>> _mapping =
std::make_unique<MappingFE<dim>>(FE_SimplexP<dim>(fe_degree_fine));

mappings[l] = _mapping->clone();

// set up triangulation
Triangulation<dim> tria_tmp;
GridGenerator::hyper_cube(tria_tmp, -1.0, 1.0);
tria_tmp.refine_global(l);

GridGenerator::convert_hypercube_to_simplex_mesh(tria_tmp, tria);
GridTools::distort_random(factor,
tria,
true,
boost::random::mt19937::default_seed);

// set up dofhandler
dof_handler.reinit(tria);
dof_handler.distribute_dofs(*fe);

// set up constraints
IndexSet locally_relevant_dofs;
DoFTools::extract_locally_relevant_dofs(dof_handler,
locally_relevant_dofs);
constraint.reinit(locally_relevant_dofs);
VectorTools::interpolate_boundary_values(*_mapping,
dof_handler,
0,
Functions::ZeroFunction<dim>(
n_components),
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] =
std::make_shared<MGTwoLevelTransferNonNested<dim, VectorType>>();
transfers[l + 1]->reinit(dof_handlers[l + 1],
dof_handlers[l],
*mappings[l + 1],
*mappings[l],
constraints[l + 1],
constraints[l]);
}

MGTransferGlobalCoarsening<dim, VectorType> transfer(
transfers,
[&](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_handlers[max_level],
operators[max_level],
operators,
transfer);

deallog << dim << ' ' << fe_degree_fine << ' ' << n_refinements << ' '
<< n_components << ' ' << "simplex" << ' '
<< solver_control.last_step() << std::endl;
}

int
main(int argc, char **argv)
{
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
MPILogInitAll all;

deallog.precision(8);
static constexpr unsigned int dim = 2;
static constexpr unsigned int max_simplex_degree = 2;

const double factor = .36;
for (unsigned int n_refinements = 2; n_refinements <= 3; ++n_refinements)
for (unsigned int degree = 1; degree <= max_simplex_degree; ++degree)
test<dim, 1>(n_refinements, degree, factor);

for (unsigned int n_refinements = 2; n_refinements <= 3; ++n_refinements)
for (unsigned int degree = 1; degree <= max_simplex_degree; ++degree)
test<dim, dim>(n_refinements, degree, factor);
}

0 comments on commit c66425a

Please sign in to comment.