Skip to content

Commit

Permalink
Merge pull request #13603 from peterrum/gc_interpolate_to_mg_fix
Browse files Browse the repository at this point in the history
Fix MGTransferGlobalCoarsening::interpolate_to_mg()
  • Loading branch information
peterrum committed Apr 11, 2022
2 parents 14207e3 + 19cf04c commit ef04483
Show file tree
Hide file tree
Showing 5 changed files with 471 additions and 9 deletions.
7 changes: 5 additions & 2 deletions include/deal.II/matrix_free/constraint_info.h
Original file line number Diff line number Diff line change
Expand Up @@ -125,8 +125,11 @@ namespace internal
this->plain_dof_indices_per_cell.resize(n_cells);
this->constraint_indicator_per_cell.resize(n_cells);

if (use_fast_hanging_node_algorithm &&
dof_handler.get_triangulation().has_hanging_nodes())
// note: has_hanging_nodes() is a global operatrion
const bool has_hanging_nodes =
dof_handler.get_triangulation().has_hanging_nodes();

if (use_fast_hanging_node_algorithm && has_hanging_nodes)
{
hanging_nodes = std::make_unique<HangingNodes<dim>>(
dof_handler.get_triangulation());
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -1481,7 +1481,8 @@ namespace internal

// indices
{
transfer.level_dof_indices_fine.resize(n_dof_indices_fine.back());
transfer.level_dof_indices_fine.resize(n_dof_indices_fine.back(),
numbers::invalid_unsigned_int);

std::vector<types::global_dof_index> local_dof_indices(
transfer.schemes[0].n_dofs_per_cell_coarse);
Expand Down Expand Up @@ -1525,7 +1526,9 @@ namespace internal
dof_handler_coarse,
transfer.schemes[0].n_coarse_cells +
transfer.schemes[1].n_coarse_cells,
use_fast_hanging_node_algorithm(dof_handler_coarse, mg_level_coarse));
constraints_coarse.n_constraints() > 0 &&
use_fast_hanging_node_algorithm(dof_handler_coarse,
mg_level_coarse));

process_cells(
[&](const auto &cell_coarse, const auto &cell_fine) {
Expand Down Expand Up @@ -1574,9 +1577,20 @@ namespace internal
for (unsigned int i = 0;
i < transfer.schemes[1].n_dofs_per_cell_coarse;
++i)
level_dof_indices_fine_1[cell_local_children_indices[c][i]] =
transfer.partitioner_fine->global_to_local(
{
const auto index = transfer.partitioner_fine->global_to_local(
local_dof_indices[lexicographic_numbering_fine[i]]);

Assert(level_dof_indices_fine_1
[cell_local_children_indices[c][i]] ==
numbers::invalid_unsigned_int ||
level_dof_indices_fine_1
[cell_local_children_indices[c][i]] == index,
ExcInternalError());

level_dof_indices_fine_1[cell_local_children_indices[c][i]] =
index;
}
}

// move pointers (only once at the end)
Expand Down Expand Up @@ -1653,7 +1667,7 @@ namespace internal
for (unsigned int j = 0; j < fe->n_dofs_per_cell(); ++j)
transfer.schemes[1]
.restriction_matrix_1d[i * n_child_dofs_1d + j +
c * shift] =
c * shift] +=
matrix(renumbering[i], renumbering[j]);
}
}
Expand Down Expand Up @@ -1694,7 +1708,7 @@ namespace internal
transfer.schemes[1].restriction_matrix
[i * n_dofs_per_cell *
GeometryInfo<dim>::max_children_per_cell +
j + c * n_dofs_per_cell] = matrix(i, j);
j + c * n_dofs_per_cell] += matrix(i, j);
}
}
}
Expand Down Expand Up @@ -2014,7 +2028,9 @@ namespace internal
transfer.constraint_info.reinit(
dof_handler_coarse,
cell_no.back(),
use_fast_hanging_node_algorithm(dof_handler_coarse, mg_level_coarse));
constraints_coarse.n_constraints() > 0 &&
use_fast_hanging_node_algorithm(dof_handler_coarse,
mg_level_coarse));

