]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement class NonMatching::MappingInfo 13835/head
authorMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Thu, 19 May 2022 11:58:06 +0000 (13:58 +0200)
committerMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Thu, 2 Jun 2022 12:17:07 +0000 (14:17 +0200)
include/deal.II/matrix_free/fe_point_evaluation.h
include/deal.II/non_matching/mapping_info.h [new file with mode: 0644]
tests/matrix_free/point_evaluation_15.cc
tests/matrix_free/point_evaluation_16.cc
tests/matrix_free/point_evaluation_17.cc [new file with mode: 0644]
tests/matrix_free/point_evaluation_17.output [new file with mode: 0644]

index e1f7e5e474c6692e4699d1041cf62da6842205ab..c722e0e1586f91a1a61dad915eef1b4e8c4945aa 100644 (file)
@@ -33,6 +33,8 @@
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/tensor_product_kernels.h>
 
+#include <deal.II/non_matching/mapping_info.h>
+
 DEAL_II_NAMESPACE_OPEN
 
 namespace internal
@@ -364,48 +366,6 @@ namespace internal
     template <int dim, int spacedim>
     std::vector<Polynomials::Polynomial<double>>
     get_polynomial_space(const FiniteElement<dim, spacedim> &fe);
-
-    /**
-     * Compute the mapping related data for the given @p mapping,
-     * @p cell and @p unit_points that is required by the FEPointEvaluation
-     * class.
-     */
-    template <int dim, int spacedim>
-    void
-    compute_mapping_data_for_generic_points(
-      const Mapping<dim> &                                        mapping,
-      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-      const ArrayView<const Point<dim>> &                         unit_points,
-      const UpdateFlags                                           update_flags,
-      internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
-        &mapping_data)
-    {
-      UpdateFlags update_flags_mapping = update_default;
-      // translate update flags
-      if (update_flags & update_jacobians)
-        update_flags_mapping |= update_jacobians;
-      if (update_flags & update_gradients ||
-          update_flags & update_inverse_jacobians)
-        update_flags_mapping |= update_inverse_jacobians;
-      if (update_flags & update_quadrature_points)
-        update_flags_mapping |= update_quadrature_points;
-
-      if (const MappingQ<dim, spacedim> *mapping_q =
-            dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
-        {
-          mapping_q->fill_mapping_data_for_generic_points(cell,
-                                                          unit_points,
-                                                          update_flags_mapping,
-                                                          mapping_data);
-        }
-      else if (const MappingCartesian<dim, spacedim> *mapping_cartesian =
-                 dynamic_cast<const MappingCartesian<dim, spacedim> *>(
-                   &mapping))
-        {
-          mapping_cartesian->fill_mapping_data_for_generic_points(
-            cell, unit_points, update_flags_mapping, mapping_data);
-        }
-    }
   } // namespace FEPointEvaluation
 } // namespace internal
 
@@ -476,55 +436,39 @@ public:
                     const unsigned int        first_selected_component = 0);
 
   /**
-   * Set up the mapping information for the given cell, e.g., by computing the
-   * Jacobian of the mapping for the given points if gradients of the functions
-   * are requested.
+   * Constructor to make the present class able to re-use the geometry
+   * data also used by other `FEPointEvaluation` objects.
    *
-   * @param[in] cell An iterator to the current cell
+   * @param mapping_info The MappingInfo class describes the geometry-related
+   * data for evaluating finite-element solutions. This object enables to
+   * construct such an object on the outside, possibly re-using it between
+   * several objects or between several calls to the same cell and unit points.
    *
-   * @param[in] unit_points List of points in the reference locations of the
-   * current cell where the FiniteElement object should be
-   * evaluated/integrated in the evaluate() and integrate() functions.
+   * @param fe The FiniteElement object that is used for the evaluation, which
+   * is typically the same on all cells to be evaluated.
+   *
+   * @param first_selected_component For multi-component FiniteElement
+   * objects, this parameter allows to select a range of `n_components`
+   * components starting from this parameter.
    */
-  void
-  reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-         const ArrayView<const Point<dim>> &unit_points);
+  FEPointEvaluation(NonMatching::MappingInfo<dim, spacedim> &mapping_info,
+                    const FiniteElement<dim> &               fe,
+                    const unsigned int first_selected_component = 0);
 
   /**
-   * Set up the mapping information for the given cell. This function is
-   * an alternative to the function with the same name above and uses
-   * a precomputed @p mapping_data object. This function can be used
-   * to avoid duplicated evaluation of the mapping if multiple
-   * FEPointEvaluation objects for the different components of the same FESystem
-   * are used. You can get the mapping data from an initialized
-   * FEPointEvaluation object by calling the function get_mapping_data().
+   * Set up the mapping information for the given cell, e.g., by computing the
+   * Jacobian of the mapping for the given points if gradients of the functions
+   * are requested.
    *
    * @param[in] cell An iterator to the current cell
    *
    * @param[in] unit_points List of points in the reference locations of the
    * current cell where the FiniteElement object should be
    * evaluated/integrated in the evaluate() and integrate() functions.
-   *
-   * @param[in] mapping_data A mapping data object that is precomputed
-   * for the given cell and positions.
    */
   void
   reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-         const ArrayView<const Point<dim>> &unit_points,
-         const dealii::internal::FEValuesImplementation::
-           MappingRelatedData<dim, spacedim> &mapping_data);
-
-  /**
-   * Returns the mapping data that was computed during the last call to
-   * the reinit() function. This can be useful if multiple FEPointEvaluation
-   * objects are used for multiple components of a FESystem. The mapping data
-   * can be retrieved from the first FEPointEvaluation object and subsequently
-   * passed to the other objects, avoiding the need to recompute the mapping
-   * data.
-   */
-  const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
-                                                                     spacedim> &
-  get_mapping_data() const;
+         const ArrayView<const Point<dim>> &unit_points);
 
   /**
    * This function interpolates the finite element solution, represented by
@@ -656,6 +600,17 @@ public:
   unit_point(const unsigned int point_index) const;
 
 private:
+  /**
+   * Common setup function for both constructors. Does the setup for both fast
+   * and slow path.
+   *
+   * @param first_selected_component For multi-component FiniteElement
+   * objects, this parameter allows to select a range of `n_components`
+   * components starting from this parameter.
+   */
+  void
+  setup(const unsigned int first_selected_component);
+
   /**
    * Pointer to the Mapping object passed to the constructor.
    */
@@ -739,29 +694,34 @@ private:
   /**
    * The desired update flags for the evaluation.
    */
-  UpdateFlags update_flags;
+  const UpdateFlags update_flags;
 
   /**
-   * The update flags specific for the mapping in the fast evaluation path.
+   * The FEValues object underlying the slow evaluation path.
    */
-  UpdateFlags update_flags_mapping;
+  std::shared_ptr<FEValues<dim, spacedim>> fe_values;
 
   /**
-   * The FEValues object underlying the slow evaluation path.
+   * Pointer to mapping info on the fly computed during reinit.
    */
-  std::shared_ptr<FEValues<dim, spacedim>> fe_values;
+  std::unique_ptr<NonMatching::MappingInfo<dim, spacedim>>
+    mapping_info_on_the_fly;
 
   /**
-   * Array to store temporary data computed by the mapping for the fast
-   * evaluation path.
+   * Pointer to currently used mapping info (either on the fly or external
+   * precomputed).
    */
-  internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
-    mapping_data;
+  SmartPointer<NonMatching::MappingInfo<dim, spacedim>> mapping_info;
 
   /**
    * The reference points specified at reinit().
    */
   std::vector<Point<dim>> unit_points;
+
+  /**
+   * Bool indicating if fast path is chosen.
+   */
+  bool fast_path;
 };
 
 // ----------------------- template and inline function ----------------------
