]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make partial initialization work properly, including a minimum functionality in case...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 30 Apr 2013 14:56:37 +0000 (14:56 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 30 Apr 2013 14:56:37 +0000 (14:56 +0000)
git-svn-id: https://svn.dealii.org/trunk@29408 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/matrix_free/dof_info.h
deal.II/include/deal.II/matrix_free/dof_info.templates.h
deal.II/include/deal.II/matrix_free/fe_evaluation.h
deal.II/include/deal.II/matrix_free/matrix_free.templates.h
tests/matrix_free/get_functions_common.h
tests/matrix_free/no_index_initialize.cc [new file with mode: 0644]
tests/matrix_free/no_index_initialize/cmp/generic [new file with mode: 0644]
tests/matrix_free/update_mapping_only.cc [new file with mode: 0644]
tests/matrix_free/update_mapping_only/cmp/generic [new file with mode: 0644]

index 90c7fccdabb9a907fc3f8f56580f7891f3d4344f..b7a791dd2bc4b75e8fd2249d35f5f5f5c676571a 100644 (file)
@@ -354,6 +354,13 @@ namespace internal
        */
       std::vector<unsigned int> plain_dof_indices;
 
+      /**
+       * Stores the dimension of the underlying DoFHandler. Since the indices
+       * are not templated, this is the variable that makes the dimension
+       * accessible in the (rare) cases it is needed inside this class.
+       */
+      unsigned int dimension;
+
       /**
        * Stores the number of components in the DoFHandler where the indices
        * have been read from.
index 2430dcecdb386efb36376fbc97a1d57b175538f2..64623405d56d36f5231f3866a8d3658f5c33c070 100644 (file)
@@ -131,6 +131,7 @@ namespace internal
       constrained_dofs (dof_info_in.constrained_dofs),
       row_starts_plain_indices (dof_info_in.row_starts_plain_indices),
       plain_dof_indices (dof_info_in.plain_dof_indices),
+      dimension (dof_info_in.dimension),
       n_components (dof_info_in.n_components),
       dofs_per_cell (dof_info_in.dofs_per_cell),
       dofs_per_face (dof_info_in.dofs_per_face),
@@ -153,6 +154,7 @@ namespace internal
       ghost_dofs.clear();
       dofs_per_cell.clear();
       dofs_per_face.clear();
+      dimension = 2;
       n_components = 0;
       row_starts_plain_indices.clear();
       plain_dof_indices.clear();
index 9acb29a336fab92456f3ea4384bf6cd12b25ec1e..187716af64e9a9e3fb2ab7bac8f7453071014419 100644 (file)
@@ -687,22 +687,6 @@ protected:
    */
   unsigned int cell_data_number;
 
-  /**
-   * If the present cell chunk for vectorization is not completely filled up
-   * with data, this field stores how many physical cells are underlying. Is
-   * between 1 and VectorizedArray<Number>::n_array_elements-1 (inclusive).
-   */
-  unsigned int n_irreg_components_filled;
-
-  /**
-   * Stores whether the present cell chunk used in vectorization is not
-   * completely filled up with physical cells. E.g. if vectorization dictates
-   * that four cells should be worked with but only three physical cells are
-   * left, this flag will be set to true, otherwise to false. Mainly used for
-   * internal checking when reading from vectors or writing to vectors.
-   */
-  bool at_irregular_cell;
-
   /**
    * Debug information to track whether dof values have been initialized
    * before accessed. Used to control exceptions when uninitialized data is
@@ -1465,17 +1449,14 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
   jacobian_grad_upper(0),
   cell               (numbers::invalid_unsigned_int),
   cell_type          (internal::MatrixFreeFunctions::undefined),
-  cell_data_number   (0),
-  n_irreg_components_filled (0),
-  at_irregular_cell  (false)
+  cell_data_number   (0)
 {
-  Assert (matrix_info.indices_initialized() == true,
-          ExcNotInitialized());
   Assert (matrix_info.mapping_initialized() == true,
           ExcNotInitialized());
   AssertDimension (matrix_info.get_size_info().vectorization_length,
                    VectorizedArray<Number>::n_array_elements);
   Assert (n_fe_components == 1 ||
+          n_components == 1 ||
           n_components == n_fe_components,
           ExcMessage ("The underlying FE is vector-valued. In this case, the "
                       "template argument n_components must be a the same "
@@ -1553,9 +1534,6 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
         }
     }
 
-  n_irreg_components_filled =
-    std_cxx1x::get<2>(dof_info.row_starts[cell_in]);
-  at_irregular_cell = n_irreg_components_filled > 0;
 #ifdef DEBUG
   dof_values_initialized      = false;
   values_quad_initialized     = false;
@@ -1832,6 +1810,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
   // into the local data field or write local data into the vector. Certain
   // operations are no-ops for the given use case.
 
+  Assert (matrix_info.indices_initialized() == true,
+          ExcNotInitialized());
   Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
 
   // loop over all local dofs. ind_local holds local number on cell, index
@@ -1844,6 +1824,10 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
     dof_info.end_indicators(cell);
   unsigned int ind_local = 0;
 
+  const unsigned int n_irreg_components_filled =
+    std_cxx1x::get<2>(dof_info.row_starts[cell]);
+  const bool at_irregular_cell = n_irreg_components_filled > 0;
+
   // scalar case (or case when all components have the same degrees of freedom
   // and sit on a different vector each)
   if (n_fe_components == 1)
@@ -2447,6 +2431,8 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
 {
   // this is different from the other three operations because we do not use
   // constraints here, so this is a separate function.
+  Assert (matrix_info.indices_initialized() == true,
+          ExcNotInitialized());
   Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
   Assert (dof_info.store_plain_indices == true, ExcNotInitialized());
 
@@ -2455,6 +2441,10 @@ FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
   // points to the global indices stored in index_local_to_global
   const unsigned int *dof_indices = dof_info.begin_indices_plain(cell);
 
+  const unsigned int n_irreg_components_filled =
+    std_cxx1x::get<2>(dof_info.row_starts[cell]);
+  const bool at_irregular_cell = n_irreg_components_filled > 0;
+
   // scalar case (or case when all components have the same degrees of freedom
   // and sit on a different vector each)
   if (n_fe_components == 1)
index ccb397894f85200ba25d66535e006386ffb49ff9..2014ab843504df9a9ab51d2db854aeeb16b13423 100644 (file)
@@ -159,6 +159,32 @@ internal_reinit(const Mapping<dim>                          &mapping,
       initialize_indices (constraint, locally_owned_set);
     }
 
+  // initialize bare structures
+  else if (dof_info.size() != dof_handler.size())
+    {
+      initialize_dof_handlers(dof_handler, additional_data.level_mg_handler);
+      std::vector<unsigned int> dummy;
+      size_info.make_layout (cell_level_index.size(),
+                             VectorizedArray<Number>::n_array_elements,
+                             dummy, dummy);
+      for (unsigned int i=0; i<dof_info.size(); ++i)
+        {
+          dof_info[i].dimension    = dim;
+          dof_info[i].n_components = dof_handler[i]->get_fe().element_multiplicity(0);
+          dof_info[i].dofs_per_cell.push_back(dof_handler[i]->get_fe().dofs_per_cell);
+          dof_info[i].row_starts.resize(size_info.n_macro_cells+1);
+          std_cxx1x::get<2>(dof_info[i].row_starts.back()) =
+            cell_level_index.size() % VectorizedArray<Number>::n_array_elements;
+
+          // if indices are not initialized, the cell_level_index might not be
+          // divisible by the vectorization length. But it must be for
+          // mapping_info...
+          while (cell_level_index.size() % VectorizedArray<Number>::n_array_elements
+                 != 0)
+            cell_level_index.push_back(cell_level_index.back());
+        }
+    }
+
   // Reads out the FE information and stores the shape function values,
   // gradients and Hessians for quadrature points.
   const unsigned int n_fe   = dof_handler.size();
@@ -254,6 +280,33 @@ internal_reinit(const Mapping<dim>                            &mapping,
       initialize_indices (constraint, locally_owned_set);
     }
 
+  // initialize bare structures
+  else if (dof_info.size() != dof_handler.size())
+    {
+      initialize_dof_handlers(dof_handler, additional_data.level_mg_handler);
+      std::vector<unsigned int> dummy;
+      size_info.make_layout (cell_level_index.size(),
+                             VectorizedArray<Number>::n_array_elements,
+                             dummy, dummy);
+      for (unsigned int i=0; i<dof_info.size(); ++i)
+        {
+          Assert(dof_handler[i]->get_fe().size() == 1, ExcNotImplemented());
+          dof_info[i].dimension    = dim;
+          dof_info[i].n_components = dof_handler[i]->get_fe()[0].element_multiplicity(0);
+          dof_info[i].dofs_per_cell.push_back(dof_handler[i]->get_fe()[0].dofs_per_cell);
+          dof_info[i].row_starts.resize(size_info.n_macro_cells+1);
+          std_cxx1x::get<2>(dof_info[i].row_starts.back()) =
+            cell_level_index.size() % VectorizedArray<Number>::n_array_elements;
+
+          // if indices are not initialized, the cell_level_index might not be
+          // divisible by the vectorization length. But it must be for
+          // mapping_info...
+          while (cell_level_index.size() % VectorizedArray<Number>::n_array_elements
+                 != 0)
+            cell_level_index.push_back(cell_level_index.back());
+        }
+    }
+
   // Reads out the FE information and stores the shape function values,
   // gradients and Hessians for quadrature points.
   const unsigned int n_components = dof_handler.size();
@@ -460,8 +513,8 @@ void MatrixFree<dim,Number>::initialize_indices
           // cache number of finite elements and dofs_per_cell
           dof_info[no].dofs_per_cell.push_back (fe.dofs_per_cell);
           dof_info[no].dofs_per_face.push_back (fe.dofs_per_face);
-          dof_info[no].n_components  = n_fe_components;
-
+          dof_info[no].dimension    = dim;
+          dof_info[no].n_components = n_fe_components;
 
           // get permutation that gives lexicographic renumbering of the cell
           // dofs renumber (this is necessary for FE_Q, for example, since
index b859b07386ba30edb35fc376aee97aaa7bb7415e..f3e66e8c2522094680258bf6a743bd0e2e9f044a 100644 (file)
@@ -59,7 +59,7 @@ class MatrixFreeTest
 
     virtual ~MatrixFreeTest ()
       {}
-    
+
                                 // make function virtual to allow derived
                                 // classes to define a different function
   virtual void
@@ -269,4 +269,3 @@ int main ()
     deallog.pop();
   }
 }
-
diff --git a/tests/matrix_free/no_index_initialize.cc b/tests/matrix_free/no_index_initialize.cc
new file mode 100644 (file)
index 0000000..f3736b8
--- /dev/null
@@ -0,0 +1,183 @@
+//------------------  no_index_initialize.cc  ------------------------
+//    $Id$
+//    Version: $Name$
+//
+//------------------  no_index_initialize.cc  ------------------------
+
+
+// check that FEEvaluation can be evaluated without indices initialized (and
+// throws an exception when trying to read/write from/to vectors)
+
+#include "../tests.h"
+
+
+std::ofstream logfile("no_index_initialize/output");
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+
+// forward declare this function
+template <int dim, int fe_degree>
+void test ();
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d=fe_degree+1, typename Number=double>
+class MatrixFreeTest
+{
+ public:
+  bool read_vector;
+
+  MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+    read_vector(false),
+    data       (data_in)
+  {};
+
+  void
+  operator () (const MatrixFree<dim,Number> &data,
+               Vector<Number> &,
+               const Vector<Number> &src,
+               const std::pair<unsigned int,unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim,fe_degree,n_q_points_1d,1,Number> fe_eval (data);
+
+    for(unsigned int cell=cell_range.first;cell<cell_range.second;++cell)
+      {
+        fe_eval.reinit (cell);
+
+        // do not read from vector as that is disallowed
+        if (read_vector)
+          fe_eval.read_dof_values(src);
+        else
+          for (unsigned int e=0; e<fe_eval.dofs_per_cell; ++e)
+            fe_eval.submit_dof_value(make_vectorized_array<Number>(1.), e);
+        fe_eval.evaluate (true,true,true);
+
+        // values should evaluate to one, derivatives to zero
+        for (unsigned int q=0; q<fe_eval.n_q_points; ++q)
+          for (unsigned int d=0; d<VectorizedArray<Number>::n_array_elements; ++d)
+            {
+              Assert(fe_eval.get_value(q)[d] == 1., ExcInternalError());
+              for (unsigned int e=0; e<dim; ++e)
+                Assert(fe_eval.get_gradient(q)[e][d] == 0., ExcInternalError());
+              Assert(fe_eval.get_laplacian(q)[d] == 0., ExcInternalError());
+            }
+      }
+  }
+
+
+
+  void test_functions (const Vector<Number> &src) const
+  {
+    Vector<Number> dst_dummy;
+    data.cell_loop (&MatrixFreeTest::operator(), this, dst_dummy, src);
+    deallog << "OK" << std::endl;
+  };
+
+protected:
+  const MatrixFree<dim,Number> &data;
+};
+
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof,
+              const ConstraintMatrix&constraints)
+{
+  // use this for info on problem
+  //std::cout << "Number of cells: " << dof.get_tria().n_active_cells()
+  //          << std::endl;
+  //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
+  //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl;
+
+  Vector<number> solution (dof.n_dofs());
+
+                                // create vector with random entries
+  for (unsigned int i=0; i<dof.n_dofs(); ++i)
+    {
+      if(constraints.is_constrained(i))
+        continue;
+      const double entry = rand()/(double)RAND_MAX;
+      solution(i) = entry;
+    }
+
+  constraints.distribute(solution);
+  MatrixFree<dim,number> mf_data;
+  {
+    const QGauss<1> quad (fe_degree+1);
+    typename MatrixFree<dim,number>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim,number>::AdditionalData::none;
+    data.mapping_update_flags = update_gradients | update_second_derivatives;
+    data.initialize_indices = false;
+    mf_data.reinit (dof, constraints, quad, data);
+  }
+
+  deallog << "Testing " << dof.get_fe().get_name() << " without read" << std::endl;
+  MatrixFreeTest<dim,fe_degree,fe_degree+1,number> mf (mf_data);
+  mf.test_functions(solution);
+
+  deallog << "Testing " << dof.get_fe().get_name() << " with read" << std::endl;
+  try
+    {
+      mf.read_vector = true;
+      mf.test_functions(solution);
+    }
+  catch (ExceptionBase &e)
+    {
+      deallog << e.get_exc_name() << std::endl;
+    }
+}
+
+
+int main ()
+{
+  deal_II_exceptions::disable_abort_on_exception();
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  deallog << std::setprecision (3);
+  test<2,1>();
+}
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_ball (tria);
+  static const HyperBallBoundary<dim> boundary;
+  tria.set_boundary (0, boundary);
+                                // refine first and last cell
+  tria.begin(tria.n_levels()-1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  tria.refine_global (4-dim);
+
+  FE_Q<dim> fe (fe_degree);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  constraints.close();
+
+  do_test <dim, fe_degree, double> (dof, constraints);
+}
diff --git a/tests/matrix_free/no_index_initialize/cmp/generic b/tests/matrix_free/no_index_initialize/cmp/generic
new file mode 100644 (file)
index 0000000..1a9723e
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::Testing FE_Q<2>(1) without read
+DEAL::OK
+DEAL::Testing FE_Q<2>(1) with read
+DEAL::ExcNotInitialized()
diff --git a/tests/matrix_free/update_mapping_only.cc b/tests/matrix_free/update_mapping_only.cc
new file mode 100644 (file)
index 0000000..d0a619f
--- /dev/null
@@ -0,0 +1,84 @@
+//------------------  update_mapping_only.cc  ------------------------
+//    $Id$
+//    Version: $Name$
+//
+//------------------  update_mapping_only.cc  ------------------------
+
+
+// correctness of matrix-free initialization when the mapping changes and the
+// indices are not updated, only the mapping info
+
+#include "../tests.h"
+
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q_eulerian.h>
+
+
+std::ofstream logfile("update_mapping_only/output");
+
+#include "get_functions_common.h"
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube (tria);
+  tria.refine_global(1);
+
+  FE_Q<dim> fe (fe_degree);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+
+  ConstraintMatrix constraints;
+  constraints.close();
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+  // use this for info on problem
+  //std::cout << "Number of cells: " << dof.get_tria().n_active_cells()
+  //          << std::endl;
+  //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
+  //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl;
+
+  Vector<double> solution (dof.n_dofs());
+
+                                // create vector with random entries
+  for (unsigned int i=0; i<dof.n_dofs(); ++i)
+    {
+      if(constraints.is_constrained(i))
+        continue;
+      const double entry = rand()/(double)RAND_MAX;
+      solution(i) = entry;
+    }
+
+  constraints.distribute(solution);
+
+  FESystem<dim> fe_sys(dof.get_fe(), dim);
+  DoFHandler<dim> dofh_eulerian(dof.get_tria());
+  dofh_eulerian.distribute_dofs(fe_sys);
+
+  MatrixFree<dim,double> mf_data;
+  Vector<double> shift(dofh_eulerian.n_dofs());
+  for (unsigned int i=0; i<2; ++i)
+    {
+      if (i == 1)
+        {
+          shift(0) = 0.121;
+          shift(1) = -0.005;
+        }
+      MappingQEulerian<dim> mapping(2,shift,dofh_eulerian);
+      
+      {
+        const QGauss<1> quad (fe_degree+1);
+        typename MatrixFree<dim,double>::AdditionalData data;
+        data.tasks_parallel_scheme = MatrixFree<dim,double>::AdditionalData::none;
+        data.mapping_update_flags = update_gradients | update_second_derivatives;
+        if (i==1)
+          data.initialize_indices = false;
+        mf_data.reinit (mapping, dof, constraints, quad, data);
+      }
+      
+      MatrixFreeTest<dim,fe_degree,fe_degree+1> mf (mf_data, mapping);
+      mf.test_functions(solution);
+    }
+}
diff --git a/tests/matrix_free/update_mapping_only/cmp/generic b/tests/matrix_free/update_mapping_only/cmp/generic
new file mode 100644 (file)
index 0000000..8586990
--- /dev/null
@@ -0,0 +1,79 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(2)
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(4)
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Error function values: 0
+DEAL:3d::Error function gradients: 0
+DEAL:3d::Error function Laplacians: 0
+DEAL:3d::Error function diagonal of Hessian: 0
+DEAL:3d::Error function Hessians: 0
+DEAL:3d::
+DEAL:3d::Error function values: 0
+DEAL:3d::Error function gradients: 0
+DEAL:3d::Error function Laplacians: 0
+DEAL:3d::Error function diagonal of Hessian: 0
+DEAL:3d::Error function Hessians: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Error function values: 0
+DEAL:3d::Error function gradients: 0
+DEAL:3d::Error function Laplacians: 0
+DEAL:3d::Error function diagonal of Hessian: 0
+DEAL:3d::Error function Hessians: 0
+DEAL:3d::
+DEAL:3d::Error function values: 0
+DEAL:3d::Error function gradients: 0
+DEAL:3d::Error function Laplacians: 0
+DEAL:3d::Error function diagonal of Hessian: 0
+DEAL:3d::Error function Hessians: 0
+DEAL:3d::

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.