process_cells([&](const auto &cell_coarse, const auto &cell_fine) {
const auto fe_pair_no =
Expand Down
201 changes: 201 additions & 0 deletions tests/multigrid-global-coarsening/interpolate_01.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,201 @@
// ---------------------------------------------------------------------
//
// Copyright (C) 2020 - 2021 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 MGTransferGlobalCoarsening::interpolate_to_mg() for FE_Q and FE_DGQ.
*/

#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/quadrature_lib.h>

#include <deal.II/distributed/tria.h>

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/mapping_q.h>

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>

#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>

#include <deal.II/multigrid/mg_constrained_dofs.h>
#include <deal.II/multigrid/mg_transfer_global_coarsening.h>

#include <deal.II/numerics/vector_tools.h>

#include "../tests.h"


template <int dim>
void
create_quadrant(Triangulation<dim> &tria, const unsigned int n_refinements)
{
GridGenerator::hyper_cube(tria, -1.0, +1.0);

if (n_refinements == 0)
return;

tria.refine_global(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.0)
flag = false;
if (flag)
cell->set_refine_flag();
}
tria.execute_coarsening_and_refinement();
}

AssertDimension(tria.n_global_levels() - 1, n_refinements);
}



template <int dim>
std::shared_ptr<Utilities::MPI::Partitioner>
create_partitioner(const DoFHandler<dim> &dof_handler)
{
return std::make_shared<Utilities::MPI::Partitioner>(
dof_handler.locally_owned_dofs(),
DoFTools::extract_locally_active_dofs(dof_handler),
dof_handler.get_communicator());
}



template <int dim>
class RightHandSideFunction : public Function<dim>
{
public:
RightHandSideFunction()
: Function<dim>(1)
{}

virtual double
value(const Point<dim> &p, const unsigned int component = 0) const
{
(void)component;

return p[0];
}
};



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

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

parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
create_quadrant(tria, n_refinements);

const auto trias =
MGTransferGlobalCoarseningTools::create_geometric_coarsening_sequence(tria);

MGLevelObject<DoFHandler<dim>> dof_handlers(min_level, max_level);
MGLevelObject<MGTwoLevelTransfer<dim, VectorType>> transfers(min_level,
max_level);

for (auto l = min_level; l <= max_level; ++l)
{
auto &dof_handler = dof_handlers[l];

dof_handler.reinit(*trias[l]);

if (do_dg)
dof_handler.distribute_dofs(FE_DGQ<dim>{fe_degree_fine});
else
dof_handler.distribute_dofs(FE_Q<dim>{fe_degree_fine});
}

for (unsigned int l = min_level; l < max_level; ++l)
transfers[l + 1].reinit(dof_handlers[l + 1], dof_handlers[l]);

MGTransferGlobalCoarsening<dim, VectorType> transfer(
transfers, [&](const auto l, auto &vec) {
vec.reinit(create_partitioner(dof_handlers[l]));
});

VectorType vec;
vec.reinit(create_partitioner(dof_handlers[max_level]));

RightHandSideFunction<dim> fu;

VectorTools::interpolate(dof_handlers[max_level], fu, vec);

MGLevelObject<VectorType> results(min_level, max_level);

transfer.interpolate_to_mg(dof_handlers[max_level], results, vec);

for (unsigned int l = min_level; l <= max_level; ++l)
{
results[l].update_ghost_values();
Vector<float> norm_per_cell(trias[l]->n_active_cells());
VectorTools::integrate_difference(dof_handlers[l],
results[l],
fu,
norm_per_cell,
QGauss<dim>(fe_degree_fine + 2),
VectorTools::L2_norm);
const double error_L2_norm =
VectorTools::compute_global_error(*trias[l],
norm_per_cell,
VectorTools::L2_norm);


deallog << l << " " << results[l].l2_norm() << " " << error_L2_norm
<< std::endl;
}

deallog << std::endl;
}

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

deallog.precision(8);

for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
for (unsigned int degree = 1; degree <= 4; ++degree)
{
test<2>(n_refinements, degree, true);
test<2>(n_refinements, degree, false);
}
}

0 comments on commit ef04483

Please sign in to comment.