]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Introduce inverse quadratic approximation for MappingQ::real_to_unit_cell
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 22 Oct 2020 14:37:50 +0000 (16:37 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 28 Oct 2020 13:10:21 +0000 (14:10 +0100)
include/deal.II/fe/mapping_q_internal.h
source/fe/mapping_q_generic.cc
tests/mappings/mapping_q_inverse_quadratic_approximation.cc [new file with mode: 0644]
tests/mappings/mapping_q_inverse_quadratic_approximation.output [new file with mode: 0644]

index e0971d82588ee8e77e808b5c545025c5badd7b01..88f2efbb527b97ef4b71f8aa0d9805475240ff11 100644 (file)
@@ -31,6 +31,8 @@
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping_q_generic.h>
 
+#include <deal.II/grid/grid_tools.h>
+
 #include <deal.II/matrix_free/evaluation_flags.h>
 #include <deal.II/matrix_free/evaluation_template_factory.h>
 #include <deal.II/matrix_free/shape_info.h>
@@ -811,6 +813,202 @@ namespace internal
       return p_unit;
     }
 
+
+
+    /**
+     * A class to compute a quadratic approximation to the inverse map from
+     * real to unit points by a least-squares fit along the mapping support
+     * points. The least squares fit is special in the sense that the
+     * approximation is constructed for the inverse function of a
+     * MappingQGeneric, which is generally a rational function. This allows
+     * for a very cheap evaluation of the inverse map by a simple polynomial
+     * interpolation, which can be used as a better initial guess for
+     * transforming points from real to unit coordinates than an affine
+     * approximation.
+     *
+     * Far away outside the unit cell, this approximation can become
+     * inaccurate for non-affine cell shapes. This must be expected from a
+     * fit of a polynomial to a rational function, and due to the fact that
+     * the region of the least squares fit, the unit cell, is left. Hence,
+     * use this function with care in those situations.
+     */
+    template <int dim, int spacedim>
+    class InverseQuadraticApproximation
+    {
+    public:
+      /**
+       * Number of basis functions in the quadratic approximation.
+       */
+      static constexpr unsigned int n_functions =
+        (spacedim == 1 ? 3 : (spacedim == 2 ? 6 : 10));
+
+      /**
+       * Constructor.
+       */
+      InverseQuadraticApproximation(
+        const std::vector<Point<spacedim>> &mapping_support_points,
+        const std::vector<Point<1>> &       line_support_points,
+        const std::vector<unsigned int> &   renumber)
+        : p_shift(mapping_support_points[0])
+        , scale(1. /
+                mapping_support_points[0].distance(mapping_support_points[1]))
+      {
+        AssertDimension(mapping_support_points.size(), renumber.size());
+        AssertDimension(mapping_support_points.size(),
+                        Utilities::pow(line_support_points.size(), dim));
+
+        const unsigned int n1 = line_support_points.size();
+
+        // For the bi-/trilinear approximation, we cannot build a quadratic
+        // polynomial due to a lack of points (interpolation matrix would get
+        // singular), so pick the affine approximation. Similarly, it is not
+        // entirely clear how to gather enough information for the case dim <
+        // spacedim
+        if (n1 == 2 || dim < spacedim)
+          {
+            const auto affine = GridTools::affine_cell_approximation<dim>(
+              make_array_view(mapping_support_points));
+            DerivativeForm<1, spacedim, dim> A_inv =
+              affine.first.covariant_form().transpose();
+            coefficients[0] = apply_transformation(A_inv, affine.second);
+            for (unsigned int d = 0; d < spacedim; ++d)
+              for (unsigned int e = 0; e < dim; ++e)
+                coefficients[1 + d][e] = A_inv[e][d];
+            is_affine = true;
+            return;
+          }
+
+        SymmetricTensor<2, n_functions> matrix;
+        std::array<double, n_functions> shape_values;
+        for (unsigned int q2 = 0, q = 0; q2 < (dim > 2 ? n1 : 1); ++q2)
+          for (unsigned int q1 = 0; q1 < (dim > 1 ? n1 : 1); ++q1)
+            for (unsigned int q0 = 0; q0 < n1; ++q0, ++q)
+              {
+                // Evaluate quadratic shape functions in point, shifted to the
+                // first support point and scaled by the length between the
+                // first two support points to avoid roundoff issues with
+                // scaling far away from 1.
+                shape_values[0] = 1.;
+                const Tensor<1, spacedim> p_scaled =
+                  (mapping_support_points[renumber[q]] - p_shift) * scale;
+                for (unsigned int d = 0; d < spacedim; ++d)
+                  shape_values[1 + d] = p_scaled[d];
+                for (unsigned int d = 0, c = 0; d < spacedim; ++d)
+                  for (unsigned int e = 0; e <= d; ++e, ++c)
+                    shape_values[1 + spacedim + c] = p_scaled[d] * p_scaled[e];
+
+                // Build lower diagonal of least squares matrix and rhs, the
+                // essential part being that we construct the matrix with the
+                // real points and the right hand side by comparing to the
+                // reference point positions which sets up an inverse
+                // interpolation.
+                for (unsigned int i = 0; i < n_functions; ++i)
+                  for (unsigned int j = 0; j <= i; ++j)
+                    matrix[i][j] += shape_values[i] * shape_values[j];
+                Point<dim> reference_point;
+                reference_point[0] = line_support_points[q0][0];
+                if (dim > 1)
+                  reference_point[1] = line_support_points[q1][0];
+                if (dim > 2)
+                  reference_point[2] = line_support_points[q2][0];
+                for (unsigned int i = 0; i < n_functions; ++i)
+                  coefficients[i] += shape_values[i] * reference_point;
+              }
+
+        // Factorize the matrix A = L * L^T in-place with the
+        // Cholesky-Banachiewicz algorithm. The implementation is similar to
+        // FullMatrix::cholesky() but re-implemented to avoid memory
+        // allocations and some unnecessary divisions which we can do here as
+        // we only need to solve with dim right hand sides.
+        for (unsigned int i = 0; i < n_functions; ++i)
+          {
+            double Lij_sum = 0;
+            for (unsigned int j = 0; j < i; ++j)
+              {
+                double Lik_Ljk_sum = 0;
+                for (unsigned int k = 0; k < j; ++k)
+                  Lik_Ljk_sum += matrix[i][k] * matrix[j][k];
+                matrix[i][j] = matrix[j][j] * (matrix[i][j] - Lik_Ljk_sum);
+                Lij_sum += matrix[i][j] * matrix[i][j];
+              }
+            AssertThrow(matrix[i][i] - Lij_sum >= 0,
+                        ExcMessage("Matrix not positive definite"));
+
+            // Store the inverse in the diagonal since that is the quantity
+            // needed later in the factorization as well as the forward and
+            // backward substitution, minimizing the number of divisions.
+            matrix[i][i] = 1. / std::sqrt(matrix[i][i] - Lij_sum);
+          }
+
+        // Solve lower triangular part, L * y = rhs.
+        for (unsigned int i = 0; i < n_functions; ++i)
+          {
+            Point<dim> sum = coefficients[i];
+            for (unsigned int j = 0; j < i; ++j)
+              sum -= matrix[i][j] * coefficients[j];
+            coefficients[i] = sum * matrix[i][i];
+          }
+
+        // Solve upper triangular part, L^T * x = y (i.e., x = A^{-1} * rhs)
+        for (unsigned int i = n_functions; i > 0;)
+          {
+            --i;
+            Point<dim> sum = coefficients[i];
+            for (unsigned int j = i + 1; j < n_functions; ++j)
+              sum -= matrix[j][i] * coefficients[j];
+            coefficients[i] = sum * matrix[i][i];
+          }
+
+        // Check whether the approximation is indeed affine, allowing to
+        // skip the quadratic terms.
+        is_affine = true;
+        for (unsigned int i = dim + 1; i < n_functions; ++i)
+          if (coefficients[i].norm_square() > 1e-20)
+            {
+              is_affine = false;
+              break;
+            }
+      }
+
+      /**
+       * Evaluate the quadratic approximation.
+       */
+      template <typename Number>
+      Point<dim, Number>
+      compute(const Point<spacedim, Number> &p)
+      {
+        Point<dim, Number> result;
+        for (unsigned int d = 0; d < dim; ++d)
+          result[d] = coefficients[0][d];
+
+        // Shift point to avoid roundoff problems. Spell out the loop
+        // because Number might be a vectorized array.
+        Point<spacedim, Number> p_scaled;
+        for (unsigned int d = 0; d < spacedim; ++d)
+          p_scaled[d] = (p[d] - p_shift[d]) * scale;
+
+        for (unsigned int d = 0; d < spacedim; ++d)
+          result += coefficients[1 + d] * p_scaled[d];
+
+        if (!is_affine)
+          {
+            for (unsigned int d = 0, c = 0; d < spacedim; ++d)
+              for (unsigned int e = 0; e <= d; ++e, ++c)
+                result +=
+                  coefficients[1 + spacedim + c] * (p_scaled[d] * p_scaled[e]);
+          }
+        return result;
+      }
+
+    private:
+      const Point<spacedim>               p_shift;
+      const double                        scale;
+      std::array<Point<dim>, n_functions> coefficients;
+      bool                                is_affine;
+    };
+
+
+
     /**
      * In case the quadrature formula is a tensor product, this is a
      * replacement for maybe_compute_q_points(), maybe_update_Jacobians() and
index 3cd4cbf975dc09edd0eeca1a1c971a37c7a6f31d..4a2d3d3292cb0cdd4a994583d7081ea33e7e7d69 100644 (file)
@@ -723,18 +723,27 @@ MappingQGeneric<dim, spacedim>::transform_points_real_to_unit_cell(
   const ArrayView<const Point<spacedim>> &                    real_points,
   const ArrayView<Point<dim>> &                               unit_points) const
 {
+  // Go to base class functions for dim < spacedim because it is not yet
+  // implemented with optimized code.
+  if (dim < spacedim)
+    {
+      Mapping<dim, spacedim>::transform_points_real_to_unit_cell(cell,
+                                                                 real_points,
+                                                                 unit_points);
+      return;
+    }
+
   AssertDimension(real_points.size(), unit_points.size());
   const std::vector<Point<spacedim>> support_points =
     this->compute_mapping_support_points(cell);
 
   // From the given (high-order) support points, now only pick the first
   // 2^dim points and construct an affine approximation from those.
-  const std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
-    affine_factors = GridTools::affine_cell_approximation<dim>(
-      ArrayView<const Point<spacedim>>(support_points.data(),
-                                       GeometryInfo<dim>::vertices_per_cell));
-  const DerivativeForm<1, spacedim, dim> A_inv =
-    affine_factors.first.covariant_form().transpose();
+  internal::MappingQGenericImplementation::
+    InverseQuadraticApproximation<dim, spacedim>
+      inverse_approximation(support_points,
+                            line_support_points,
+                            renumber_lexicographic_to_hierarchic);
 
   const unsigned int n_points = real_points.size();
   const unsigned int n_lanes  = VectorizedArray<double>::size();
@@ -753,25 +762,11 @@ MappingQGeneric<dim, spacedim>::transform_points_real_to_unit_cell(
             for (unsigned int d = 0; d < spacedim; ++d)
               p_vec[d][j] = real_points[i][d];
 
-        // Compute an initial guess by inverting the affine approximation
-        // A * x_unit + b = x_real
-        Tensor<1, spacedim, VectorizedArray<double>> rhs;
-        for (unsigned int d = 0; d < spacedim; ++d)
-          rhs[d] = affine_factors.second[d];
-        rhs = p_vec - rhs;
-
-        Point<dim, VectorizedArray<double>> initial_guess;
-        for (unsigned int d = 0; d < dim; ++d)
-          {
-            initial_guess[d] = A_inv[d][0] * rhs[0];
-            for (unsigned int e = 1; e < spacedim; ++e)
-              initial_guess[d] += A_inv[d][e] * rhs[e];
-          }
         Point<dim, VectorizedArray<double>> unit_point =
           internal::MappingQGenericImplementation::
             do_transform_real_to_unit_cell_internal<dim, spacedim>(
               p_vec,
-              initial_guess,
+              inverse_approximation.compute(p_vec),
               support_points,
               polynomials_1d,
               renumber_lexicographic_to_hierarchic);
@@ -786,8 +781,7 @@ MappingQGeneric<dim, spacedim>::transform_points_real_to_unit_cell(
             unit_points[i + j] = internal::MappingQGenericImplementation::
               do_transform_real_to_unit_cell_internal<dim, spacedim>(
                 real_points[i + j],
-                Point<dim>(apply_transformation(
-                  A_inv, real_points[i + j] - affine_factors.second)),
+                inverse_approximation.compute(real_points[i + j]),
                 support_points,
                 polynomials_1d,
                 renumber_lexicographic_to_hierarchic);
@@ -799,8 +793,7 @@ MappingQGeneric<dim, spacedim>::transform_points_real_to_unit_cell(
       unit_points[i] = internal::MappingQGenericImplementation::
         do_transform_real_to_unit_cell_internal<dim, spacedim>(
           real_points[i],
-          Point<dim>(apply_transformation(
-            A_inv, real_points[i] - affine_factors.second)),
+          inverse_approximation.compute(real_points[i]),
           support_points,
           polynomials_1d,
           renumber_lexicographic_to_hierarchic);
diff --git a/tests/mappings/mapping_q_inverse_quadratic_approximation.cc b/tests/mappings/mapping_q_inverse_quadratic_approximation.cc
new file mode 100644 (file)
index 0000000..8972746
--- /dev/null
@@ -0,0 +1,141 @@
+// ---------------------------------------------------------------------
+//
+// 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Check InverseQuadraticApproximation used for the initial guess in
+// MappingQGeneric::transform_points_real_to_unit_cell
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.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_fe_field.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/mapping_q_generic.h>
+#include <deal.II/fe/mapping_q_internal.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+template <int dim, int spacedim>
+void
+print_result(const unsigned int                  mapping_degree,
+             const Triangulation<dim, spacedim> &tria,
+             const Point<dim>                    p)
+{
+  deallog << "Testing " << dim << "D with point " << p << std::endl;
+
+  FE_Q<dim>            dummy(mapping_degree);
+  MappingQGeneric<dim> mapping(mapping_degree);
+
+  FEValues<dim> fe_values(mapping,
+                          dummy,
+                          Quadrature<dim>(dummy.get_unit_support_points()),
+                          update_quadrature_points);
+
+  std::vector<Point<1>> mapping_points(
+    QGaussLobatto<1>(mapping_degree + 1).get_points());
+  std::vector<unsigned int> renumber =
+    FETools::lexicographic_to_hierarchic_numbering<dim>(mapping_degree);
+
+  for (const auto &cell : tria.active_cell_iterators())
+    {
+      fe_values.reinit(cell);
+      try
+        {
+          const Point<dim> p_unit =
+            mapping.transform_real_to_unit_cell(cell, p);
+          deallog << "Testing on cell " << cell->id() << " with center "
+                  << cell->center(true) << std::endl;
+          deallog << "Exact inverse:                   " << p_unit << std::endl;
+          deallog << "Affine approximation:            "
+                  << cell->real_to_unit_cell_affine_approximation(p)
+                  << std::endl;
+          internal::MappingQGenericImplementation::
+            InverseQuadraticApproximation<dim, spacedim>
+              approx(fe_values.get_quadrature_points(),
+                     mapping_points,
+                     renumber);
+          deallog << "Inverse quadratic approximation: " << approx.compute(p)
+                  << std::endl
+                  << std::endl;
+        }
+      catch (typename Mapping<dim, spacedim>::ExcTransformationFailed &)
+        {}
+    }
+  deallog << std::endl;
+}
+
+
+
+template <int dim, int spacedim>
+void
+test(const unsigned mapping_degree, const unsigned int n_ref)
+{
+  Triangulation<dim, spacedim> tria;
+  GridGenerator::hyper_shell(tria, Point<dim>(), 0.8, 1., dim == 2 ? 3 : 6);
+  tria.refine_global(n_ref);
+  {
+    const double phi   = 0.56 * numbers::PI;
+    const double theta = 0.49 * numbers::PI;
+    for (unsigned int i = 0; i <= 4; ++i)
+      {
+        const double r = 0.8 + 0.2 * i / 4;
+        Point<dim>   p;
+        p[0] = std::cos(phi) * std::sin(theta) * r;
+        p[1] = std::sin(phi) * std::sin(theta) * r;
+        if (dim == 3)
+          p[2] = std::cos(theta) * r;
+        print_result(mapping_degree, tria, p);
+      }
+  }
+  if (dim == 3)
+    {
+      const double phi   = 0.3 * numbers::PI;
+      const double theta = 0.1;
+      for (unsigned int i = 0; i <= 6; ++i)
+        {
+          const double r = 0.8 + 0.2 * i / 6;
+          Point<dim>   p;
+          p[0] = std::cos(phi) * std::sin(theta) * r;
+          p[1] = std::sin(phi) * std::sin(theta) * r;
+          if (dim == 3)
+            p[2] = std::cos(theta) * r;
+          print_result(mapping_degree, tria, p);
+        }
+    }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(10);
+
+  test<2, 2>(6, 1);
+  test<3, 3>(4, 0);
+}
diff --git a/tests/mappings/mapping_q_inverse_quadratic_approximation.output b/tests/mappings/mapping_q_inverse_quadratic_approximation.output
new file mode 100644 (file)
index 0000000..d8686eb
--- /dev/null
@@ -0,0 +1,355 @@
+
+DEAL::Testing 2D with point -0.1498310826 0.7854420410
+DEAL::Testing on cell 0_1:0 with center 0.8227241336 0.4750000000
+DEAL::Exact inverse:                   1.679734433 2.008122269
+DEAL::Affine approximation:            1.294871897 6.963559088
+DEAL::Inverse quadratic approximation: 1.644896745 2.003800318
+DEAL::
+DEAL::Testing on cell 0_1:1 with center 5.817072296e-17 0.9500000000
+DEAL::Exact inverse:                   0.6800000317 2.003947473
+DEAL::Affine approximation:            0.6577169290 0.9304965238
+DEAL::Inverse quadratic approximation: 0.6639670426 1.893521583
+DEAL::
+DEAL::Testing on cell 0_1:2 with center 0.7361215932 0.4250000000
+DEAL::Exact inverse:                   1.679734435 1.008122258
+DEAL::Affine approximation:            1.388386238 5.963559088
+DEAL::Inverse quadratic approximation: 1.780446266 1.113104041
+DEAL::
+DEAL::Testing on cell 0_1:3 with center 5.204748896e-17 0.8500000000
+DEAL::Exact inverse:                   0.6800000317 1.003947473
+DEAL::Affine approximation:            0.6762718619 -0.06950347624
+DEAL::Inverse quadratic approximation: 0.6705854753 0.9955246392
+DEAL::
+DEAL::Testing on cell 1_1:0 with center -0.8227241336 0.4750000000
+DEAL::Exact inverse:                   -0.3199797939 2.004105574
+DEAL::Affine approximation:            -0.1371549678 3.966937436
+DEAL::Inverse quadratic approximation: -0.2909091648 1.929190316
+DEAL::
+DEAL::Testing on cell 1_1:1 with center -0.8227241336 -0.4750000000
+DEAL::Exact inverse:                   -1.321431976 2.123538334
+DEAL::Affine approximation:            -0.2948718968 13.03644091
+DEAL::Inverse quadratic approximation: -0.9655241914 2.275064795
+DEAL::
+DEAL::Testing on cell 1_1:2 with center -0.7361215932 0.4250000000
+DEAL::Exact inverse:                   -0.3199797942 1.004105575
+DEAL::Affine approximation:            -0.2121143758 2.966937436
+DEAL::Inverse quadratic approximation: -0.3588753056 1.026324172
+DEAL::
+DEAL::Testing on cell 1_1:3 with center -0.7361215932 -0.4250000000
+DEAL::Exact inverse:                   -1.321431983 1.123539113
+DEAL::Affine approximation:            -0.3883862377 12.03644091
+DEAL::Inverse quadratic approximation: -1.203941726 1.460747557
+DEAL::
+DEAL::
+DEAL::Testing 2D with point -0.1591955252 0.8345321686
+DEAL::Testing on cell 0_1:0 with center 0.8227241336 0.4750000000
+DEAL::Exact inverse:                   1.679734435 1.508629900
+DEAL::Affine approximation:            1.344551390 6.773781531
+DEAL::Inverse quadratic approximation: 1.705806959 1.564309847
+DEAL::
+DEAL::Testing on cell 0_1:1 with center 5.817072296e-17 0.9500000000
+DEAL::Exact inverse:                   0.6800000317 1.504194190
+DEAL::Affine approximation:            0.6675742371 0.3636525565
+DEAL::Inverse quadratic approximation: 0.6679057071 1.457711038
+DEAL::
+DEAL::Testing on cell 0_1:2 with center 0.7361215932 0.4250000000
+DEAL::Exact inverse:                   1.679734435 0.5086298996
+DEAL::Affine approximation:            1.443910378 5.773781531
+DEAL::Inverse quadratic approximation: 1.846412785 0.6214152525
+DEAL::
+DEAL::Testing on cell 0_1:3 with center 5.204748896e-17 0.8500000000
+DEAL::Exact inverse:                   0.6800000317 0.5041941900
+DEAL::Affine approximation:            0.6872888532 -0.6363474435
+DEAL::Inverse quadratic approximation: 0.6729135650 0.5116134621
+DEAL::
+DEAL::Testing on cell 1_1:0 with center -0.8227241336 0.4750000000
+DEAL::Exact inverse:                   -0.3199797942 1.504362173
+DEAL::Affine approximation:            -0.1769771533 3.589871026
+DEAL::Inverse quadratic approximation: -0.3233858788 1.488970865
+DEAL::
+DEAL::Testing on cell 1_1:1 with center -0.8227241336 -0.4750000000
+DEAL::Exact inverse:                   -1.321431983 1.631260238
+DEAL::Affine approximation:            -0.3445513904 13.22621847
+DEAL::Inverse quadratic approximation: -1.067765287 1.852528354
+DEAL::
+DEAL::Testing on cell 1_1:2 with center -0.7361215932 0.4250000000
+DEAL::Exact inverse:                   -0.3199797942 0.5043621731
+DEAL::Affine approximation:            -0.2566215243 2.589871026
+DEAL::Inverse quadratic approximation: -0.3901601371 0.5348403961
+DEAL::
+DEAL::Testing on cell 1_1:3 with center -0.7361215932 -0.4250000000
+DEAL::Exact inverse:                   -1.321431983 0.6312603084
+DEAL::Affine approximation:            -0.4439103775 12.22621847
+DEAL::Inverse quadratic approximation: -1.324499457 0.9907864881
+DEAL::
+DEAL::
+DEAL::Testing 2D with point -0.1685599679 0.8836222961
+DEAL::Testing on cell 0_1:0 with center 0.8227241336 0.4750000000
+DEAL::Exact inverse:                   1.679734435 1.009137541
+DEAL::Affine approximation:            1.394230884 6.584003974
+DEAL::Inverse quadratic approximation: 1.765464721 1.098697349
+DEAL::
+DEAL::Testing on cell 0_1:1 with center 5.817072296e-17 0.9500000000
+DEAL::Exact inverse:                   0.6800000317 1.004440907
+DEAL::Affine approximation:            0.6774315452 -0.2031914108
+DEAL::Inverse quadratic approximation: 0.6711021038 0.9970222699
+DEAL::
+DEAL::Testing on cell 0_1:2 with center 0.7361215932 0.4250000000
+DEAL::Exact inverse:                   1.679734435 0.009137540713
+DEAL::Affine approximation:            1.499434517 5.584003974
+DEAL::Inverse quadratic approximation: 1.910725025 0.1005855588
+DEAL::
+DEAL::Testing on cell 0_1:3 with center 5.204748896e-17 0.8500000000
+DEAL::Exact inverse:                   0.6800000317 0.004440907068
+DEAL::Affine approximation:            0.6983058446 -1.203191411
+DEAL::Inverse quadratic approximation: 0.6742612427 0.0003409477391
+DEAL::
+DEAL::Testing on cell 1_1:0 with center -0.8227241336 0.4750000000
+DEAL::Exact inverse:                   -0.3199797942 1.004618772
+DEAL::Affine approximation:            -0.2167993388 3.212804616
+DEAL::Inverse quadratic approximation: -0.3538678741 1.023092227
+DEAL::
+DEAL::Testing on cell 1_1:1 with center -0.8227241336 -0.4750000000
+DEAL::Exact inverse:                   -1.321431983 1.138981502
+DEAL::Affine approximation:            -0.3942308840 13.41599603
+DEAL::Inverse quadratic approximation: -1.171258833 1.403869886
+DEAL::
+DEAL::Testing on cell 1_1:2 with center -0.7361215932 0.4250000000
+DEAL::Exact inverse:                   -0.3199797942 0.004618771511
+DEAL::Affine approximation:            -0.3011286728 2.212804616
+DEAL::Inverse quadratic approximation: -0.4188102773 0.01487792182
+DEAL::
+DEAL::Testing on cell 1_1:3 with center -0.7361215932 -0.4250000000
+DEAL::Exact inverse:                   -1.321431983 0.1389815030
+DEAL::Affine approximation:            -0.4994345174 12.41599603
+DEAL::Inverse quadratic approximation: -1.446711466 0.4916845141
+DEAL::
+DEAL::
+DEAL::Testing 2D with point -0.1779244106 0.9327124237
+DEAL::Testing on cell 0_1:0 with center 0.8227241336 0.4750000000
+DEAL::Exact inverse:                   1.679734435 0.5096451819
+DEAL::Affine approximation:            1.443910378 6.394226417
+DEAL::Inverse quadratic approximation: 1.823870032 0.6069628230
+DEAL::
+DEAL::Testing on cell 0_1:1 with center 5.817072296e-17 0.9500000000
+DEAL::Exact inverse:                   0.6800000317 0.5046876241
+DEAL::Affine approximation:            0.6872888532 -0.7700353780
+DEAL::Inverse quadratic approximation: 0.6735562329 0.5114552797
+DEAL::
+DEAL::Testing on cell 0_1:2 with center 0.7361215932 0.4250000000
+DEAL::Exact inverse:                   1.679734434 -0.4903548181
+DEAL::Affine approximation:            1.554958657 5.394226417
+DEAL::Inverse quadratic approximation: 1.973382985 -0.4493850397
+DEAL::
+DEAL::Testing on cell 0_1:3 with center 5.204748896e-17 0.8500000000
+DEAL::Exact inverse:                   0.6800000317 -0.4953123759
+DEAL::Affine approximation:            0.7093228359 -1.770035378
+DEAL::Inverse quadratic approximation: 0.6746285084 -0.5382929038
+DEAL::
+DEAL::Testing on cell 1_1:0 with center -0.8227241336 0.4750000000
+DEAL::Exact inverse:                   -0.3199797942 0.5048753699
+DEAL::Affine approximation:            -0.2566215243 2.835738205
+DEAL::Inverse quadratic approximation: -0.3823551506 0.5315544008
+DEAL::
+DEAL::Testing on cell 1_1:1 with center -0.8227241336 -0.4750000000
+DEAL::Exact inverse:                   -1.321431983 0.6467026977
+DEAL::Affine approximation:            -0.4439103775 13.60577358
+DEAL::Inverse quadratic approximation: -1.276004830 0.9290893899
+DEAL::
+DEAL::Testing on cell 1_1:2 with center -0.7361215932 0.4250000000
+DEAL::Exact inverse:                   -0.3199797942 -0.4951246301
+DEAL::Affine approximation:            -0.3456358213 1.835738205
+DEAL::Inverse quadratic approximation: -0.4448257264 -0.5335632508
+DEAL::
+DEAL::Testing on cell 1_1:3 with center -0.7361215932 -0.4250000000
+DEAL::Exact inverse:                   -1.321431982 -0.3532973076
+DEAL::Affine approximation:            -0.5549586572 12.60577358
+DEAL::Inverse quadratic approximation: -1.570577755 -0.03655836469
+DEAL::
+DEAL::
+DEAL::Testing 2D with point -0.1872888532 0.9818025513
+DEAL::Testing on cell 0_1:0 with center 0.8227241336 0.4750000000
+DEAL::Exact inverse:                   1.679734435 0.01015282302
+DEAL::Affine approximation:            1.493589871 6.204448860
+DEAL::Inverse quadratic approximation: 1.881022893 0.08910626957
+DEAL::
+DEAL::Testing on cell 0_1:1 with center 5.817072296e-17 0.9500000000
+DEAL::Exact inverse:                   0.6800000317 0.004934341186
+DEAL::Affine approximation:            0.6971461613 -1.336879345
+DEAL::Inverse quadratic approximation: 0.6752680942 0.001010067356
+DEAL::
+DEAL::Testing on cell 0_1:2 with center 0.7361215932 0.4250000000
+DEAL::Exact inverse:                   1.679734433 -0.9898471766
+DEAL::Affine approximation:            1.610482797 5.204448860
+DEAL::Inverse quadratic approximation: 2.034386667 -1.028496543
+DEAL::
+DEAL::Testing on cell 0_1:3 with center 5.204748896e-17 0.8500000000
+DEAL::Exact inverse:                   0.6800000317 -0.9950656588
+DEAL::Affine approximation:            0.7203398273 -2.336879345
+DEAL::Inverse quadratic approximation: 0.6740153620 -1.104288093
+DEAL::
+DEAL::Testing on cell 1_1:0 with center -0.8227241336 0.4750000000
+DEAL::Exact inverse:                   -0.3199797942 0.005131968346
+DEAL::Affine approximation:            -0.2964437098 2.458671795
+DEAL::Inverse quadratic approximation: -0.4088477085 0.01435738807
+DEAL::
+DEAL::Testing on cell 1_1:1 with center -0.8227241336 -0.4750000000
+DEAL::Exact inverse:                   -1.321431983 0.1544238922
+DEAL::Affine approximation:            -0.4935898711 13.79555114
+DEAL::Inverse quadratic approximation: -1.382003278 0.4281868663
+DEAL::
+DEAL::Testing on cell 1_1:2 with center -0.7361215932 0.4250000000
+DEAL::Exact inverse:                   -0.3199797940 -0.9948680322
+DEAL::Affine approximation:            -0.3901429698 1.458671795
+DEAL::Inverse quadratic approximation: -0.4682064843 -1.110483122
+DEAL::
+DEAL::Testing on cell 1_1:3 with center -0.7361215932 -0.4250000000
+DEAL::Exact inverse:                   -1.321431978 -0.8455762513
+DEAL::Affine approximation:            -0.6104827971 12.79555114
+DEAL::Inverse quadratic approximation: -1.696098323 -0.5939421483
+DEAL::
+DEAL::
+DEAL::Testing 3D with point -0.1498310826 0.7854420410 0.02512860726
+DEAL::Testing on cell 4_0: with center -1.334478803e-17 -0.9000000000 0.000000000
+DEAL::Exact inverse:                   0.6173118684 9.000405545 0.4801846852
+DEAL::Affine approximation:            0.3558249736 11.80212761 0.5241800136
+DEAL::Inverse quadratic approximation: 0.1534157163 1.178704386 0.5581266597
+DEAL::
+DEAL::Testing on cell 5_0: with center -5.358086410e-17 0.9000000000 0.000000000
+DEAL::Exact inverse:                   0.9995944549 0.3826881316 0.5198153148
+DEAL::Affine approximation:            -1.802127607 0.3558249736 0.5241800136
+DEAL::Inverse quadratic approximation: 0.9864508017 0.3984630214 0.5170290624
+DEAL::
+DEAL::
+DEAL::Testing 3D with point -0.1591955252 0.8345321686 0.02669914522
+DEAL::Testing on cell 4_0: with center -1.334478803e-17 -0.9000000000 0.000000000
+DEAL::Exact inverse:                   0.6173118684 9.250430892 0.4801846852
+DEAL::Affine approximation:            0.3468140344 12.22726058 0.5256912645
+DEAL::Inverse quadratic approximation: 0.1236178622 0.9629093195 0.5631241446
+DEAL::
+DEAL::Testing on cell 5_0: with center -5.358086410e-17 0.9000000000 0.000000000
+DEAL::Exact inverse:                   0.7495691083 0.3826881316 0.5198153148
+DEAL::Affine approximation:            -2.227260583 0.3468140344 0.5256912645
+DEAL::Inverse quadratic approximation: 0.7586398860 0.4002532966 0.5167288102
+DEAL::
+DEAL::
+DEAL::Testing 3D with point -0.1685599679 0.8836222961 0.02826968317
+DEAL::Testing on cell 4_0: with center -1.334478803e-17 -0.9000000000 0.000000000
+DEAL::Exact inverse:                   0.6173118684 9.500456238 0.4801846852
+DEAL::Affine approximation:            0.3378030953 12.65239356 0.5272025153
+DEAL::Inverse quadratic approximation: 0.09286279217 0.7336716472 0.5682821670
+DEAL::
+DEAL::Testing on cell 5_0: with center -5.358086410e-17 0.9000000000 0.000000000
+DEAL::Exact inverse:                   0.4995437617 0.3826881316 0.5198153148
+DEAL::Affine approximation:            -2.652393558 0.3378030953 0.5272025153
+DEAL::Inverse quadratic approximation: 0.5173863647 0.4030007878 0.5162680204
+DEAL::
+DEAL::
+DEAL::Testing 3D with point -0.1779244106 0.9327124237 0.02984022112
+DEAL::Testing on cell 4_0: with center -1.334478803e-17 -0.9000000000 0.000000000
+DEAL::Exact inverse:                   0.6173118684 9.750481585 0.4801846852
+DEAL::Affine approximation:            0.3287921561 13.07752653 0.5287137662
+DEAL::Inverse quadratic approximation: 0.06115050607 0.4909913695 0.5736007268
+DEAL::
+DEAL::Testing on cell 5_0: with center -5.358086410e-17 0.9000000000 0.000000000
+DEAL::Exact inverse:                   0.2495184151 0.3826881316 0.5198153148
+DEAL::Affine approximation:            -3.077526534 0.3287921561 0.5287137662
+DEAL::Inverse quadratic approximation: 0.2626902379 0.4067054950 0.5156466932
+DEAL::
+DEAL::
+DEAL::Testing 3D with point -0.1872888532 0.9818025513 0.03141075908
+DEAL::Testing on cell 4_0: with center -1.334478803e-17 -0.9000000000 0.000000000
+DEAL::Exact inverse:                   0.6173118684 10.00050693 0.4801846852
+DEAL::Affine approximation:            0.3197812170 13.50265951 0.5302250170
+DEAL::Inverse quadratic approximation: 0.02848100392 0.2348684862 0.5790798242
+DEAL::
+DEAL::Testing on cell 5_0: with center -5.358086410e-17 0.9000000000 0.000000000
+DEAL::Exact inverse:                   -0.0005069314277 0.3826881316 0.5198153148
+DEAL::Affine approximation:            -3.502659509 0.3197812170 0.5302250170
+DEAL::Inverse quadratic approximation: -0.005448494407 0.4113674182 0.5148648284
+DEAL::
+DEAL::
+DEAL::Testing 3D with point 0.04694448799 0.06461354454 0.7960033322
+DEAL::Testing on cell 0_0: with center -9.723960037e-18 -9.723960037e-18 -0.9000000000
+DEAL::Exact inverse:                   0.4633821073 0.4496318287 9.000130204
+DEAL::Affine approximation:            0.5451723546 0.5621744122 11.89359107
+DEAL::Inverse quadratic approximation: 0.6091066156 0.6501723730 1.181455366
+DEAL::
+DEAL::Testing on cell 2_0: with center -1.696561603e-17 0.000000000 0.9000000000
+DEAL::Exact inverse:                   0.5366178927 0.9998697959 0.5503681713
+DEAL::Affine approximation:            0.5451723546 -1.893591072 0.5621744122
+DEAL::Inverse quadratic approximation: 0.5312969833 0.9866166821 0.5430766020
+DEAL::
+DEAL::
+DEAL::Testing 3D with point 0.04890050833 0.06730577556 0.8291701377
+DEAL::Testing on cell 0_0: with center -9.723960037e-18 -9.723960037e-18 -0.9000000000
+DEAL::Exact inverse:                   0.4633821073 0.4496318287 9.166802296
+DEAL::Affine approximation:            0.5470545361 0.5647650127 12.18082403
+DEAL::Inverse quadratic approximation: 0.6153413016 0.6587536822 1.039263536
+DEAL::
+DEAL::Testing on cell 2_0: with center -1.696561603e-17 0.000000000 0.9000000000
+DEAL::Exact inverse:                   0.5366178927 0.8331977041 0.5503681713
+DEAL::Affine approximation:            0.5470545361 -2.180824033 0.5647650127
+DEAL::Inverse quadratic approximation: 0.5309124472 0.8363065727 0.5425473335
+DEAL::
+DEAL::
+DEAL::Testing 3D with point 0.05085652866 0.06999800658 0.8623369432
+DEAL::Testing on cell 0_0: with center -9.723960037e-18 -9.723960037e-18 -0.9000000000
+DEAL::Exact inverse:                   0.4633821073 0.4496318287 9.333474387
+DEAL::Affine approximation:            0.5489367175 0.5673556132 12.46805699
+DEAL::Inverse quadratic approximation: 0.6217110739 0.6675209216 0.8911022773
+DEAL::
+DEAL::Testing on cell 2_0: with center -1.696561603e-17 0.000000000 0.9000000000
+DEAL::Exact inverse:                   0.5366178927 0.6665256123 0.5503681713
+DEAL::Affine approximation:            0.5489367175 -2.468056995 0.5673556132
+DEAL::Inverse quadratic approximation: 0.5303928249 0.6800270360 0.5418321347
+DEAL::
+DEAL::
+DEAL::Testing 3D with point 0.05281254899 0.07269023761 0.8955037488
+DEAL::Testing on cell 0_0: with center -9.723960037e-18 -9.723960037e-18 -0.9000000000
+DEAL::Exact inverse:                   0.4633821073 0.4496318287 9.500146479
+DEAL::Affine approximation:            0.5508188990 0.5699462137 12.75528996
+DEAL::Inverse quadratic approximation: 0.6282159323 0.6764740911 0.7369715917
+DEAL::
+DEAL::Testing on cell 2_0: with center -1.696561603e-17 0.000000000 0.9000000000
+DEAL::Exact inverse:                   0.5366178927 0.4998535204 0.5503681713
+DEAL::Affine approximation:            0.5508188990 -2.755289956 0.5699462137
+DEAL::Inverse quadratic approximation: 0.5297381165 0.5177780719 0.5409310059
+DEAL::
+DEAL::
+DEAL::Testing 3D with point 0.05476856932 0.07538246863 0.9286705543
+DEAL::Testing on cell 0_0: with center -9.723960037e-18 -9.723960037e-18 -0.9000000000
+DEAL::Exact inverse:                   0.4633821073 0.4496318287 9.666818571
+DEAL::Affine approximation:            0.5527010804 0.5725368143 13.04252292
+DEAL::Inverse quadratic approximation: 0.6348558768 0.6856131907 0.5768714788
+DEAL::
+DEAL::Testing on cell 2_0: with center -1.696561603e-17 0.000000000 0.9000000000
+DEAL::Exact inverse:                   0.5366178927 0.3331814286 0.5503681713
+DEAL::Affine approximation:            0.5527010804 -3.042522917 0.5725368143
+DEAL::Inverse quadratic approximation: 0.5289483219 0.3495596805 0.5398439468
+DEAL::
+DEAL::
+DEAL::Testing 3D with point 0.05672458966 0.07807469965 0.9618373598
+DEAL::Testing on cell 0_0: with center -9.723960037e-18 -9.723960037e-18 -0.9000000000
+DEAL::Exact inverse:                   0.4633821073 0.4496318287 9.833490663
+DEAL::Affine approximation:            0.5545832618 0.5751274148 13.32975588
+DEAL::Inverse quadratic approximation: 0.6416309076 0.6949382206 0.4108019385
+DEAL::
+DEAL::Testing on cell 2_0: with center -1.696561603e-17 0.000000000 0.9000000000
+DEAL::Exact inverse:                   0.5366178927 0.1665093367 0.5503681713
+DEAL::Affine approximation:            0.5545832618 -3.329755879 0.5751274148
+DEAL::Inverse quadratic approximation: 0.5280234411 0.1753718617 0.5385709576
+DEAL::
+DEAL::
+DEAL::Testing 3D with point 0.05868060999 0.08076693067 0.9950041653
+DEAL::Testing on cell 0_0: with center -9.723960037e-18 -9.723960037e-18 -0.9000000000
+DEAL::Exact inverse:                   0.4633821073 0.4496318287 10.00016275
+DEAL::Affine approximation:            0.5564654433 0.5777180153 13.61698884
+DEAL::Inverse quadratic approximation: 0.6485410245 0.7044491805 0.2387629709
+DEAL::
+DEAL::Testing on cell 2_0: with center -1.696561603e-17 0.000000000 0.9000000000
+DEAL::Exact inverse:                   0.5366178927 -0.0001627550900 0.5503681713
+DEAL::Affine approximation:            0.5564654433 -3.616988840 0.5777180153
+DEAL::Inverse quadratic approximation: 0.5269634741 -0.004785384458 0.5371120383
+DEAL::
+DEAL::

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.