@@ -776,73 +736,97 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
   : mapping(&mapping)
   , fe(&fe)
   , update_flags(update_flags)
-  , update_flags_mapping(update_default)
+  , mapping_info_on_the_fly(
+      std::make_unique<NonMatching::MappingInfo<dim, spacedim>>(mapping,
+                                                                update_flags))
+  , mapping_info(mapping_info_on_the_fly.get())
+{
+  setup(first_selected_component);
+}
+
+
+
+template <int n_components, int dim, int spacedim, typename Number>
+FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
+  NonMatching::MappingInfo<dim, spacedim> &mapping_info,
+  const FiniteElement<dim> &               fe,
+  const unsigned int                       first_selected_component)
+  : mapping(&mapping_info.get_mapping())
+  , fe(&fe)
+  , update_flags(mapping_info.get_update_flags())
+  , mapping_info(&mapping_info)
+{
+  setup(first_selected_component);
+}
+
+
+
+template <int n_components, int dim, int spacedim, typename Number>
+void
+FEPointEvaluation<n_components, dim, spacedim, Number>::setup(
+  const unsigned int first_selected_component)
 {
   AssertIndexRange(first_selected_component + n_components,
-                   fe.n_components() + 1);
+                   fe->n_components() + 1);
 
   bool         same_base_element   = true;
   unsigned int base_element_number = 0;
   component_in_base_element        = 0;
   unsigned int component           = 0;
-  for (; base_element_number < fe.n_base_elements(); ++base_element_number)
-    if (component + fe.element_multiplicity(base_element_number) >
+  for (; base_element_number < fe->n_base_elements(); ++base_element_number)
+    if (component + fe->element_multiplicity(base_element_number) >
         first_selected_component)
       {
         if (first_selected_component + n_components >
-            component + fe.element_multiplicity(base_element_number))
+            component + fe->element_multiplicity(base_element_number))
           same_base_element = false;
         component_in_base_element = first_selected_component - component;
         break;
       }
     else
-      component += fe.element_multiplicity(base_element_number);
+      component += fe->element_multiplicity(base_element_number);
 
-  if (internal::FEPointEvaluation::is_fast_path_supported(mapping) &&
+  if (internal::FEPointEvaluation::is_fast_path_supported(*mapping) &&
       internal::FEPointEvaluation::is_fast_path_supported(
-        fe, base_element_number) &&
+        *fe, base_element_number) &&
       same_base_element)
     {
       internal::MatrixFreeFunctions::ShapeInfo<double> shape_info;
 
-      shape_info.reinit(QMidpoint<1>(), fe, base_element_number);
+      shape_info.reinit(QMidpoint<1>(), *fe, base_element_number);
       renumber           = shape_info.lexicographic_numbering;
       dofs_per_component = shape_info.dofs_per_component_on_cell;
       poly               = internal::FEPointEvaluation::get_polynomial_space(
-        fe.base_element(base_element_number));
+        fe->base_element(base_element_number));
 
       polynomials_are_hat_functions =
         (poly.size() == 2 && poly[0].value(0.) == 1. &&
          poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
          poly[1].value(1.) == 1.);
+
+      fast_path = true;
     }
   else
     {
-      nonzero_shape_function_component.resize(fe.n_dofs_per_cell());
+      nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
       for (unsigned int d = 0; d < n_components; ++d)
         {
           const unsigned int component = first_selected_component + d;
-          for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
+          for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
             {
-              const bool is_primitive = fe.is_primitive() || fe.is_primitive(i);
+              const bool is_primitive =
+                fe->is_primitive() || fe->is_primitive(i);
               if (is_primitive)
                 nonzero_shape_function_component[i][d] =
-                  (component == fe.system_to_component_index(i).first);
+                  (component == fe->system_to_component_index(i).first);
               else
                 nonzero_shape_function_component[i][d] =
-                  (fe.get_nonzero_components(i)[component] == true);
+                  (fe->get_nonzero_components(i)[component] == true);
             }
         }
-    }
 
-  // translate update flags
-  if (update_flags & update_jacobians)
-    update_flags_mapping |= update_jacobians;
-  if (update_flags & update_gradients ||
-      update_flags & update_inverse_jacobians)
-    update_flags_mapping |= update_inverse_jacobians;
-  if (update_flags & update_quadrature_points)
-    update_flags_mapping |= update_quadrature_points;
+      fast_path = false;
+    }
 }
 
 
@@ -853,73 +837,25 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
   const typename Triangulation<dim, spacedim>::cell_iterator &cell,
   const ArrayView<const Point<dim>> &                         unit_points)
 {
-  // If using the fast path, we need to precompute the mapping data.
-  if (!poly.empty())
-    internal::FEPointEvaluation::compute_mapping_data_for_generic_points(
-      *mapping, cell, unit_points, update_flags_mapping, mapping_data);
-
-  // then call the other version of this function with the precomputed data
-  reinit(cell, unit_points, mapping_data);
-}
-
+  // reinit is only allowed for mapping computation on the fly
+  AssertThrow(mapping_info_on_the_fly.get() != nullptr, ExcNotImplemented());
 
+  mapping_info->reinit(cell, unit_points);
 
-template <int n_components, int dim, int spacedim, typename Number>
-void
-FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
-  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-  const ArrayView<const Point<dim>> &                         unit_points,
-  const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
-                                                                     spacedim>
-    &precomputed_mapping_data)
-{
-  this->unit_points.resize(unit_points.size());
-  std::copy(unit_points.begin(), unit_points.end(), this->unit_points.begin());
-
-  if (!poly.empty())
-    {
-      // Check the mapping data for consistency.
-      if ((update_flags_mapping & update_jacobians) != 0u)
-        Assert(precomputed_mapping_data.jacobians.size() == unit_points.size(),
-               ExcDimensionMismatch(precomputed_mapping_data.jacobians.size(),
-                                    unit_points.size()));
-      if ((update_flags_mapping & update_inverse_jacobians) != 0u)
-        Assert(precomputed_mapping_data.inverse_jacobians.size() ==
-                 unit_points.size(),
-               ExcDimensionMismatch(
-                 precomputed_mapping_data.inverse_jacobians.size(),
-                 unit_points.size()));
-
-      if ((update_flags_mapping & update_quadrature_points) != 0u)
-        Assert(precomputed_mapping_data.inverse_jacobians.size() ==
-                 unit_points.size(),
-               ExcDimensionMismatch(
-                 precomputed_mapping_data.quadrature_points.size(),
-                 unit_points.size()));
-
-      mapping_data = precomputed_mapping_data;
-    }
-  else
+  if (!fast_path)
     {
       fe_values = std::make_shared<FEValues<dim, spacedim>>(
         *mapping,
         *fe,
         Quadrature<dim>(
           std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
-        update_flags | update_flags_mapping);
+        update_flags);
       fe_values->reinit(cell);
-      mapping_data.initialize(unit_points.size(), update_flags_mapping);
-      if (update_flags_mapping & update_jacobians)
-        for (unsigned int q = 0; q < unit_points.size(); ++q)
-          mapping_data.jacobians[q] = fe_values->jacobian(q);
-      if (update_flags_mapping & update_inverse_jacobians)
-        for (unsigned int q = 0; q < unit_points.size(); ++q)
-          mapping_data.inverse_jacobians[q] = fe_values->inverse_jacobian(q);
-      if (update_flags_mapping & update_quadrature_points)
-        for (unsigned int q = 0; q < unit_points.size(); ++q)
-          mapping_data.quadrature_points[q] = fe_values->quadrature_point(q);
     }
 
+  this->unit_points =
+    std::vector<Point<dim>>(unit_points.begin(), unit_points.end());
+
   if (update_flags & update_values)
     values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
   if (update_flags & update_gradients)
@@ -929,29 +865,31 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
 
 
 
-template <int n_components, int dim, int spacedim, typename Number>
-const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
-                                                                   spacedim> &
-FEPointEvaluation<n_components, dim, spacedim, Number>::get_mapping_data() const
-{
-  return mapping_data;
-}
-
-
-
 template <int n_components, int dim, int spacedim, typename Number>
 void
 FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
   const ArrayView<const Number> &         solution_values,
   const EvaluationFlags::EvaluationFlags &evaluation_flag)
 {
+  const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr;
+  if (precomputed_mapping)
+    {
+      unit_points = mapping_info->get_unit_points();
+
+      if (update_flags & update_values)
+        values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
+      if (update_flags & update_gradients)
+        gradients.resize(unit_points.size(),
+                         numbers::signaling_nan<gradient_type>());
+    }
+
   if (unit_points.empty())
     return;
 
   AssertDimension(solution_values.size(), fe->dofs_per_cell);
   if (((evaluation_flag & EvaluationFlags::values) ||
        (evaluation_flag & EvaluationFlags::gradients)) &&
-      !poly.empty())
+      fast_path)
     {
       // fast path with tensor product evaluation
       if (solution_renumbered.size() != dofs_per_component)
@@ -1002,20 +940,21 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
                      ExcNotInitialized());
               for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
                 {
-                  Assert(update_flags_mapping & update_inverse_jacobians,
-                         ExcNotInitialized());
                   internal::FEPointEvaluation::EvaluatorTypeTraits<
                     dim,
                     n_components,
                     Number>::set_gradient(val_and_grad.second,
                                           j,
                                           unit_gradients[i + j]);
-                  gradients[i + j] = static_cast<
-                    typename internal::FEPointEvaluation::
-                      EvaluatorTypeTraits<dim, n_components, double>::
-                        gradient_type>(apply_transformation(
-                    mapping_data.inverse_jacobians[i + j].transpose(),
-                    unit_gradients[i + j]));
+                  gradients[i + j] =
+                    static_cast<typename internal::FEPointEvaluation::
+                                  EvaluatorTypeTraits<dim,
+                                                      n_components,
+                                                      double>::gradient_type>(
+                      apply_transformation(mapping_info->get_mapping_data()
+                                             .inverse_jacobians[i + j]
+                                             .transpose(),
+                                           unit_gradients[i + j]));
                 }
             }
         }
