--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Solve Poisson problem and Helmholtz problem on a simplex mesh with
+// continuous elements and compare results between matrix-free and matrix-based
+// implementations.
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+
+template <int dim>
+class PoissonOperator
+{
+public:
+ using VectorType = LinearAlgebra::distributed::Vector<double>;
+
+ PoissonOperator(const MatrixFree<dim, double> &matrix_free,
+ const bool do_helmholtz)
+ : matrix_free(matrix_free)
+ , do_helmholtz(do_helmholtz)
+ {}
+
+ void
+ initialize_dof_vector(VectorType &vec)
+ {
+ matrix_free.initialize_dof_vector(vec);
+ }
+
+ void
+ rhs(VectorType &vec) const
+ {
+ const int dummy = 0;
+
+ matrix_free.template cell_loop<VectorType, int>(
+ [&](const auto &, auto &dst, const auto &, const auto cells) {
+ FEEvaluation<dim, -1, 0, 1, double> phi(matrix_free);
+ for (unsigned int cell = cells.first; cell < cells.second; ++cell)
+ {
+ phi.reinit(cell);
+ for (unsigned int q = 0; q < phi.n_q_points; ++q)
+ phi.submit_value(1.0, q);
+
+ phi.integrate_scatter(true, false, dst);
+ }
+ },
+ vec,
+ dummy,
+ true);
+ }
+
+
+ void
+ vmult(VectorType &dst, const VectorType &src) const
+ {
+ matrix_free.template cell_loop<VectorType, VectorType>(
+ [&](const auto &, auto &dst, const auto &src, const auto cells) {
+ FEEvaluation<dim, -1, 0, 1, double> phi(matrix_free);
+ for (unsigned int cell = cells.first; cell < cells.second; ++cell)
+ {
+ phi.reinit(cell);
+ phi.gather_evaluate(src, do_helmholtz, true);
+
+ for (unsigned int q = 0; q < phi.n_q_points; ++q)
+ {
+ if (do_helmholtz)
+ phi.submit_value(phi.get_value(q), q);
+
+ phi.submit_gradient(phi.get_gradient(q), q);
+ }
+
+ phi.integrate_scatter(do_helmholtz, true, dst);
+ }
+ },
+ dst,
+ src,
+ true);
+ }
+
+private:
+ const MatrixFree<dim, double> &matrix_free;
+ const bool do_helmholtz;
+};
+
+template <int dim>
+void
+test(const unsigned int degree, const bool do_helmholtz)
+{
+ Triangulation<dim> tria;
+
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, dim == 2 ? 16 : 8);
+
+ Simplex::FE_P<dim> fe(degree);
+ Simplex::QGauss<dim> quad(degree + 1);
+ MappingFE<dim> mapping(Simplex::FE_P<dim>(1));
+
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ AffineConstraints<double> constraints;
+ VectorTools::interpolate_boundary_values(
+ mapping, dof_handler, 0, Functions::ZeroFunction<dim>(), constraints);
+ constraints.close();
+
+ const auto solve_and_postprocess =
+ [&](const auto &poisson_operator,
+ auto & x,
+ auto & b) -> std::pair<unsigned int, double> {
+ ReductionControl reduction_control;
+ SolverCG<typename std::remove_reference<decltype(x)>::type> solver(
+ reduction_control);
+ solver.solve(poisson_operator, x, b, PreconditionIdentity());
+
+ if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+ printf("Solved in %d iterations.\n", reduction_control.last_step());
+
+ constraints.distribute(x);
+
+#if 0
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler(dof_handler);
+ x.update_ghost_values();
+ data_out.add_data_vector(dof_handler, x, "solution");
+ data_out.build_patches(mapping, 2);
+ data_out.write_vtu_with_pvtu_record("./", "result", 0, MPI_COMM_WORLD);
+#endif
+
+ return {reduction_control.last_step(), reduction_control.last_value()};
+ };
+
+ const auto mf_algo = [&]() {
+ typename MatrixFree<dim, double>::AdditionalData additional_data;
+ additional_data.mapping_update_flags = update_gradients | update_values;
+
+ MatrixFree<dim, double> matrix_free;
+ matrix_free.reinit(
+ mapping, dof_handler, constraints, quad, additional_data);
+
+ PoissonOperator<dim> poisson_operator(matrix_free, do_helmholtz);
+
+ LinearAlgebra::distributed::Vector<double> x, b;
+ poisson_operator.initialize_dof_vector(x);
+ poisson_operator.initialize_dof_vector(b);
+
+ poisson_operator.rhs(b);
+
+ return solve_and_postprocess(poisson_operator, x, b);
+ };
+
+ const auto mb_algo = [&]() {
+ Vector<double> x, b;
+
+ x.reinit(dof_handler.n_dofs());
+ b.reinit(dof_handler.n_dofs());
+
+ SparseMatrix<double> A;
+
+ DynamicSparsityPattern dsp(dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
+
+ SparsityPattern sparsity_pattern;
+ sparsity_pattern.copy_from(dsp);
+ A.reinit(sparsity_pattern);
+
+ const auto flags = update_values | update_gradients | update_JxW_values;
+
+ FEValues<dim> fe_values(mapping, fe, quad, flags);
+
+ FullMatrix<double> cell_matrix;
+ Vector<double> cell_rhs;
+ std::vector<types::global_dof_index> local_dof_indices;
+
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ {
+ if (cell->is_locally_owned() == false)
+ continue;
+
+ fe_values.reinit(cell);
+
+ const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+ cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+ cell_rhs.reinit(dofs_per_cell);
+
+ for (const auto q : fe_values.quadrature_point_indices())
+ {
+ for (const auto i : fe_values.dof_indices())
+ for (const auto j : fe_values.dof_indices())
+ cell_matrix(i, j) += (fe_values.shape_grad(i, q) * //
+ fe_values.shape_grad(j, q) + //
+ static_cast<double>(do_helmholtz) * //
+ fe_values.shape_value(i, q) * //
+ fe_values.shape_value(j, q)) * //
+ fe_values.JxW(q); //
+
+ for (const unsigned int i : fe_values.dof_indices())
+ cell_rhs(i) += (fe_values.shape_value(i, q) * //
+ 1. * //
+ fe_values.JxW(q)); //
+ }
+
+ local_dof_indices.resize(cell->get_fe().dofs_per_cell);
+ cell->get_dof_indices(local_dof_indices);
+
+ constraints.distribute_local_to_global(
+ cell_matrix, cell_rhs, local_dof_indices, A, b);
+ }
+
+ return solve_and_postprocess(A, x, b);
+ };
+
+ const auto compare = [&](const auto result_mf, const auto result_mb) {
+ AssertDimension(result_mf.first, result_mb.first);
+ Assert(std::abs(result_mf.second - result_mb.second) < 1e-8,
+ ExcNotImplemented());
+
+ deallog << "dim=" << dim << " ";
+ deallog << "degree=" << degree << " ";
+ deallog << "Type=";
+
+ if (do_helmholtz)
+ deallog << "Helmholtz"
+ << " : ";
+ else
+ deallog << "Possion "
+ << " : ";
+
+ deallog << "Convergence step " << result_mf.first << " value "
+ << result_mf.second << "." << std::endl;
+ };
+
+ compare(mf_algo(), mb_algo());
+}
+
+
+int
+main(int argc, char **argv)
+{
+ initlog();
+
+ deallog.depth_file(1);
+
+ Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+
+ test<2>(/*degree=*/1, /*do_helmholtz*/ false);
+ test<2>(/*degree=*/1, /*do_helmholtz*/ true);
+ test<2>(/*degree=*/2, /*do_helmholtz*/ false);
+ test<2>(/*degree=*/2, /*do_helmholtz*/ true);
+
+ test<3>(/*degree=*/1, /*do_helmholtz*/ false);
+ test<3>(/*degree=*/1, /*do_helmholtz*/ true);
+ test<3>(/*degree=*/2, /*do_helmholtz*/ false);
+ test<3>(/*degree=*/2, /*do_helmholtz*/ true);
+}
projects / dealii.git / commitdiff