--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// 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) with deformed
+// elements. This tests whether non-affine geometries with more complicated
+// terms for the Jacobian 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/manifold_lib.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.5 * component) * p[0] * p[0] + 0.4 + component;
+ for (unsigned int d = 1; d < dim; ++d)
+ value -= (2.7 - 0.6 * component) * d * p[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.5 * component) * p[0] * 2;
+ for (unsigned int d = 1; d < dim; ++d)
+ grad[d] = -(2.7 - 0.6 * component) * d * p[d] * 2;
+ 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)
+ {
+ 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)) *
+ fe_eval.JxW(q)[j];
+ 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]) *
+ fe_eval.JxW(q)[j];
+ }
+ 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) *
+ fe_eval.JxW(q)[j];
+ }
+ }
+ }
+
+ 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)
+ {
+ 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)) *
+ fe_evalm.JxW(q)[j];
+ 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]) *
+ fe_evalm.JxW(q)[j];
+ }
+ }
+ 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) *
+ fe_evalm.JxW(q)[j];
+
+ // 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)) *
+ fe_evalm.JxW(q)[j];
+ 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]) *
+ fe_evalm.JxW(q)[j];
+ }
+ facep_errors[2][0] +=
+ std::abs(fe_evalp.get_divergence(q)[j] - divergence) *
+ fe_evalm.JxW(q)[j];
+ }
+ }
+ }
+ }
+
+ 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)
+ {
+ 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)) *
+ fe_evalm.JxW(q)[j];
+ 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]) *
+ fe_evalm.JxW(q)[j];
+ }
+ 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) *
+ fe_evalm.JxW(q)[j];
+ }
+ }
+ }
+ }
+
+
+
+ static void
+ print_error(const std::string &text,
+ const dealii::ndarray<double, 3, dim> &array)
+ {
+ deallog << "Error " << std::left << std::setw(6) << text << " values: ";
+ for (unsigned int d = 0; d < dim; ++d)
+ deallog << array[0][d] << " ";
+ deallog << std::endl;
+ deallog << "Error " << std::left << std::setw(6) << text << " gradients: ";
+ for (unsigned int d = 0; d < dim; ++d)
+ deallog << array[1][d] << " ";
+ deallog << std::endl;
+ deallog << "Error " << std::left << std::setw(6) << text << " divergence: ";
+ deallog << array[2][0] << " ";
+ 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;
+ }
+
+ Vector<Number> dst_dummy;
+ data.loop(&MatrixFreeTest::cell,
+ &MatrixFreeTest::face,
+ &MatrixFreeTest::boundary,
+ this,
+ dst_dummy,
+ src);
+
+ print_error("cell", cell_errors);
+ print_error("face-", facem_errors);
+ print_error("face+", facep_errors);
+ print_error("face b", boundary_errors);
+ deallog << std::endl;
+ };
+
+protected:
+ const MatrixFree<dim, Number> &data;
+ mutable dealii::ndarray<double, 3, dim> cell_errors, facem_errors,
+ facep_errors, boundary_errors;
+};
+
+
+
+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;
+
+ MappingQ<dim> mapping(dof.get_fe().degree);
+ Vector<number> interpolated(dof.n_dofs());
+ VectorTools::interpolate(mapping, 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(mapping, dof, constraints, quad, data);
+ }
+
+ MatrixFreeTest<dim, fe_degree, fe_degree + 1, number> mf(mf_data);
+ mf.test_functions(interpolated);
+}
+
+
+
+template <int dim>
+class DeformedCubeManifold : public dealii::ChartManifold<dim, dim, dim>
+{
+public:
+ DeformedCubeManifold(const double left,
+ const double right,
+ const double deformation,
+ const unsigned int frequency = 1)
+ : left(left)
+ , right(right)
+ , deformation(deformation)
+ , frequency(frequency)
+ {}
+
+ dealii::Point<dim>
+ push_forward(const dealii::Point<dim> &chart_point) const override
+ {
+ double sinval = deformation;
+ for (unsigned int d = 0; d < dim; ++d)
+ sinval *= std::sin(frequency * dealii::numbers::PI *
+ (chart_point(d) - left) / (right - left));
+ dealii::Point<dim> space_point;
+ for (unsigned int d = 0; d < dim; ++d)
+ space_point(d) = chart_point(d) + sinval;
+ return space_point;
+ }
+
+ dealii::Point<dim>
+ pull_back(const dealii::Point<dim> &space_point) const override
+ {
+ dealii::Point<dim> x = space_point;
+ dealii::Point<dim> one;
+ for (unsigned int d = 0; d < dim; ++d)
+ one(d) = 1.;
+
+ // Newton iteration to solve the nonlinear equation given by the point
+ dealii::Tensor<1, dim> sinvals;
+ for (unsigned int d = 0; d < dim; ++d)
+ sinvals[d] = std::sin(frequency * dealii::numbers::PI * (x(d) - left) /
+ (right - left));
+
+ double sinval = deformation;
+ for (unsigned int d = 0; d < dim; ++d)
+ sinval *= sinvals[d];
+ dealii::Tensor<1, dim> residual = space_point - x - sinval * one;
+ unsigned int its = 0;
+ while (residual.norm() > 1e-12 && its < 100)
+ {
+ dealii::Tensor<2, dim> jacobian;
+ for (unsigned int d = 0; d < dim; ++d)
+ jacobian[d][d] = 1.;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ double sinval_der = deformation * frequency / (right - left) *
+ dealii::numbers::PI *
+ std::cos(frequency * dealii::numbers::PI *
+ (x(d) - left) / (right - left));
+ for (unsigned int e = 0; e < dim; ++e)
+ if (e != d)
+ sinval_der *= sinvals[e];
+ for (unsigned int e = 0; e < dim; ++e)
+ jacobian[e][d] += sinval_der;
+ }
+
+ x += invert(jacobian) * residual;
+
+ for (unsigned int d = 0; d < dim; ++d)
+ sinvals[d] = std::sin(frequency * dealii::numbers::PI *
+ (x(d) - left) / (right - left));
+
+ sinval = deformation;
+ for (unsigned int d = 0; d < dim; ++d)
+ sinval *= sinvals[d];
+ residual = space_point - x - sinval * one;
+ ++its;
+ }
+ AssertThrow(residual.norm() < 1e-12,
+ dealii::ExcMessage("Newton for point did not converge."));
+ return x;
+ }
+
+ std::unique_ptr<dealii::Manifold<dim>>
+ clone() const override
+ {
+ return std::make_unique<DeformedCubeManifold<dim>>(left,
+ right,
+ deformation,
+ frequency);
+ }
+
+private:
+ const double left;
+ const double right;
+ const double deformation;
+ const unsigned int frequency;
+};
+
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ Triangulation<dim> tria;
+ Point<dim> left;
+ Point<dim> right;
+ right[0] = 2.;
+ right[1] = 1.;
+ if (dim > 2)
+ right[2] = 1.;
+ GridGenerator::hyper_rectangle(tria, left, right);
+ DeformedCubeManifold<dim> manifold(0, 1.0, 0.05, 1.);
+ tria.set_all_manifold_ids(1);
+ tria.set_manifold(1, manifold);
+ tria.refine_global(1);
+
+ {
+ FE_RaviartThomasNodal<dim> fe(fe_degree - 1);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+
+ AffineConstraints<double> constraints;
+ constraints.close();
+ deallog.push("L1");
+ if (fe_degree > 2)
+ do_test<dim, -1, double>(dof, constraints);
+ else
+ do_test<dim, fe_degree, double>(dof, constraints);
+ deallog.pop();
+
+ tria.refine_global(1);
+ dof.distribute_dofs(fe);
+ deallog.push("L2");
+ if (fe_degree > 2)
+ do_test<dim, -1, double>(dof, constraints);
+ else
+ do_test<dim, fe_degree, double>(dof, constraints);
+ deallog.pop();
+ }
+}
+
+
+
+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:L1::Testing FE_RaviartThomasNodal<2>(0)
+DEAL:2d:L1::Error cell values: 0.84082 0.80829
+DEAL:2d:L1::Error cell gradients: 6.7856 4.0124
+DEAL:2d:L1::Error cell divergence: 1.3856
+DEAL:2d:L1::Error face- values: 1.6772 1.5958
+DEAL:2d:L1::Error face- gradients: 10.686 6.9062
+DEAL:2d:L1::Error face- divergence: 3.3000
+DEAL:2d:L1::Error face+ values: 1.8022 1.4208
+DEAL:2d:L1::Error face+ gradients: 10.686 6.9062
+DEAL:2d:L1::Error face+ divergence: 3.3000
+DEAL:2d:L1::Error face b values: 3.6044 3.0166
+DEAL:2d:L1::Error face b gradients: 21.371 13.812
+DEAL:2d:L1::Error face b divergence: 6.6000
+DEAL:2d:L1::
+DEAL:2d:L2::Testing FE_RaviartThomasNodal<2>(0)
+DEAL:2d:L2::Error cell values: 0.39688 0.40699
+DEAL:2d:L2::Error cell gradients: 6.1965 3.5312
+DEAL:2d:L2::Error cell divergence: 0.72001
+DEAL:2d:L2::Error face- values: 2.5592 2.4054
+DEAL:2d:L2::Error face- gradients: 29.152 17.800
+DEAL:2d:L2::Error face- divergence: 4.9440
+DEAL:2d:L2::Error face+ values: 2.6998 2.2721
+DEAL:2d:L2::Error face+ gradients: 28.810 17.584
+DEAL:2d:L2::Error face+ divergence: 4.9697
+DEAL:2d:L2::Error face b values: 1.7815 1.4453
+DEAL:2d:L2::Error face b gradients: 19.336 11.084
+DEAL:2d:L2::Error face b divergence: 3.4879
+DEAL:2d:L2::
+DEAL:2d:L1::Testing FE_RaviartThomasNodal<2>(1)
+DEAL:2d:L1::Error cell values: 0.21645 0.22161
+DEAL:2d:L1::Error cell gradients: 1.4387 1.0370
+DEAL:2d:L1::Error cell divergence: 0.42416
+DEAL:2d:L1::Error face- values: 0.13278 0.23717
+DEAL:2d:L1::Error face- gradients: 3.8120 1.3753
+DEAL:2d:L1::Error face- divergence: 0.47783
+DEAL:2d:L1::Error face+ values: 0.13874 0.23808
+DEAL:2d:L1::Error face+ gradients: 3.8007 2.4672
+DEAL:2d:L1::Error face+ divergence: 0.70132
+DEAL:2d:L1::Error face b values: 0.34714 0.49328
+DEAL:2d:L1::Error face b gradients: 7.5692 3.9064
+DEAL:2d:L1::Error face b divergence: 2.1365
+DEAL:2d:L1::
+DEAL:2d:L2::Testing FE_RaviartThomasNodal<2>(1)
+DEAL:2d:L2::Error cell values: 0.054207 0.054789
+DEAL:2d:L2::Error cell gradients: 0.63920 0.39444
+DEAL:2d:L2::Error cell divergence: 0.090741
+DEAL:2d:L2::Error face- values: 0.095632 0.16809
+DEAL:2d:L2::Error face- gradients: 5.2642 2.1852
+DEAL:2d:L2::Error face- divergence: 0.39297
+DEAL:2d:L2::Error face+ values: 0.095714 0.16806
+DEAL:2d:L2::Error face+ gradients: 5.1413 2.6374
+DEAL:2d:L2::Error face+ divergence: 0.39710
+DEAL:2d:L2::Error face b values: 0.072933 0.12250
+DEAL:2d:L2::Error face b gradients: 3.5245 1.5231
+DEAL:2d:L2::Error face b divergence: 0.46954
+DEAL:2d:L2::
+DEAL:2d:L1::Testing FE_RaviartThomasNodal<2>(2)
+DEAL:2d:L1::Error cell values: 0.0047119 0.0023984
+DEAL:2d:L1::Error cell gradients: 0.085717 0.029296
+DEAL:2d:L1::Error cell divergence: 0.037673
+DEAL:2d:L1::Error face- values: 0.00065558 0.0035632
+DEAL:2d:L1::Error face- gradients: 0.16207 0.086326
+DEAL:2d:L1::Error face- divergence: 0.072374
+DEAL:2d:L1::Error face+ values: 0.00068398 0.0035694
+DEAL:2d:L1::Error face+ gradients: 0.27606 0.11372
+DEAL:2d:L1::Error face+ divergence: 0.070448
+DEAL:2d:L1::Error face b values: 0.0075049 0.00027532
+DEAL:2d:L1::Error face b gradients: 0.33903 0.075394
+DEAL:2d:L1::Error face b divergence: 0.19221
+DEAL:2d:L1::
+DEAL:2d:L2::Testing FE_RaviartThomasNodal<2>(2)
+DEAL:2d:L2::Error cell values: 0.00078526 0.0035443
+DEAL:2d:L2::Error cell gradients: 0.025118 0.056502
+DEAL:2d:L2::Error cell divergence: 0.010753
+DEAL:2d:L2::Error face- values: 0.0022416 0.013546
+DEAL:2d:L2::Error face- gradients: 0.21582 0.39209
+DEAL:2d:L2::Error face- divergence: 0.063865
+DEAL:2d:L2::Error face+ values: 0.0022387 0.013547
+DEAL:2d:L2::Error face+ gradients: 0.21748 0.36902
+DEAL:2d:L2::Error face+ divergence: 0.064165
+DEAL:2d:L2::Error face b values: 0.00085314 1.8183e-05
+DEAL:2d:L2::Error face b gradients: 0.089290 0.12059
+DEAL:2d:L2::Error face b divergence: 0.052536
+DEAL:2d:L2::
+DEAL:2d:L1::Testing FE_RaviartThomasNodal<2>(3)
+DEAL:2d:L1::Error cell values: 0.0017407 0.0076702
+DEAL:2d:L1::Error cell gradients: 0.039706 0.095855
+DEAL:2d:L1::Error cell divergence: 0.022800
+DEAL:2d:L1::Error face- values: 0.0010959 0.012127
+DEAL:2d:L1::Error face- gradients: 0.11291 0.26940
+DEAL:2d:L1::Error face- divergence: 0.055434
+DEAL:2d:L1::Error face+ values: 0.0010897 0.012127
+DEAL:2d:L1::Error face+ gradients: 0.081715 0.25146
+DEAL:2d:L1::Error face+ divergence: 0.048271
+DEAL:2d:L1::Error face b values: 0.0028746 3.3089e-05
+DEAL:2d:L1::Error face b gradients: 0.23518 0.25632
+DEAL:2d:L1::Error face b divergence: 0.12323
+DEAL:2d:L1::
+DEAL:2d:L2::Testing FE_RaviartThomasNodal<2>(3)
+DEAL:2d:L2::Error cell values: 7.0146e-05 0.00037059
+DEAL:2d:L2::Error cell gradients: 0.0032332 0.0080000
+DEAL:2d:L2::Error cell divergence: 0.0012671
+DEAL:2d:L2::Error face- values: 0.00021264 0.0014135
+DEAL:2d:L2::Error face- gradients: 0.028685 0.071013
+DEAL:2d:L2::Error face- divergence: 0.010209
+DEAL:2d:L2::Error face+ values: 0.00021264 0.0014135
+DEAL:2d:L2::Error face+ gradients: 0.026697 0.056646
+DEAL:2d:L2::Error face+ divergence: 0.0080636
+DEAL:2d:L2::Error face b values: 8.3179e-05 1.0251e-06
+DEAL:2d:L2::Error face b gradients: 0.023934 0.022109
+DEAL:2d:L2::Error face b divergence: 0.0077525
+DEAL:2d:L2::
+DEAL:3d:L1::Testing FE_RaviartThomasNodal<3>(0)
+DEAL:3d:L1::Error cell values: 1.7432 1.3640 0.49183
+DEAL:3d:L1::Error cell gradients: 17.586 12.412 5.5321
+DEAL:3d:L1::Error cell divergence: 2.1651
+DEAL:3d:L1::Error face- values: 5.3910 5.1728 1.9268
+DEAL:3d:L1::Error face- gradients: 44.471 31.919 15.098
+DEAL:3d:L1::Error face- divergence: 6.9450
+DEAL:3d:L1::Error face+ values: 5.3535 3.3645 1.0625
+DEAL:3d:L1::Error face+ gradients: 44.471 31.919 15.098
+DEAL:3d:L1::Error face+ divergence: 6.9450
+DEAL:3d:L1::Error face b values: 12.220 8.9571 3.0912
+DEAL:3d:L1::Error face b gradients: 88.943 63.837 30.196
+DEAL:3d:L1::Error face b divergence: 13.890
+DEAL:3d:L1::
+DEAL:3d:L2::Testing FE_RaviartThomasNodal<3>(0)
+DEAL:3d:L2::Error cell values: 0.84880 0.69655 0.24846
+DEAL:3d:L2::Error cell gradients: 17.166 11.906 4.9397
+DEAL:3d:L2::Error cell divergence: 1.0855
+DEAL:3d:L2::Error face- values: 8.1676 7.2237 2.7157
+DEAL:3d:L2::Error face- gradients: 131.44 91.011 39.859
+DEAL:3d:L2::Error face- divergence: 10.563
+DEAL:3d:L2::Error face+ values: 8.1842 6.1136 2.1595
+DEAL:3d:L2::Error face+ gradients: 129.73 91.230 39.150
+DEAL:3d:L2::Error face+ divergence: 10.562
+DEAL:3d:L2::Error face b values: 6.0673 4.6128 1.5397
+DEAL:3d:L2::Error face b gradients: 86.159 60.569 25.724
+DEAL:3d:L2::Error face b divergence: 7.1112
+DEAL:3d:L2::
+DEAL:3d:L1::Testing FE_RaviartThomasNodal<3>(1)
+DEAL:3d:L1::Error cell values: 0.66818 0.27115 0.26053
+DEAL:3d:L1::Error cell gradients: 4.1680 3.5791 2.5768
+DEAL:3d:L1::Error cell divergence: 0.94875
+DEAL:3d:L1::Error face- values: 1.0170 0.81232 0.76712
+DEAL:3d:L1::Error face- gradients: 15.419 11.678 7.9102
+DEAL:3d:L1::Error face- divergence: 2.6275
+DEAL:3d:L1::Error face+ values: 0.99721 0.82913 0.73940
+DEAL:3d:L1::Error face+ gradients: 14.487 11.520 7.6165
+DEAL:3d:L1::Error face+ divergence: 2.7748
+DEAL:3d:L1::Error face b values: 2.0660 1.3680 0.71188
+DEAL:3d:L1::Error face b gradients: 29.285 22.080 13.242
+DEAL:3d:L1::Error face b divergence: 4.8361
+DEAL:3d:L1::
+DEAL:3d:L2::Testing FE_RaviartThomasNodal<3>(1)
+DEAL:3d:L2::Error cell values: 0.16658 0.067645 0.063432
+DEAL:3d:L2::Error cell gradients: 1.9130 1.5030 0.92443
+DEAL:3d:L2::Error cell divergence: 0.24685
+DEAL:3d:L2::Error face- values: 0.75475 0.54934 0.44999
+DEAL:3d:L2::Error face- gradients: 21.501 15.407 8.7615
+DEAL:3d:L2::Error face- divergence: 1.7590
+DEAL:3d:L2::Error face+ values: 0.75342 0.55080 0.44932
+DEAL:3d:L2::Error face+ gradients: 20.810 14.938 8.6456
+DEAL:3d:L2::Error face+ divergence: 1.7340
+DEAL:3d:L2::Error face b values: 0.51279 0.32049 0.17243
+DEAL:3d:L2::Error face b gradients: 13.912 9.8701 5.0657
+DEAL:3d:L2::Error face b divergence: 1.1911
+DEAL:3d:L2::
+DEAL:3d:L1::Testing FE_RaviartThomasNodal<3>(2)
+DEAL:3d:L1::Error cell values: 0.011206 0.017434 0.020916
+DEAL:3d:L1::Error cell gradients: 0.25129 0.25923 0.28027
+DEAL:3d:L1::Error cell divergence: 0.052910
+DEAL:3d:L1::Error face- values: 0.023217 0.033124 0.052602
+DEAL:3d:L1::Error face- gradients: 0.65036 0.67093 0.89354
+DEAL:3d:L1::Error face- divergence: 0.18134
+DEAL:3d:L1::Error face+ values: 0.023066 0.032185 0.051233
+DEAL:3d:L1::Error face+ gradients: 1.1868 0.83573 0.93673
+DEAL:3d:L1::Error face+ divergence: 0.21746
+DEAL:3d:L1::Error face b values: 0.023435 0.057775 0.060372
+DEAL:3d:L1::Error face b gradients: 1.2755 1.2702 1.2949
+DEAL:3d:L1::Error face b divergence: 0.33420
+DEAL:3d:L1::
+DEAL:3d:L2::Testing FE_RaviartThomasNodal<3>(2)
+DEAL:3d:L2::Error cell values: 0.0016697 0.0055041 0.0073614
+DEAL:3d:L2::Error cell gradients: 0.066790 0.11750 0.14010
+DEAL:3d:L2::Error cell divergence: 0.023259
+DEAL:3d:L2::Error face- values: 0.011054 0.039495 0.052916
+DEAL:3d:L2::Error face- gradients: 0.77293 1.1694 1.3631
+DEAL:3d:L2::Error face- divergence: 0.18844
+DEAL:3d:L2::Error face+ values: 0.011043 0.039429 0.053023
+DEAL:3d:L2::Error face+ gradients: 0.77385 1.1778 1.3552
+DEAL:3d:L2::Error face+ divergence: 0.18832
+DEAL:3d:L2::Error face b values: 0.0030012 0.011224 0.010566
+DEAL:3d:L2::Error face b gradients: 0.37621 0.46520 0.49611
+DEAL:3d:L2::Error face b divergence: 0.11592
+DEAL:3d:L2::