]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Matrix-free for Raviart-Thomas: Test correct interpolation (non-affine)
authorMartin Kronbichler <martin.kronbichler@uni-a.de>
Thu, 7 Sep 2023 18:53:25 +0000 (20:53 +0200)
committerMartin Kronbichler <martin.kronbichler@uni-a.de>
Thu, 7 Sep 2023 18:53:25 +0000 (20:53 +0200)
tests/matrix_free/interpolate_rt_deformed.cc [new file with mode: 0644]
tests/matrix_free/interpolate_rt_deformed.output [new file with mode: 0644]

diff --git a/tests/matrix_free/interpolate_rt_deformed.cc b/tests/matrix_free/interpolate_rt_deformed.cc
new file mode 100644 (file)
index 0000000..8b0374c
--- /dev/null
@@ -0,0 +1,529 @@
+// ---------------------------------------------------------------------
+//
+// 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();
+  }
+}
diff --git a/tests/matrix_free/interpolate_rt_deformed.output b/tests/matrix_free/interpolate_rt_deformed.output
new file mode 100644 (file)
index 0000000..6f47535
--- /dev/null
@@ -0,0 +1,197 @@
+
+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::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.