--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this function tests the correctness of the implementation of matrix free
+// operations in getting the function values, the function gradients, and the
+// function Laplacians on a cartesian mesh (hyper cube). This tests whether
+// cartesian meshes are treated correctly. The test case is without any
+// constraints
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+template <int dim>
+class CompareFunction : public Function<dim>
+{
+public:
+ CompareFunction()
+ : Function<dim>(dim)
+ {}
+
+ virtual double
+ value(const Point<dim> &p, const unsigned int component) const
+ {
+ double value = (1.2 - 0.1 * component) * p[0] + 0.4;
+ for (unsigned int d = 1; d < dim; ++d)
+ value -= (2.7 - 0.2 * component) * d * p[d];
+ return value;
+ }
+
+ virtual Tensor<1, dim>
+ gradient(const Point<dim> &p, const unsigned int component) const
+ {
+ Tensor<1, dim> grad;
+ grad[0] = 1.2 - 0.1 * component;
+ for (unsigned int d = 1; d < dim; ++d)
+ grad[d] = -(2.7 - 0.2 * component) * d;
+ return grad;
+ }
+};
+
+
+
+template <int dim,
+ int fe_degree,
+ int n_q_points_1d = fe_degree + 1,
+ typename Number = double>
+class MatrixFreeTest
+{
+public:
+ MatrixFreeTest(const MatrixFree<dim, Number> &data_in)
+ : data(data_in){};
+
+ MatrixFreeTest(const MatrixFree<dim, Number> &data_in,
+ const Mapping<dim> &mapping)
+ : data(data_in){};
+
+ virtual ~MatrixFreeTest()
+ {}
+
+ // make function virtual to allow derived classes to define a different
+ // function
+ virtual void
+ cell(const MatrixFree<dim, Number> &data,
+ Vector<Number> &,
+ const Vector<Number> &src,
+ const std::pair<unsigned int, unsigned int> &cell_range) const
+ {
+ FEEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_eval(data);
+
+ CompareFunction<dim> function;
+
+ for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+ {
+ fe_eval.reinit(cell);
+ fe_eval.read_dof_values(src);
+ fe_eval.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+
+ for (unsigned int j = 0; j < data.n_active_entries_per_cell_batch(cell);
+ ++j)
+ for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
+ {
+ ++cell_times;
+ Point<dim> p;
+ for (unsigned int d = 0; d < dim; ++d)
+ p[d] = fe_eval.quadrature_point(q)[d][j];
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ cell_errors[0][d] +=
+ std::abs(fe_eval.get_value(q)[d][j] - function.value(p, d));
+ for (unsigned int e = 0; e < dim; ++e)
+ cell_errors[1][d] +=
+ std::abs(fe_eval.get_gradient(q)[d][e][j] -
+ function.gradient(p, d)[e]);
+ }
+ double divergence = 0;
+ for (unsigned int d = 0; d < dim; ++d)
+ divergence += function.gradient(p, d)[d];
+ cell_errors[2][0] +=
+ std::abs(fe_eval.get_divergence(q)[j] - divergence);
+ }
+ }
+ }
+
+ virtual void
+ face(const MatrixFree<dim, Number> &data,
+ Vector<Number> &,
+ const Vector<Number> &src,
+ const std::pair<unsigned int, unsigned int> &face_range) const
+ {
+ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_evalm(data,
+ true);
+ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_evalp(
+ data, false);
+
+ CompareFunction<dim> function;
+
+ for (unsigned int face = face_range.first; face < face_range.second; ++face)
+ {
+ fe_evalm.reinit(face);
+ fe_evalm.read_dof_values(src);
+ fe_evalm.evaluate(EvaluationFlags::values | EvaluationFlags::gradients |
+ EvaluationFlags::hessians);
+ fe_evalp.reinit(face);
+ fe_evalp.read_dof_values(src);
+ fe_evalp.evaluate(EvaluationFlags::values | EvaluationFlags::gradients |
+ EvaluationFlags::hessians);
+
+ for (unsigned int j = 0; j < VectorizedArray<Number>::size(); ++j)
+ {
+ // skip empty components in VectorizedArray
+ if (data.get_face_info(face).cells_interior[j] ==
+ numbers::invalid_unsigned_int)
+ break;
+ for (unsigned int q = 0; q < fe_evalm.n_q_points; ++q)
+ {
+ ++facem_times;
+ ++facep_times;
+ Point<dim> p;
+
+ // interior face
+ for (unsigned int d = 0; d < dim; ++d)
+ p[d] = fe_evalm.quadrature_point(q)[d][j];
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ facem_errors[0][d] += std::abs(fe_evalm.get_value(q)[d][j] -
+ function.value(p, d));
+ for (unsigned int e = 0; e < dim; ++e)
+ {
+ facem_errors[1][d] +=
+ std::abs(fe_evalm.get_gradient(q)[d][e][j] -
+ function.gradient(p, d)[e]);
+ }
+ }
+ double divergence = 0;
+ for (unsigned int d = 0; d < dim; ++d)
+ divergence += function.gradient(p, d)[d];
+ facem_errors[2][0] +=
+ std::abs(fe_evalm.get_divergence(q)[j] - divergence);
+
+ // exterior face
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ facep_errors[0][d] += std::abs(fe_evalp.get_value(q)[d][j] -
+ function.value(p, d));
+ for (unsigned int e = 0; e < dim; ++e)
+ facep_errors[1][d] +=
+ std::abs(fe_evalp.get_gradient(q)[d][e][j] -
+ function.gradient(p, d)[e]);
+ }
+ facep_errors[2][0] +=
+ std::abs(fe_evalp.get_divergence(q)[j] - divergence);
+ }
+ }
+ }
+ }
+
+ virtual void
+ boundary(const MatrixFree<dim, Number> &data,
+ Vector<Number> &,
+ const Vector<Number> &src,
+ const std::pair<unsigned int, unsigned int> &face_range) const
+ {
+ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_evalm(data,
+ true);
+
+ CompareFunction<dim> function;
+
+ for (unsigned int face = face_range.first; face < face_range.second; ++face)
+ {
+ fe_evalm.reinit(face);
+ fe_evalm.read_dof_values(src);
+ fe_evalm.evaluate(EvaluationFlags::values | EvaluationFlags::gradients |
+ EvaluationFlags::hessians);
+
+ for (unsigned int j = 0; j < VectorizedArray<Number>::size(); ++j)
+ {
+ // skip empty components in VectorizedArray
+ if (data.get_face_info(face).cells_interior[j] ==
+ numbers::invalid_unsigned_int)
+ break;
+ for (unsigned int q = 0; q < fe_evalm.n_q_points; ++q)
+ {
+ ++boundary_times;
+ Point<dim> p;
+ for (unsigned int d = 0; d < dim; ++d)
+ p[d] = fe_evalm.quadrature_point(q)[d][j];
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ boundary_errors[0][d] += std::abs(
+ fe_evalm.get_value(q)[d][j] - function.value(p, d));
+ for (unsigned int e = 0; e < dim; ++e)
+ boundary_errors[1][d] +=
+ std::abs(fe_evalm.get_gradient(q)[d][e][j] -
+ function.gradient(p, d)[e]);
+ }
+ double divergence = 0;
+ for (unsigned int d = 0; d < dim; ++d)
+ divergence += function.gradient(p, d)[d];
+ boundary_errors[2][0] +=
+ std::abs(fe_evalm.get_divergence(q)[j] - divergence);
+ }
+ }
+ }
+ }
+
+
+
+ static void
+ print_error(const std::string &text,
+ const dealii::ndarray<double, 3, dim> &array,
+ const unsigned long long n_times)
+ {
+ deallog << "Error " << std::left << std::setw(6) << text << " values: ";
+ for (unsigned int d = 0; d < dim; ++d)
+ deallog << array[0][d] / n_times << " ";
+ deallog << std::endl;
+ deallog << "Error " << std::left << std::setw(6) << text << " gradients: ";
+ for (unsigned int d = 0; d < dim; ++d)
+ deallog << array[1][d] / n_times << " ";
+ deallog << std::endl;
+ deallog << "Error " << std::left << std::setw(6) << text << " divergence: ";
+ deallog << array[2][0] / n_times << " ";
+ deallog << std::endl;
+ }
+
+ void
+ test_functions(const Vector<Number> &src) const
+ {
+ for (unsigned int d = 0; d < dim; ++d)
+ for (unsigned int i = 0; i < 3; ++i)
+ {
+ cell_errors[i][d] = 0;
+ facem_errors[i][d] = 0;
+ facep_errors[i][d] = 0;
+ boundary_errors[i][d] = 0;
+ }
+ cell_times = facem_times = facep_times = boundary_times = 0;
+
+ Vector<Number> dst_dummy;
+ data.loop(&MatrixFreeTest::cell,
+ &MatrixFreeTest::face,
+ &MatrixFreeTest::boundary,
+ this,
+ dst_dummy,
+ src);
+
+ print_error("cell", cell_errors, cell_times);
+ print_error("face-", facem_errors, facem_times);
+ print_error("face+", facep_errors, facep_times);
+ print_error("face b", boundary_errors, boundary_times);
+ deallog << std::endl;
+ };
+
+protected:
+ const MatrixFree<dim, Number> &data;
+ mutable dealii::ndarray<double, 3, dim> cell_errors, facem_errors,
+ facep_errors, boundary_errors;
+ mutable unsigned long long cell_times, facem_times, facep_times,
+ boundary_times;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void
+do_test(const DoFHandler<dim> &dof,
+ const AffineConstraints<double> &constraints)
+{
+ deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+ // use this for info on problem
+ // std::cout << "Number of cells: " <<
+ // dof.get_triangulation().n_active_cells()
+ // << std::endl;
+ // std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
+ // std::cout << "Number of constraints: " << constraints.n_constraints() <<
+ // std::endl;
+
+ Vector<number> interpolated(dof.n_dofs());
+ VectorTools::interpolate(dof, CompareFunction<dim>(), interpolated);
+
+ constraints.distribute(interpolated);
+ MatrixFree<dim, number> mf_data;
+ {
+ const QGauss<1> quad(dof.get_fe().degree + 1);
+ typename MatrixFree<dim, number>::AdditionalData data;
+ data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+ data.mapping_update_flags = update_gradients | update_quadrature_points;
+ data.mapping_update_flags_boundary_faces =
+ update_gradients | update_quadrature_points;
+ data.mapping_update_flags_inner_faces =
+ update_gradients | update_quadrature_points;
+ mf_data.reinit(MappingQ1<dim>{}, dof, constraints, quad, data);
+ }
+
+ MatrixFreeTest<dim, fe_degree, fe_degree + 1, number> mf(mf_data);
+ mf.test_functions(interpolated);
+}
+
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria, -1, 1);
+ // 1 refinement will make sure that the Piola transform is 1
+ tria.refine_global(1);
+
+ {
+ FE_RaviartThomasNodal<dim> fe(fe_degree - 1);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+
+ AffineConstraints<double> constraints;
+ constraints.close();
+ if (fe_degree > 2)
+ do_test<dim, -1, double>(dof, constraints);
+ else
+ do_test<dim, fe_degree, double>(dof, constraints);
+ }
+}
+
+
+
+int
+main()
+{
+ initlog();
+
+ deallog << std::setprecision(5);
+ {
+ deallog.push("2d");
+ test<2, 1>();
+ test<2, 2>();
+ test<2, 3>();
+ test<2, 4>();
+ deallog.pop();
+ deallog.push("3d");
+ test<3, 1>();
+ test<3, 2>();
+ test<3, 3>();
+ deallog.pop();
+ }
+}
--- /dev/null
+
+DEAL:2d::Testing FE_RaviartThomasNodal<2>(0)
+DEAL:2d::Error cell values: 0.77942 0.31754
+DEAL:2d::Error cell gradients: 2.7000 1.1000
+DEAL:2d::Error cell divergence: 2.2204e-16
+DEAL:2d::Error face- values: 1.0647 0.43377
+DEAL:2d::Error face- gradients: 2.7000 1.1000
+DEAL:2d::Error face- divergence: 2.2204e-16
+DEAL:2d::Error face+ values: 1.0647 0.43377
+DEAL:2d::Error face+ gradients: 2.7000 1.1000
+DEAL:2d::Error face+ divergence: 2.2204e-16
+DEAL:2d::Error face b values: 1.0647 0.43377
+DEAL:2d::Error face b gradients: 2.7000 1.1000
+DEAL:2d::Error face b divergence: 2.2204e-16
+DEAL:2d::
+DEAL:2d::Testing FE_RaviartThomasNodal<2>(1)
+DEAL:2d::Error cell values: 7.8641e-17 1.7193e-16
+DEAL:2d::Error cell gradients: 3.6391e-16 3.5774e-16
+DEAL:2d::Error cell divergence: 3.3923e-16
+DEAL:2d::Error face- values: 3.2382e-17 8.7893e-17
+DEAL:2d::Error face- gradients: 2.2204e-16 4.6259e-16
+DEAL:2d::Error face- divergence: 4.4409e-16
+DEAL:2d::Error face+ values: 3.2382e-17 8.7893e-17
+DEAL:2d::Error face+ gradients: 4.6259e-16 3.1456e-16
+DEAL:2d::Error face+ divergence: 4.6259e-16
+DEAL:2d::Error face b values: 1.0640e-16 7.8641e-17
+DEAL:2d::Error face b gradients: 5.8287e-16 5.3661e-16
+DEAL:2d::Error face b divergence: 6.9389e-16
+DEAL:2d::
+DEAL:2d::Testing FE_RaviartThomasNodal<2>(2)
+DEAL:2d::Error cell values: 2.6498e-16 1.9418e-16
+DEAL:2d::Error cell gradients: 1.9949e-15 1.3982e-15
+DEAL:2d::Error cell divergence: 1.9186e-15
+DEAL:2d::Error face- values: 1.8757e-16 1.8388e-16
+DEAL:2d::Error face- gradients: 3.2058e-15 1.3323e-15
+DEAL:2d::Error face- divergence: 3.3445e-15
+DEAL:2d::Error face+ values: 9.9096e-17 1.2750e-16
+DEAL:2d::Error face+ gradients: 2.2482e-15 2.3870e-15
+DEAL:2d::Error face+ divergence: 3.4417e-15
+DEAL:2d::Error face b values: 3.0358e-16 3.7643e-16
+DEAL:2d::Error face b gradients: 4.0870e-15 2.9213e-15
+DEAL:2d::Error face b divergence: 5.1972e-15
+DEAL:2d::
+DEAL:2d::Testing FE_RaviartThomasNodal<2>(3)
+DEAL:2d::Error cell values: 2.5709e-16 2.3020e-16
+DEAL:2d::Error cell gradients: 3.3573e-15 3.1575e-15
+DEAL:2d::Error cell divergence: 2.8955e-15
+DEAL:2d::Error face- values: 1.7764e-16 2.3315e-16
+DEAL:2d::Error face- gradients: 2.7867e-15 3.4972e-15
+DEAL:2d::Error face- divergence: 2.2982e-15
+DEAL:2d::Error face+ values: 1.4502e-16 1.7902e-16
+DEAL:2d::Error face+ gradients: 3.9191e-15 3.1086e-15
+DEAL:2d::Error face+ divergence: 2.8533e-15
+DEAL:2d::Error face b values: 4.0662e-16 2.8276e-16
+DEAL:2d::Error face b gradients: 6.0452e-15 6.1173e-15
+DEAL:2d::Error face b divergence: 5.0460e-15
+DEAL:2d::
+DEAL:3d::Testing FE_RaviartThomasNodal<3>(0)
+DEAL:3d::Error cell values: 1.5588 1.4434 0.66395
+DEAL:3d::Error cell gradients: 8.1000 6.1000 3.3000
+DEAL:3d::Error cell divergence: 4.4409e-16
+DEAL:3d::Error face- values: 1.9392 1.7956 0.82597
+DEAL:3d::Error face- gradients: 8.1000 6.1000 3.3000
+DEAL:3d::Error face- divergence: 3.7007e-16
+DEAL:3d::Error face+ values: 1.9392 1.7956 0.82597
+DEAL:3d::Error face+ gradients: 8.1000 6.1000 3.3000
+DEAL:3d::Error face+ divergence: 4.4409e-16
+DEAL:3d::Error face b values: 1.9392 1.7956 0.82597
+DEAL:3d::Error face b gradients: 8.1000 6.1000 3.3000
+DEAL:3d::Error face b divergence: 4.0708e-16
+DEAL:3d::
+DEAL:3d::Testing FE_RaviartThomasNodal<3>(1)
+DEAL:3d::Error cell values: 2.9323e-16 3.6905e-16 2.5109e-16
+DEAL:3d::Error cell gradients: 1.0465e-15 1.0043e-15 1.2002e-15
+DEAL:3d::Error cell divergence: 1.1801e-15
+DEAL:3d::Error face- values: 1.8478e-16 2.6355e-16 1.8799e-16
+DEAL:3d::Error face- gradients: 7.1342e-16 7.8744e-16 1.0352e-15
+DEAL:3d::Error face- divergence: 1.0938e-15
+DEAL:3d::Error face+ values: 1.8478e-16 2.6355e-16 1.8799e-16
+DEAL:3d::Error face+ gradients: 9.7864e-16 7.4015e-16 1.1164e-15
+DEAL:3d::Error face+ divergence: 1.3158e-15
+DEAL:3d::Error face b values: 2.9555e-16 3.3975e-16 2.6985e-16
+DEAL:3d::Error face b gradients: 1.1472e-15 1.1030e-15 1.4006e-15
+DEAL:3d::Error face b divergence: 1.6324e-15
+DEAL:3d::
+DEAL:3d::Testing FE_RaviartThomasNodal<3>(2)
+DEAL:3d::Error cell values: 8.0941e-16 8.5099e-16 8.2424e-16
+DEAL:3d::Error cell gradients: 5.5828e-15 5.0810e-15 4.6376e-15
+DEAL:3d::Error cell divergence: 5.7940e-15
+DEAL:3d::Error face- values: 8.3182e-16 7.5558e-16 7.5558e-16
+DEAL:3d::Error face- gradients: 6.7793e-15 5.4089e-15 4.6999e-15
+DEAL:3d::Error face- divergence: 7.8224e-15
+DEAL:3d::Error face+ values: 5.1581e-16 4.8728e-16 6.1290e-16
+DEAL:3d::Error face+ gradients: 5.9964e-15 5.3950e-15 5.5135e-15
+DEAL:3d::Error face+ divergence: 9.3305e-15
+DEAL:3d::Error face b values: 1.1222e-15 1.2358e-15 1.0999e-15
+DEAL:3d::Error face b gradients: 8.5892e-15 7.4477e-15 6.4367e-15
+DEAL:3d::Error face b divergence: 1.1618e-14
+DEAL:3d::
projects / dealii.git / commitdiff