@@ -1084,6 +1023,18 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
   const ArrayView<Number> &               solution_values,
   const EvaluationFlags::EvaluationFlags &integration_flags)
 {
+  const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr;
+  if (precomputed_mapping)
+    {
+      unit_points = mapping_info->get_unit_points();
+
+      if (update_flags & update_values)
+        values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
+      if (update_flags & update_gradients)
+        gradients.resize(unit_points.size(),
+                         numbers::signaling_nan<gradient_type>());
+    }
+
   if (unit_points.size() == 0) // no evaluation points provided
     {
       std::fill(solution_values.begin(), solution_values.end(), 0.0);
@@ -1093,7 +1044,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
   AssertDimension(solution_values.size(), fe->dofs_per_cell);
   if (((integration_flags & EvaluationFlags::values) ||
        (integration_flags & EvaluationFlags::gradients)) &&
-      !poly.empty())
+      fast_path)
     {
       // fast path with tensor product integration
 
@@ -1139,14 +1090,12 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
           if (integration_flags & EvaluationFlags::gradients)
             for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
               {
-                Assert(update_flags_mapping & update_inverse_jacobians,
-                       ExcNotInitialized());
                 gradients[i + j] =
                   static_cast<typename internal::FEPointEvaluation::
                                 EvaluatorTypeTraits<dim, n_components, double>::
-                                  gradient_type>(
-                    apply_transformation(mapping_data.inverse_jacobians[i + j],
-                                         gradients[i + j]));
+                                  gradient_type>(apply_transformation(
+                    mapping_info->get_mapping_data().inverse_jacobians[i + j],
+                    gradients[i + j]));
                 internal::FEPointEvaluation::
                   EvaluatorTypeTraits<dim, n_components, Number>::get_gradient(
                     gradient, j, gradients[i + j]);
@@ -1269,7 +1218,7 @@ inline const typename FEPointEvaluation<n_components, dim, spacedim, Number>::
   FEPointEvaluation<n_components, dim, spacedim, Number>::get_unit_gradient(
     const unsigned int point_index) const
 {
-  Assert(!poly.empty(),
+  Assert(fast_path,
          ExcMessage("Unit gradients are currently only implemented for tensor "
                     "product finite elements combined with MappingQ "
                     "mappings"));
@@ -1308,9 +1257,9 @@ inline DerivativeForm<1, dim, spacedim>
 FEPointEvaluation<n_components, dim, spacedim, Number>::jacobian(
   const unsigned int point_index) const
 {
-  Assert(update_flags_mapping & update_jacobians, ExcNotInitialized());
-  AssertIndexRange(point_index, mapping_data.jacobians.size());
-  return mapping_data.jacobians[point_index];
+  AssertIndexRange(point_index,
+                   mapping_info->get_mapping_data().jacobians.size());
+  return mapping_info->get_mapping_data().jacobians[point_index];
 }
 
 
@@ -1320,11 +1269,9 @@ inline DerivativeForm<1, spacedim, dim>
 FEPointEvaluation<n_components, dim, spacedim, Number>::inverse_jacobian(
   const unsigned int point_index) const
 {
-  Assert(update_flags_mapping & update_inverse_jacobians ||
-           update_flags_mapping & update_gradients,
-         ExcNotInitialized());
-  AssertIndexRange(point_index, mapping_data.inverse_jacobians.size());
-  return mapping_data.inverse_jacobians[point_index];
+  AssertIndexRange(point_index,
+                   mapping_info->get_mapping_data().inverse_jacobians.size());
+  return mapping_info->get_mapping_data().inverse_jacobians[point_index];
 }
 
 
@@ -1334,9 +1281,9 @@ inline Point<spacedim>
 FEPointEvaluation<n_components, dim, spacedim, Number>::real_point(
   const unsigned int point_index) const
 {
-  Assert(update_flags_mapping & update_quadrature_points, ExcNotInitialized());
-  AssertIndexRange(point_index, mapping_data.quadrature_points.size());
-  return mapping_data.quadrature_points[point_index];
+  AssertIndexRange(point_index,
+                   mapping_info->get_mapping_data().quadrature_points.size());
+  return mapping_info->get_mapping_data().quadrature_points[point_index];
 }
 
 
diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h
new file mode 100644 (file)
index 0000000..1a70680
--- /dev/null
@@ -0,0 +1,252 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2011 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_non_matching_mapping_info_h
+#define dealii_non_matching_mapping_info_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/aligned_vector.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/vectorization.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_update_flags.h>
+#include <deal.II/fe/mapping_cartesian.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <memory>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace NonMatching
+{
+  /**
+   * This class provides the mapping information computation and mapping data
+   * storage to be used together with FEPointEvaluation.
+   */
+  template <int dim, int spacedim = dim>
+  class MappingInfo : public Subscriptor
+  {
+  public:
+    /**
+     * Constructor.
+     *
+     * @param mapping The Mapping class describing the geometry of a cell.
+     *
+     * @param update_flags Specify the quantities to be computed by the mapping
+     * during the call of reinit(). These update flags are also handed to a
+     * FEEvaluation object if you construct it with this MappingInfo object.
+     */
+    MappingInfo(const Mapping<dim> &mapping, const UpdateFlags update_flags);
+
+    /**
+     * Reinitialize the mapping information for the incoming cell and unit
+     * points.
+     */
+    void
+    reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+           const ArrayView<const Point<dim>> &unit_points);
+
+    /**
+     * Getter function for current unit points.
+     */
+    const std::vector<Point<dim>> &
+    get_unit_points() const;
+
+    /**
+     * Getter function for computed mapping data. This function accesses
+     * internal data and is therefore not a stable interface.
+     */
+    const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
+                                                                       spacedim>
+      &
+      get_mapping_data() const;
+
+    /**
+     * Getter function for underlying mapping.
+     */
+    const Mapping<dim, spacedim> &
+    get_mapping() const;
+
+    /**
+     * Getter function for the update flags.
+     */
+    UpdateFlags
+    get_update_flags() const;
+
+  private:
+    /**
+     * Compute the mapping related data for the given @p mapping,
+     * @p cell and @p unit_points that is required by the FEPointEvaluation
+     * class.
+     */
+    void
+    compute_mapping_data_for_generic_points(
+      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+      const ArrayView<const Point<dim>> &                         unit_points);
+
+    /**
+     * The reference points specified at reinit().
+     */
+    std::vector<Point<dim>> unit_points;
+
+    /**
+     * A pointer to the underlying mapping.
+     */
+    const SmartPointer<const Mapping<dim, spacedim>> mapping;
+
+    /**
+     * The desired update flags for the evaluation.
+     */
+    const UpdateFlags update_flags;
+
+    /**
+     * The update flags for the desired mapping information.
+     */
+    UpdateFlags update_flags_mapping;
+
+    /**
+     * The internal data container for mapping information. The implementation
+     * is subject to future changes.
+     */
+    dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+      mapping_data;
+  };
+
+  // ----------------------- template functions ----------------------
+
+
+  template <int dim, int spacedim>
+  MappingInfo<dim, spacedim>::MappingInfo(const Mapping<dim> &mapping,
+                                          const UpdateFlags   update_flags)
+    : mapping(&mapping)
+    , update_flags(update_flags)
+  {
+    update_flags_mapping = update_default;
+    // translate update flags
+    if (update_flags & update_jacobians)
+      update_flags_mapping |= update_jacobians;
+    if (update_flags & update_gradients ||
+        update_flags & update_inverse_jacobians)
+      update_flags_mapping |= update_inverse_jacobians;
+
+    // always save quadrature points for now
+    update_flags_mapping |= update_quadrature_points;
+  }
+
+
+
+  template <int dim, int spacedim>
+  void
+  MappingInfo<dim, spacedim>::reinit(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const ArrayView<const Point<dim>> &                         unit_points)
+  {
+    this->unit_points =
+      std::vector<Point<dim>>(unit_points.begin(), unit_points.end());
+    compute_mapping_data_for_generic_points(cell, unit_points);
+  }
+
+
+
+  template <int dim, int spacedim>
+  const std::vector<Point<dim>> &
+  MappingInfo<dim, spacedim>::get_unit_points() const
+  {
+    return unit_points;
+  }
+
+
+
+  template <int dim, int spacedim>
+  const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
+                                                                     spacedim> &
+  MappingInfo<dim, spacedim>::get_mapping_data() const
+  {
+    return mapping_data;
+  }
+
+
+
+  template <int dim, int spacedim>
+  const Mapping<dim, spacedim> &
+  MappingInfo<dim, spacedim>::get_mapping() const
+  {
+    return *mapping;
+  }
+
+
+
+  template <int dim, int spacedim>
+  UpdateFlags
+  MappingInfo<dim, spacedim>::get_update_flags() const
+  {
+    return update_flags;
+  }
+
+
+
+  template <int dim, int spacedim>
+  void
+  MappingInfo<dim, spacedim>::compute_mapping_data_for_generic_points(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const ArrayView<const Point<dim>> &                         unit_points)
+  {
+    if (const MappingQ<dim, spacedim> *mapping_q =
+          dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+      {
+        mapping_q->fill_mapping_data_for_generic_points(cell,
+                                                        unit_points,
+                                                        update_flags_mapping,
+                                                        mapping_data);
+      }
+    else if (const MappingCartesian<dim, spacedim> *mapping_cartesian =
+               dynamic_cast<const MappingCartesian<dim, spacedim> *>(
+                 &(*mapping)))
+      {
+        mapping_cartesian->fill_mapping_data_for_generic_points(
+          cell, unit_points, update_flags_mapping, mapping_data);
+      }
+    else
+      {
+        FE_DGQ<dim, spacedim>           dummy_fe(1);
+        dealii::FEValues<dim, spacedim> fe_values(
+          *mapping,
+          dummy_fe,
+          Quadrature<dim>(
+            std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
+          update_flags_mapping);
+        fe_values.reinit(cell);
+        mapping_data.initialize(unit_points.size(), update_flags_mapping);
+        if (update_flags_mapping & update_jacobians)
+          for (unsigned int q = 0; q < unit_points.size(); ++q)
+            mapping_data.jacobians[q] = fe_values.jacobian(q);
+        if (update_flags_mapping & update_inverse_jacobians)
+          for (unsigned int q = 0; q < unit_points.size(); ++q)
+            mapping_data.inverse_jacobians[q] = fe_values.inverse_jacobian(q);
+        if (update_flags_mapping & update_quadrature_points)
+          for (unsigned int q = 0; q < unit_points.size(); ++q)
+            mapping_data.quadrature_points[q] = fe_values.quadrature_point(q);
+      }
+  }
+} // namespace NonMatching
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 18c4476e264d6a85a00bea81a277ba9bace96b3d..8dc2df4be12edc2cfc622fc0100494bf3b1ecdb4 100644 (file)
@@ -84,16 +84,15 @@ test()
   dof_handler.distribute_dofs(fe);
   Vector<double> vector(dof_handler.n_dofs());
 
-  // Use the first evaluator to compute the values of the function
-  FEPointEvaluation<dim, dim> evaluator1(mapping,
-                                         fe,
-                                         update_values | update_gradients);
+  // mapping info object for precomputed mapping
+  NonMatching::MappingInfo<dim, dim> mapping_info(mapping,
+                                                  update_values |
+                                                    update_gradients);
+  // First evaluator to use precomputed mapping
+  FEPointEvaluation<dim, dim> evaluator1(mapping_info, fe);
 
-  // The second evaluator will use the precomputed mapping information
-  // from evaluator1 to compute the values of the function
-  FEPointEvaluation<dim, dim> evaluator2(mapping,
-                                         fe,
-                                         update_values | update_gradients);
+  // Second evaluator to use precomputed mapping
+  FEPointEvaluation<dim, dim> evaluator2(mapping_info, fe);
 
   VectorTools::interpolate(mapping, dof_handler, MyFunction<dim>(), vector);
 
@@ -108,11 +107,10 @@ test()
                            solution_values.begin(),
                            solution_values.end());
 
-      evaluator1.reinit(cell, unit_points);
+      mapping_info.reinit(cell, unit_points);
+
       evaluator1.evaluate(solution_values,
                           EvaluationFlags::values | EvaluationFlags::gradients);
-
-      evaluator2.reinit(cell, unit_points, evaluator1.get_mapping_data());
       evaluator2.evaluate(solution_values,
                           EvaluationFlags::values | EvaluationFlags::gradients);
 
index 9f9fbaaa2ef781025474fb6b81706a99a6420828..e2c12e07660908c530c5f7f1d86191db7436644e 100644 (file)
@@ -88,12 +88,12 @@ test(const unsigned int degree)
   dof_handler.distribute_dofs(fe);
   Vector<double> vector(dof_handler.n_dofs());
 
-  FEPointEvaluation<dim, dim> evaluator(mapping,
-                                        fe,
-                                        update_values | update_gradients);
-  FEPointEvaluation<dim, dim> evaluator2(mapping,
-                                         fe,
-                                         update_values | update_gradients);
+  NonMatching::MappingInfo<dim, dim> mapping_info(mapping,
+                                                  update_values |
+                                                    update_gradients);
+
+  FEPointEvaluation<dim, dim> evaluator(mapping_info, fe);
+  FEPointEvaluation<dim, dim> evaluator2(mapping_info, fe);
 
   VectorTools::interpolate(mapping, dof_handler, MyFunction<dim>(), vector);
 
@@ -108,11 +108,11 @@ test(const unsigned int degree)
                            solution_values.begin(),
                            solution_values.end());
 
-      evaluator.reinit(cell, unit_points);
+      mapping_info.reinit(cell, unit_points);
+
       evaluator.evaluate(solution_values,
                          EvaluationFlags::values | EvaluationFlags::gradients);
 
-      evaluator2.reinit(cell, unit_points, evaluator.get_mapping_data());
       evaluator2.evaluate(solution_values,
                           EvaluationFlags::values | EvaluationFlags::gradients);
 
diff --git a/tests/matrix_free/point_evaluation_17.cc b/tests/matrix_free/point_evaluation_17.cc
new file mode 100644 (file)
index 0000000..1a2c5ad
--- /dev/null
@@ -0,0 +1,163 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Check precomputed mapping information against mapping
+// information on the fly for a vector valued FE and MappingCartesian.
+// This is a modified version of point_evaluation_13.cc.
+
+#include <deal.II/base/function_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_cartesian.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+class MyFunction : public Function<dim>
+{
+public:
+  MyFunction()
+    : Function<dim>(dim)
+  {}
+
+  double
+  value(const Point<dim> &p, const unsigned int component) const override
+  {
+    return component + p[component];
+  }
+};
+
+
+
+template <int dim>
+void
+test()
+{
+  using namespace dealii;
+  Triangulation<dim> tria;
+
+  GridGenerator::subdivided_hyper_cube(tria, 2, 0, 1);
+
+  MappingCartesian<dim> mapping;
+  deallog << "Cartesian linear mapping" << std::endl;
+
+  std::vector<Point<dim>> unit_points;
+  for (unsigned int i = 0; i < 13; ++i)
+    {
+      Point<dim> p;
+      for (unsigned int d = 0; d < dim; ++d)
+        p[d] = static_cast<double>(i) / 17. + 0.015625 * d;
+      unit_points.push_back(p);
+    }
+
+  FESystem<dim>   fe(FE_Q<dim>(2), dim);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+  Vector<double> vector(dof_handler.n_dofs());
+
+  // mapping info object for precomputed mapping
+  NonMatching::MappingInfo<dim, dim> mapping_info(mapping,
+                                                  update_values |
+                                                    update_gradients);
+  // The first evaluator will use the precomputed mapping
+  FEPointEvaluation<dim, dim> evaluator1(mapping_info, fe);
+
+  // The second evaluator will compute the mapping on the fly
+  FEPointEvaluation<dim, dim> evaluator2(mapping,
+                                         fe,
+                                         update_values | update_gradients);
+
+  VectorTools::interpolate(mapping, dof_handler, MyFunction<dim>(), vector);
+
+  std::vector<double> solution_values(fe.dofs_per_cell);
+  std::vector<double> solution_values2(fe.dofs_per_cell);
+
+  FEValuesExtractors::Vector extractor(0);
+
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      cell->get_dof_values(vector,
+                           solution_values.begin(),
+                           solution_values.end());
+
+      mapping_info.reinit(cell, unit_points);
+      evaluator1.evaluate(solution_values,
+                          EvaluationFlags::values | EvaluationFlags::gradients);
+
+      evaluator2.reinit(cell, unit_points);
+      evaluator2.evaluate(solution_values,
+                          EvaluationFlags::values | EvaluationFlags::gradients);
+
+      deallog << "Cell with center " << cell->center(true) << std::endl;
+      for (unsigned int i = 0; i < unit_points.size(); ++i)
+        deallog
+          << mapping.transform_unit_to_real_cell(cell, unit_points[i]) << ": "
+          << evaluator1.get_value(i) << " error value "
+          << (evaluator1.get_value(i) - evaluator2.get_value(i)).norm()
+          << " error grad "
+          << (evaluator1.get_gradient(i) - evaluator2.get_gradient(i)).norm()
+          << std::endl;
+      deallog << std::endl;
+
+      for (unsigned int i = 0; i < unit_points.size(); ++i)
+        {
+          evaluator1.submit_value(evaluator1.get_value(i), i);
+          evaluator1.submit_gradient(evaluator1.get_gradient(i), i);
+
+          evaluator2.submit_value(evaluator2.get_value(i), i);
+          evaluator2.submit_gradient(evaluator2.get_gradient(i), i);
+        }
+
+      evaluator1.integrate(solution_values,
+                           EvaluationFlags::values |
+                             EvaluationFlags::gradients);
+      evaluator2.integrate(solution_values2,
+                           EvaluationFlags::values |
+                             EvaluationFlags::gradients);
+
+      for (unsigned int i = 0; i < solution_values.size(); ++i)
+        deallog << solution_values[i] - solution_values2[i] << ' ';
+      deallog << std::endl;
+    }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(10);
+
+  test<2>();
+  test<3>();
+}
diff --git a/tests/matrix_free/point_evaluation_17.output b/tests/matrix_free/point_evaluation_17.output
new file mode 100644 (file)
index 0000000..0a2facb
--- /dev/null
@@ -0,0 +1,195 @@
+
+DEAL::Cartesian linear mapping
+DEAL::Cell with center 0.2500000000 0.2500000000
+DEAL::0.000000000 0.007812500000: 0.000000000 1.007812500 error value 0.000000000 error grad 0.000000000
+DEAL::0.02941176471 0.03722426471: 0.02941176471 1.037224265 error value 0.000000000 error grad 0.000000000
+DEAL::0.05882352941 0.06663602941: 0.05882352941 1.066636029 error value 0.000000000 error grad 0.000000000
+DEAL::0.08823529412 0.09604779412: 0.08823529412 1.096047794 error value 0.000000000 error grad 0.000000000
+DEAL::0.1176470588 0.1254595588: 0.1176470588 1.125459559 error value 0.000000000 error grad 0.000000000
+DEAL::0.1470588235 0.1548713235: 0.1470588235 1.154871324 error value 0.000000000 error grad 0.000000000
+DEAL::0.1764705882 0.1842830882: 0.1764705882 1.184283088 error value 0.000000000 error grad 0.000000000
+DEAL::0.2058823529 0.2136948529: 0.2058823529 1.213694853 error value 0.000000000 error grad 0.000000000
+DEAL::0.2352941176 0.2431066176: 0.2352941176 1.243106618 error value 0.000000000 error grad 0.000000000
+DEAL::0.2647058824 0.2725183824: 0.2647058824 1.272518382 error value 0.000000000 error grad 0.000000000
+DEAL::0.2941176471 0.3019301471: 0.2941176471 1.301930147 error value 0.000000000 error grad 0.000000000
+DEAL::0.3235294118 0.3313419118: 0.3235294118 1.331341912 error value 0.000000000 error grad 0.000000000
+DEAL::0.3529411765 0.3607536765: 0.3529411765 1.360753676 error value 0.000000000 error grad 0.000000000
+DEAL::
+DEAL::0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::Cell with center 0.7500000000 0.2500000000
+DEAL::0.5000000000 0.007812500000: 0.5000000000 1.007812500 error value 0.000000000 error grad 0.000000000
+DEAL::0.5294117647 0.03722426471: 0.5294117647 1.037224265 error value 0.000000000 error grad 0.000000000
+DEAL::0.5588235294 0.06663602941: 0.5588235294 1.066636029 error value 0.000000000 error grad 0.000000000
+DEAL::0.5882352941 0.09604779412: 0.5882352941 1.096047794 error value 0.000000000 error grad 0.000000000
+DEAL::0.6176470588 0.1254595588: 0.6176470588 1.125459559 error value 0.000000000 error grad 0.000000000
+DEAL::0.6470588235 0.1548713235: 0.6470588235 1.154871324 error value 0.000000000 error grad 0.000000000
+DEAL::0.6764705882 0.1842830882: 0.6764705882 1.184283088 error value 0.000000000 error grad 0.000000000
+DEAL::0.7058823529 0.2136948529: 0.7058823529 1.213694853 error value 0.000000000 error grad 0.000000000
+DEAL::0.7352941176 0.2431066176: 0.7352941176 1.243106618 error value 0.000000000 error grad 0.000000000
+DEAL::0.7647058824 0.2725183824: 0.7647058824 1.272518382 error value 0.000000000 error grad 0.000000000
+DEAL::0.7941176471 0.3019301471: 0.7941176471 1.301930147 error value 0.000000000 error grad 0.000000000
+DEAL::0.8235294118 0.3313419118: 0.8235294118 1.331341912 error value 0.000000000 error grad 0.000000000
+DEAL::0.8529411765 0.3607536765: 0.8529411765 1.360753676 error value 0.000000000 error grad 0.000000000
+DEAL::
+DEAL::0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::Cell with center 0.2500000000 0.7500000000
+DEAL::0.000000000 0.5078125000: 0.000000000 1.507812500 error value 0.000000000 error grad 0.000000000
+DEAL::0.02941176471 0.5372242647: 0.02941176471 1.537224265 error value 0.000000000 error grad 0.000000000
+DEAL::0.05882352941 0.5666360294: 0.05882352941 1.566636029 error value 0.000000000 error grad 0.000000000
+DEAL::0.08823529412 0.5960477941: 0.08823529412 1.596047794 error value 0.000000000 error grad 0.000000000
+DEAL::0.1176470588 0.6254595588: 0.1176470588 1.625459559 error value 0.000000000 error grad 0.000000000
+DEAL::0.1470588235 0.6548713235: 0.1470588235 1.654871324 error value 0.000000000 error grad 0.000000000
+DEAL::0.1764705882 0.6842830882: 0.1764705882 1.684283088 error value 0.000000000 error grad 0.000000000
+DEAL::0.2058823529 0.7136948529: 0.2058823529 1.713694853 error value 0.000000000 error grad 0.000000000
+DEAL::0.2352941176 0.7431066176: 0.2352941176 1.743106618 error value 0.000000000 error grad 0.000000000
+DEAL::0.2647058824 0.7725183824: 0.2647058824 1.772518382 error value 0.000000000 error grad 0.000000000
+DEAL::0.2941176471 0.8019301471: 0.2941176471 1.801930147 error value 0.000000000 error grad 0.000000000
+DEAL::0.3235294118 0.8313419118: 0.3235294118 1.831341912 error value 0.000000000 error grad 0.000000000
+DEAL::0.3529411765 0.8607536765: 0.3529411765 1.860753676 error value 0.000000000 error grad 0.000000000
+DEAL::
+DEAL::0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::Cell with center 0.7500000000 0.7500000000
+DEAL::0.5000000000 0.5078125000: 0.5000000000 1.507812500 error value 0.000000000 error grad 0.000000000
+DEAL::0.5294117647 0.5372242647: 0.5294117647 1.537224265 error value 0.000000000 error grad 0.000000000
+DEAL::0.5588235294 0.5666360294: 0.5588235294 1.566636029 error value 0.000000000 error grad 0.000000000
+DEAL::0.5882352941 0.5960477941: 0.5882352941 1.596047794 error value 0.000000000 error grad 0.000000000
+DEAL::0.6176470588 0.6254595588: 0.6176470588 1.625459559 error value 0.000000000 error grad 0.000000000
+DEAL::0.6470588235 0.6548713235: 0.6470588235 1.654871324 error value 0.000000000 error grad 0.000000000
+DEAL::0.6764705882 0.6842830882: 0.6764705882 1.684283088 error value 0.000000000 error grad 0.000000000
+DEAL::0.7058823529 0.7136948529: 0.7058823529 1.713694853 error value 0.000000000 error grad 0.000000000
+DEAL::0.7352941176 0.7431066176: 0.7352941176 1.743106618 error value 0.000000000 error grad 0.000000000
+DEAL::0.7647058824 0.7725183824: 0.7647058824 1.772518382 error value 0.000000000 error grad 0.000000000
+DEAL::0.7941176471 0.8019301471: 0.7941176471 1.801930147 error value 0.000000000 error grad 0.000000000
+DEAL::0.8235294118 0.8313419118: 0.8235294118 1.831341912 error value 0.000000000 error grad 0.000000000
+DEAL::0.8529411765 0.8607536765: 0.8529411765 1.860753676 error value 0.000000000 error grad 0.000000000
+DEAL::
+DEAL::0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::Cartesian linear mapping
+DEAL::Cell with center 0.2500000000 0.2500000000 0.2500000000
+DEAL::0.000000000 0.007812500000 0.01562500000: 0.000000000 1.007812500 2.015625000 error value 0.000000000 error grad 0.000000000
+DEAL::0.02941176471 0.03722426471 0.04503676471: 0.02941176471 1.037224265 2.045036765 error value 0.000000000 error grad 0.000000000
+DEAL::0.05882352941 0.06663602941 0.07444852941: 0.05882352941 1.066636029 2.074448529 error value 0.000000000 error grad 0.000000000
+DEAL::0.08823529412 0.09604779412 0.1038602941: 0.08823529412 1.096047794 2.103860294 error value 0.000000000 error grad 0.000000000
+DEAL::0.1176470588 0.1254595588 0.1332720588: 0.1176470588 1.125459559 2.133272059 error value 0.000000000 error grad 0.000000000
+DEAL::0.1470588235 0.1548713235 0.1626838235: 0.1470588235 1.154871324 2.162683824 error value 0.000000000 error grad 0.000000000
+DEAL::0.1764705882 0.1842830882 0.1920955882: 0.1764705882 1.184283088 2.192095588 error value 0.000000000 error grad 0.000000000
+DEAL::0.2058823529 0.2136948529 0.2215073529: 0.2058823529 1.213694853 2.221507353 error value 0.000000000 error grad 0.000000000
+DEAL::0.2352941176 0.2431066176 0.2509191176: 0.2352941176 1.243106618 2.250919118 error value 0.000000000 error grad 0.000000000
+DEAL::0.2647058824 0.2725183824 0.2803308824: 0.2647058824 1.272518382 2.280330882 error value 0.000000000 error grad 0.000000000
+DEAL::0.2941176471 0.3019301471 0.3097426471: 0.2941176471 1.301930147 2.309742647 error value 0.000000000 error grad 0.000000000
+DEAL::0.3235294118 0.3313419118 0.3391544118: 0.3235294118 1.331341912 2.339154412 error value 0.000000000 error grad 0.000000000
+DEAL::0.3529411765 0.3607536765 0.3685661765: 0.3529411765 1.360753676 2.368566176 error value 0.000000000 error grad 0.000000000
+DEAL::
+DEAL::0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::Cell with center 0.7500000000 0.2500000000 0.2500000000
+DEAL::0.5000000000 0.007812500000 0.01562500000: 0.5000000000 1.007812500 2.015625000 error value 0.000000000 error grad 0.000000000
+DEAL::0.5294117647 0.03722426471 0.04503676471: 0.5294117647 1.037224265 2.045036765 error value 0.000000000 error grad 0.000000000
+DEAL::0.5588235294 0.06663602941 0.07444852941: 0.5588235294 1.066636029 2.074448529 error value 0.000000000 error grad 0.000000000
+DEAL::0.5882352941 0.09604779412 0.1038602941: 0.5882352941 1.096047794 2.103860294 error value 0.000000000 error grad 0.000000000
+DEAL::0.6176470588 0.1254595588 0.1332720588: 0.6176470588 1.125459559 2.133272059 error value 0.000000000 error grad 0.000000000
+DEAL::0.6470588235 0.1548713235 0.1626838235: 0.6470588235 1.154871324 2.162683824 error value 0.000000000 error grad 0.000000000
+DEAL::0.6764705882 0.1842830882 0.1920955882: 0.6764705882 1.184283088 2.192095588 error value 0.000000000 error grad 0.000000000
+DEAL::0.7058823529 0.2136948529 0.2215073529: 0.7058823529 1.213694853 2.221507353 error value 0.000000000 error grad 0.000000000
+DEAL::0.7352941176 0.2431066176 0.2509191176: 0.7352941176 1.243106618 2.250919118 error value 0.000000000 error grad 0.000000000
+DEAL::0.7647058824 0.2725183824 0.2803308824: 0.7647058824 1.272518382 2.280330882 error value 0.000000000 error grad 0.000000000
+DEAL::0.7941176471 0.3019301471 0.3097426471: 0.7941176471 1.301930147 2.309742647 error value 0.000000000 error grad 0.000000000
+DEAL::0.8235294118 0.3313419118 0.3391544118: 0.8235294118 1.331341912 2.339154412 error value 0.000000000 error grad 0.000000000
+DEAL::0.8529411765 0.3607536765 0.3685661765: 0.8529411765 1.360753676 2.368566176 error value 0.000000000 error grad 0.000000000
+DEAL::
+DEAL::0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::Cell with center 0.2500000000 0.7500000000 0.2500000000
+DEAL::0.000000000 0.5078125000 0.01562500000: 0.000000000 1.507812500 2.015625000 error value 0.000000000 error grad 0.000000000
+DEAL::0.02941176471 0.5372242647 0.04503676471: 0.02941176471 1.537224265 2.045036765 error value 0.000000000 error grad 0.000000000
+DEAL::0.05882352941 0.5666360294 0.07444852941: 0.05882352941 1.566636029 2.074448529 error value 0.000000000 error grad 0.000000000
+DEAL::0.08823529412 0.5960477941 0.1038602941: 0.08823529412 1.596047794 2.103860294 error value 0.000000000 error grad 0.000000000
+DEAL::0.1176470588 0.6254595588 0.1332720588: 0.1176470588 1.625459559 2.133272059 error value 0.000000000 error grad 0.000000000
+DEAL::0.1470588235 0.6548713235 0.1626838235: 0.1470588235 1.654871324 2.162683824 error value 0.000000000 error grad 0.000000000
+DEAL::0.1764705882 0.6842830882 0.1920955882: 0.1764705882 1.684283088 2.192095588 error value 0.000000000 error grad 0.000000000
+DEAL::0.2058823529 0.7136948529 0.2215073529: 0.2058823529 1.713694853 2.221507353 error value 0.000000000 error grad 0.000000000
+DEAL::0.2352941176 0.7431066176 0.2509191176: 0.2352941176 1.743106618 2.250919118 error value 0.000000000 error grad 0.000000000
+DEAL::0.2647058824 0.7725183824 0.2803308824: 0.2647058824 1.772518382 2.280330882 error value 0.000000000 error grad 0.000000000
+DEAL::0.2941176471 0.8019301471 0.3097426471: 0.2941176471 1.801930147 2.309742647 error value 0.000000000 error grad 0.000000000
+DEAL::0.3235294118 0.8313419118 0.3391544118: 0.3235294118 1.831341912 2.339154412 error value 0.000000000 error grad 0.000000000
+DEAL::0.3529411765 0.8607536765 0.3685661765: 0.3529411765 1.860753676 2.368566176 error value 0.000000000 error grad 0.000000000
+DEAL::
+DEAL::0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::Cell with center 0.7500000000 0.7500000000 0.2500000000
+DEAL::0.5000000000 0.5078125000 0.01562500000: 0.5000000000 1.507812500 2.015625000 error value 0.000000000 error grad 0.000000000
+DEAL::0.5294117647 0.5372242647 0.04503676471: 0.5294117647 1.537224265 2.045036765 error value 0.000000000 error grad 0.000000000
+DEAL::0.5588235294 0.5666360294 0.07444852941: 0.5588235294 1.566636029 2.074448529 error value 0.000000000 error grad 0.000000000
+DEAL::0.5882352941 0.5960477941 0.1038602941: 0.5882352941 1.596047794 2.103860294 error value 0.000000000 error grad 0.000000000
+DEAL::0.6176470588 0.6254595588 0.1332720588: 0.6176470588 1.625459559 2.133272059 error value 0.000000000 error grad 0.000000000
+DEAL::0.6470588235 0.6548713235 0.1626838235: 0.6470588235 1.654871324 2.162683824 error value 0.000000000 error grad 0.000000000
+DEAL::0.6764705882 0.6842830882 0.1920955882: 0.6764705882 1.684283088 2.192095588 error value 0.000000000 error grad 0.000000000
+DEAL::0.7058823529 0.7136948529 0.2215073529: 0.7058823529 1.713694853 2.221507353 error value 0.000000000 error grad 0.000000000
+DEAL::0.7352941176 0.7431066176 0.2509191176: 0.7352941176 1.743106618 2.250919118 error value 0.000000000 error grad 0.000000000
+DEAL::0.7647058824 0.7725183824 0.2803308824: 0.7647058824 1.772518382 2.280330882 error value 0.000000000 error grad 0.000000000
+DEAL::0.7941176471 0.8019301471 0.3097426471: 0.7941176471 1.801930147 2.309742647 error value 0.000000000 error grad 0.000000000
+DEAL::0.8235294118 0.8313419118 0.3391544118: 0.8235294118 1.831341912 2.339154412 error value 0.000000000 error grad 0.000000000
+DEAL::0.8529411765 0.8607536765 0.3685661765: 0.8529411765 1.860753676 2.368566176 error value 0.000000000 error grad 0.000000000
+DEAL::
+DEAL::0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::Cell with center 0.2500000000 0.2500000000 0.7500000000
+DEAL::0.000000000 0.007812500000 0.5156250000: 0.000000000 1.007812500 2.515625000 error value 0.000000000 error grad 0.000000000
+DEAL::0.02941176471 0.03722426471 0.5450367647: 0.02941176471 1.037224265 2.545036765 error value 0.000000000 error grad 0.000000000
+DEAL::0.05882352941 0.06663602941 0.5744485294: 0.05882352941 1.066636029 2.574448529 error value 0.000000000 error grad 0.000000000
+DEAL::0.08823529412 0.09604779412 0.6038602941: 0.08823529412 1.096047794 2.603860294 error value 0.000000000 error grad 0.000000000
+DEAL::0.1176470588 0.1254595588 0.6332720588: 0.1176470588 1.125459559 2.633272059 error value 0.000000000 error grad 0.000000000
+DEAL::0.1470588235 0.1548713235 0.6626838235: 0.1470588235 1.154871324 2.662683824 error value 0.000000000 error grad 0.000000000
+DEAL::0.1764705882 0.1842830882 0.6920955882: 0.1764705882 1.184283088 2.692095588 error value 0.000000000 error grad 0.000000000
+DEAL::0.2058823529 0.2136948529 0.7215073529: 0.2058823529 1.213694853 2.721507353 error value 0.000000000 error grad 0.000000000
+DEAL::0.2352941176 0.2431066176 0.7509191176: 0.2352941176 1.243106618 2.750919118 error value 0.000000000 error grad 0.000000000
+DEAL::0.2647058824 0.2725183824 0.7803308824: 0.2647058824 1.272518382 2.780330882 error value 0.000000000 error grad 0.000000000
+DEAL::0.2941176471 0.3019301471 0.8097426471: 0.2941176471 1.301930147 2.809742647 error value 0.000000000 error grad 0.000000000
+DEAL::0.3235294118 0.3313419118 0.8391544118: 0.3235294118 1.331341912 2.839154412 error value 0.000000000 error grad 0.000000000
+DEAL::0.3529411765 0.3607536765 0.8685661765: 0.3529411765 1.360753676 2.868566176 error value 0.000000000 error grad 0.000000000
+DEAL::
+DEAL::0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::Cell with center 0.7500000000 0.2500000000 0.7500000000
+DEAL::0.5000000000 0.007812500000 0.5156250000: 0.5000000000 1.007812500 2.515625000 error value 0.000000000 error grad 0.000000000
+DEAL::0.5294117647 0.03722426471 0.5450367647: 0.5294117647 1.037224265 2.545036765 error value 0.000000000 error grad 0.000000000
+DEAL::0.5588235294 0.06663602941 0.5744485294: 0.5588235294 1.066636029 2.574448529 error value 0.000000000 error grad 0.000000000
+DEAL::0.5882352941 0.09604779412 0.6038602941: 0.5882352941 1.096047794 2.603860294 error value 0.000000000 error grad 0.000000000
+DEAL::0.6176470588 0.1254595588 0.6332720588: 0.6176470588 1.125459559 2.633272059 error value 0.000000000 error grad 0.000000000
+DEAL::0.6470588235 0.1548713235 0.6626838235: 0.6470588235 1.154871324 2.662683824 error value 0.000000000 error grad 0.000000000
+DEAL::0.6764705882 0.1842830882 0.6920955882: 0.6764705882 1.184283088 2.692095588 error value 0.000000000 error grad 0.000000000
+DEAL::0.7058823529 0.2136948529 0.7215073529: 0.7058823529 1.213694853 2.721507353 error value 0.000000000 error grad 0.000000000
+DEAL::0.7352941176 0.2431066176 0.7509191176: 0.7352941176 1.243106618 2.750919118 error value 0.000000000 error grad 0.000000000
+DEAL::0.7647058824 0.2725183824 0.7803308824: 0.7647058824 1.272518382 2.780330882 error value 0.000000000 error grad 0.000000000
+DEAL::0.7941176471 0.3019301471 0.8097426471: 0.7941176471 1.301930147 2.809742647 error value 0.000000000 error grad 0.000000000
+DEAL::0.8235294118 0.3313419118 0.8391544118: 0.8235294118 1.331341912 2.839154412 error value 0.000000000 error grad 0.000000000
+DEAL::0.8529411765 0.3607536765 0.8685661765: 0.8529411765 1.360753676 2.868566176 error value 0.000000000 error grad 0.000000000
+DEAL::
+DEAL::0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::Cell with center 0.2500000000 0.7500000000 0.7500000000
+DEAL::0.000000000 0.5078125000 0.5156250000: 0.000000000 1.507812500 2.515625000 error value 0.000000000 error grad 0.000000000
+DEAL::0.02941176471 0.5372242647 0.5450367647: 0.02941176471 1.537224265 2.545036765 error value 0.000000000 error grad 0.000000000
+DEAL::0.05882352941 0.5666360294 0.5744485294: 0.05882352941 1.566636029 2.574448529 error value 0.000000000 error grad 0.000000000
+DEAL::0.08823529412 0.5960477941 0.6038602941: 0.08823529412 1.596047794 2.603860294 error value 0.000000000 error grad 0.000000000
+DEAL::0.1176470588 0.6254595588 0.6332720588: 0.1176470588 1.625459559 2.633272059 error value 0.000000000 error grad 0.000000000
+DEAL::0.1470588235 0.6548713235 0.6626838235: 0.1470588235 1.654871324 2.662683824 error value 0.000000000 error grad 0.000000000
+DEAL::0.1764705882 0.6842830882 0.6920955882: 0.1764705882 1.684283088 2.692095588 error value 0.000000000 error grad 0.000000000
+DEAL::0.2058823529 0.7136948529 0.7215073529: 0.2058823529 1.713694853 2.721507353 error value 0.000000000 error grad 0.000000000
+DEAL::0.2352941176 0.7431066176 0.7509191176: 0.2352941176 1.743106618 2.750919118 error value 0.000000000 error grad 0.000000000
+DEAL::0.2647058824 0.7725183824 0.7803308824: 0.2647058824 1.772518382 2.780330882 error value 0.000000000 error grad 0.000000000
+DEAL::0.2941176471 0.8019301471 0.8097426471: 0.2941176471 1.801930147 2.809742647 error value 0.000000000 error grad 0.000000000
+DEAL::0.3235294118 0.8313419118 0.8391544118: 0.3235294118 1.831341912 2.839154412 error value 0.000000000 error grad 0.000000000
+DEAL::0.3529411765 0.8607536765 0.8685661765: 0.3529411765 1.860753676 2.868566176 error value 0.000000000 error grad 0.000000000
+DEAL::
+DEAL::0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::Cell with center 0.7500000000 0.7500000000 0.7500000000
+DEAL::0.5000000000 0.5078125000 0.5156250000: 0.5000000000 1.507812500 2.515625000 error value 0.000000000 error grad 0.000000000
+DEAL::0.5294117647 0.5372242647 0.5450367647: 0.5294117647 1.537224265 2.545036765 error value 0.000000000 error grad 0.000000000
+DEAL::0.5588235294 0.5666360294 0.5744485294: 0.5588235294 1.566636029 2.574448529 error value 0.000000000 error grad 0.000000000
+DEAL::0.5882352941 0.5960477941 0.6038602941: 0.5882352941 1.596047794 2.603860294 error value 0.000000000 error grad 0.000000000
+DEAL::0.6176470588 0.6254595588 0.6332720588: 0.6176470588 1.625459559 2.633272059 error value 0.000000000 error grad 0.000000000
+DEAL::0.6470588235 0.6548713235 0.6626838235: 0.6470588235 1.654871324 2.662683824 error value 0.000000000 error grad 0.000000000
+DEAL::0.6764705882 0.6842830882 0.6920955882: 0.6764705882 1.684283088 2.692095588 error value 0.000000000 error grad 0.000000000
+DEAL::0.7058823529 0.7136948529 0.7215073529: 0.7058823529 1.713694853 2.721507353 error value 0.000000000 error grad 0.000000000
+DEAL::0.7352941176 0.7431066176 0.7509191176: 0.7352941176 1.743106618 2.750919118 error value 0.000000000 error grad 0.000000000
+DEAL::0.7647058824 0.7725183824 0.7803308824: 0.7647058824 1.772518382 2.780330882 error value 0.000000000 error grad 0.000000000
+DEAL::0.7941176471 0.8019301471 0.8097426471: 0.7941176471 1.801930147 2.809742647 error value 0.000000000 error grad 0.000000000
+DEAL::0.8235294118 0.8313419118 0.8391544118: 0.8235294118 1.831341912 2.839154412 error value 0.000000000 error grad 0.000000000
+DEAL::0.8529411765 0.8607536765 0.8685661765: 0.8529411765 1.860753676 2.868566176 error value 0.000000000 error grad 0.000000000
+DEAL::
+DEAL::0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 

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.