]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add tests for AD framework: ADOL-C
authorJean-Paul Pelteret <jppelteret@gmail.com>
Fri, 22 Mar 2019 09:09:47 +0000 (10:09 +0100)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Sat, 23 Mar 2019 17:40:30 +0000 (18:40 +0100)
122 files changed:
tests/adolc/helper_retaping_01.cc [new file with mode: 0644]
tests/adolc/helper_retaping_01.output [new file with mode: 0644]
tests/adolc/helper_retaping_02.cc [new file with mode: 0644]
tests/adolc/helper_retaping_02.output [new file with mode: 0644]
tests/adolc/helper_retaping_03.cc [new file with mode: 0644]
tests/adolc/helper_retaping_03.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_01_1.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_01_1.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_01_2.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_01_2.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_02_1.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_02_1.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_02_2.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_02_2.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_03_1.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_03_1.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_03_2.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_03_2.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_04_1.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_04_1.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_04_2.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_04_2.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_05_1.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_05_1.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_05_2.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_05_2.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_06_1.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_06_1.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_06_2.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_06_2.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_07_1.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_07_1.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_07_2.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_07_2.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_08_1.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_08_1.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_08_2.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_2_components_08_2.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_3_components_01_1.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_3_components_01_1.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_3_components_01_1.with_mpi=true.mpirun=2.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_3_components_01_2.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_3_components_01_2.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_3_components_01_2.with_mpi=true.mpirun=2.output [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_3_components_01_tapeless_only_2.cc [new file with mode: 0644]
tests/adolc/helper_scalar_coupled_3_components_01_tapeless_only_2.output [new file with mode: 0644]
tests/adolc/helper_scalar_failure_01.cc [new file with mode: 0644]
tests/adolc/helper_scalar_failure_01.output [new file with mode: 0644]
tests/adolc/helper_scalar_failure_02.cc [new file with mode: 0644]
tests/adolc/helper_scalar_failure_02.output [new file with mode: 0644]
tests/adolc/helper_scalar_failure_03.cc [new file with mode: 0644]
tests/adolc/helper_scalar_failure_03.output [new file with mode: 0644]
tests/adolc/helper_scalar_failure_04.cc [new file with mode: 0644]
tests/adolc/helper_scalar_failure_04.output [new file with mode: 0644]
tests/adolc/helper_scalar_failure_05.cc [new file with mode: 0644]
tests/adolc/helper_scalar_failure_05.output [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_01_1.cc [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_01_1.output [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_01_2.cc [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_01_2.output [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_02_1.cc [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_02_1.output [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_02_2.cc [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_02_2.output [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_03_1.cc [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_03_1.output [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_03_2.cc [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_03_2.output [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_04_1.cc [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_04_1.output [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_04_2.cc [new file with mode: 0644]
tests/adolc/helper_scalar_single_component_04_2.output [new file with mode: 0644]
tests/adolc/helper_tape_index_01.cc [new file with mode: 0644]
tests/adolc/helper_tape_index_01.output [new file with mode: 0644]
tests/adolc/helper_vector_1_indep_1_dep_vars_01_1.cc [new file with mode: 0644]
tests/adolc/helper_vector_1_indep_1_dep_vars_01_1.output [new file with mode: 0644]
tests/adolc/helper_vector_1_indep_1_dep_vars_01_2.cc [new file with mode: 0644]
tests/adolc/helper_vector_1_indep_1_dep_vars_01_2.output [new file with mode: 0644]
tests/adolc/helper_vector_1_indep_2_dep_vars_01_1.cc [new file with mode: 0644]
tests/adolc/helper_vector_1_indep_2_dep_vars_01_1.output [new file with mode: 0644]
tests/adolc/helper_vector_1_indep_2_dep_vars_01_2.cc [new file with mode: 0644]
tests/adolc/helper_vector_1_indep_2_dep_vars_01_2.output [new file with mode: 0644]
tests/adolc/helper_vector_2_indep_1_dep_vars_01_1.cc [new file with mode: 0644]
tests/adolc/helper_vector_2_indep_1_dep_vars_01_1.output [new file with mode: 0644]
tests/adolc/helper_vector_2_indep_1_dep_vars_01_2.cc [new file with mode: 0644]
tests/adolc/helper_vector_2_indep_1_dep_vars_01_2.output [new file with mode: 0644]
tests/adolc/helper_vector_2_indep_2_dep_vars_01_1.cc [new file with mode: 0644]
tests/adolc/helper_vector_2_indep_2_dep_vars_01_1.output [new file with mode: 0644]
tests/adolc/helper_vector_2_indep_2_dep_vars_01_2.cc [new file with mode: 0644]
tests/adolc/helper_vector_2_indep_2_dep_vars_01_2.output [new file with mode: 0644]
tests/adolc/step-44-helper_scalar_01_1.cc [new file with mode: 0644]
tests/adolc/step-44-helper_scalar_01_1.output [new file with mode: 0644]
tests/adolc/step-44-helper_vector_01_1.cc [new file with mode: 0644]
tests/adolc/step-44-helper_vector_01_1.output [new file with mode: 0644]
tests/adolc/step-44-helper_vector_01_2.cc [new file with mode: 0644]
tests/adolc/step-44-helper_vector_01_2.output [new file with mode: 0644]
tests/adolc/step-44.cc [new file with mode: 0644]
tests/adolc/step-44.output [new file with mode: 0644]
tests/adolc/symmetric_tensor_functions_02_1.cc [new file with mode: 0644]
tests/adolc/symmetric_tensor_functions_02_1.output [new file with mode: 0644]
tests/adolc/symmetric_tensor_functions_02_2.cc [new file with mode: 0644]
tests/adolc/symmetric_tensor_functions_02_2.output [new file with mode: 0644]
tests/adolc/symmetric_tensor_functions_03_1.cc [new file with mode: 0644]
tests/adolc/symmetric_tensor_functions_03_1.with_adolc_with_advanced_branching=false.output [new file with mode: 0644]
tests/adolc/symmetric_tensor_functions_03_1.with_adolc_with_advanced_branching=true.expect=run.output [new file with mode: 0644]
tests/adolc/symmetric_tensor_functions_03_2.cc [new file with mode: 0644]
tests/adolc/symmetric_tensor_functions_03_2.output [new file with mode: 0644]
tests/adolc/symmetric_tensor_functions_04_1.cc [new file with mode: 0644]
tests/adolc/symmetric_tensor_functions_04_1.with_adolc_with_advanced_branching=false.output [new file with mode: 0644]
tests/adolc/symmetric_tensor_functions_04_1.with_adolc_with_advanced_branching=true.expect=run.output [new file with mode: 0644]
tests/adolc/symmetric_tensor_functions_04_2.cc [new file with mode: 0644]
tests/adolc/symmetric_tensor_functions_04_2.output [new file with mode: 0644]
tests/adolc/tensor_functions_02_1.cc [new file with mode: 0644]
tests/adolc/tensor_functions_02_1.output [new file with mode: 0644]
tests/adolc/tensor_functions_02_2.cc [new file with mode: 0644]
tests/adolc/tensor_functions_02_2.output [new file with mode: 0644]
tests/adolc/tensor_ops_01.cc [new file with mode: 0644]
tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=false.output [new file with mode: 0644]
tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.expect=run.debug.output [new file with mode: 0644]
tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.release.output [new file with mode: 0644]
tests/adolc/tensor_ops_02.cc [new file with mode: 0644]
tests/adolc/tensor_ops_02.output [new file with mode: 0644]

diff --git a/tests/adolc/helper_retaping_01.cc b/tests/adolc/helper_retaping_01.cc
new file mode 100644 (file)
index 0000000..a119926
--- /dev/null
@@ -0,0 +1,297 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Check that a tape can be recorded and reused. In particular, we check
+// that the state of the ADHelper can flipped from "record" to "reuse"
+// and vice versa at any point in a simulation. Note that here the ADHelper
+// will not have a persistent lifetime, so we expect that ADOL-C will keep
+// any stored tape data as this object goes out of scope and is reconstructed
+// at a later evaluation stage.
+//
+// Based off of helper_scalar_coupled_2_components_01.h, this test the
+// evaluation of a coupled system (scalar + scalar components)
+// using a helper class.
+//
+// AD number type: ADOL-C taped
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/differentiation/ad.h>
+
+#include <deal.II/fe/fe_values_extractors.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+using namespace dealii;
+namespace AD = dealii::Differentiation::AD;
+
+// Function and its derivatives
+template <int dim, typename NumberType>
+struct FunctionsTestScalarScalarCoupled
+{
+  static NumberType
+  psi(const NumberType &s1, const NumberType &s2)
+  {
+    return 2.0 * std::pow(s1, 4) * std::pow(s2, 3);
+  };
+
+  static NumberType
+  dpsi_ds1(const NumberType &s1, const NumberType &s2)
+  {
+    return 8.0 * std::pow(s1, 3) * std::pow(s2, 3);
+  };
+
+  static NumberType
+  dpsi_ds2(const NumberType &s1, const NumberType &s2)
+  {
+    return 6.0 * std::pow(s1, 4) * std::pow(s2, 2);
+  };
+
+  static NumberType
+  d2psi_ds1_ds1(const NumberType &s1, const NumberType &s2)
+  {
+    return 24.0 * std::pow(s1, 2) * std::pow(s2, 3);
+  };
+
+  static NumberType
+  d2psi_ds2_ds1(const NumberType &s1, const NumberType &s2)
+  {
+    return 24.0 * std::pow(s1, 3) * std::pow(s2, 2);
+  };
+
+  static NumberType
+  d2psi_ds1_ds2(const NumberType &s1, const NumberType &s2)
+  {
+    return d2psi_ds2_ds1(s1, s2);
+  };
+
+  static NumberType
+  d2psi_ds2_ds2(const NumberType &s1, const NumberType &s2)
+  {
+    return 12.0 * std::pow(s1, 4) * std::pow(s2, 1);
+  };
+};
+
+template <int dim, typename number_t, enum AD::NumberTypes ad_type_code>
+void
+test_scalar_scalar_coupled_one_run(const number_t s1,
+                                   const number_t s2,
+                                   const bool     overwrite_tape)
+{
+  typedef AD::ScalarFunction<dim, ad_type_code, number_t> ADHelper;
+  typedef typename ADHelper::ad_type                      ADNumberType;
+  typedef typename ADHelper::scalar_type                  ScalarNumberType;
+
+  std::cout << "*** Test variables: Scalar + Scalar (coupled), "
+            << "dim = " << Utilities::to_string(dim) << ", "
+            << "Type code: " << static_cast<int>(ad_type_code) << ". "
+            << "Overwrite tape: " << std::boolalpha << overwrite_tape
+            << std::endl;
+
+  // Values computed from the AD energy function
+  ScalarNumberType             psi;
+  Vector<ScalarNumberType>     Dpsi;
+  FullMatrix<ScalarNumberType> D2psi;
+
+  // Function and its derivatives
+  typedef FunctionsTestScalarScalarCoupled<dim, ADNumberType> func_ad;
+
+  // Setup the variable components and choose a value at which to
+  // evaluate the tape
+  const FEValuesExtractors::Scalar s1_dof(0);
+  const FEValuesExtractors::Scalar s2_dof(1);
+  const unsigned int               n_AD_components = 2;
+  ADHelper                         ad_helper(n_AD_components);
+  if (overwrite_tape)
+    ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the
+                                       // default values
+
+  // Configure tape
+  const int  tape_no = 1;
+  const bool is_recording =
+    ad_helper.start_recording_operations(tape_no /*material_id*/,
+                                         overwrite_tape,
+                                         true /*keep*/);
+  if (is_recording == true)
+    {
+      ad_helper.register_independent_variable(s1, s1_dof);
+      ad_helper.register_independent_variable(s2, s2_dof);
+
+      const ADNumberType s1_ad = ad_helper.get_sensitive_variables(s1_dof);
+      const ADNumberType s2_ad = ad_helper.get_sensitive_variables(s2_dof);
+
+      const ADNumberType psi(func_ad::psi(s1_ad, s2_ad));
+
+      ad_helper.register_dependent_variable(psi);
+      ad_helper.stop_recording_operations(false /*write_tapes_to_file*/);
+
+      std::cout << "Taped data..." << std::endl;
+      std::cout << "independent variable values: " << std::flush;
+      ad_helper.print_values(std::cout);
+      std::cout << "s1_ad: " << s1_ad << std::endl;
+      std::cout << "s2_ad: " << s2_ad << std::endl;
+      std::cout << "psi: " << psi << std::endl;
+      std::cout << std::endl;
+    }
+  else
+    {
+      std::cout
+        << "Using tape with different values for independent variables..."
+        << std::endl;
+      // Set a new evaluation point
+      ad_helper.activate_recorded_tape(tape_no);
+      ad_helper.set_independent_variable(s1, s1_dof);
+      ad_helper.set_independent_variable(s2, s2_dof);
+    }
+
+  std::cout << "independent variable values: " << std::flush;
+  ad_helper.print_values(std::cout);
+
+  // Compute the function value, gradient and hessian for the new evaluation
+  // point
+  psi = ad_helper.compute_value();
+  ad_helper.compute_gradient(Dpsi);
+  if (AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2)
+    {
+      ad_helper.compute_hessian(D2psi);
+    }
+
+  // Output the full stored function, gradient vector and hessian matrix
+  std::cout << "psi: " << psi << std::endl;
+  std::cout << "Dpsi: \n";
+  Dpsi.print(std::cout);
+  if (AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2)
+    {
+      std::cout << "D2psi: \n";
+      D2psi.print_formatted(std::cout, 3, true, 0, "0.0");
+    }
+
+  // Extract components of the solution
+  const ScalarNumberType dpsi_ds1 =
+    ad_helper.extract_gradient_component(Dpsi, s1_dof);
+  const ScalarNumberType dpsi_ds2 =
+    ad_helper.extract_gradient_component(Dpsi, s2_dof);
+  std::cout << "extracted Dpsi (s1): " << dpsi_ds1 << "\n"
+            << "extracted Dpsi (s2): " << dpsi_ds2 << "\n";
+
+  // Verify the result
+  typedef FunctionsTestScalarScalarCoupled<dim, ScalarNumberType> func;
+  static const ScalarNumberType                                   tol =
+    1e8 * std::numeric_limits<ScalarNumberType>::epsilon();
+  std::cout << "dpsi_ds1:            " << dpsi_ds1 << std::endl;
+  std::cout << "func::dpsi_ds1(s1,s2): " << func::dpsi_ds1(s1, s2) << std::endl;
+  std::cout << "diff: " << std::abs(dpsi_ds1 - func::dpsi_ds1(s1, s2))
+            << std::endl;
+  std::cout << "dpsi_ds2:            " << dpsi_ds2 << std::endl;
+  std::cout << "func::dpsi_ds2(s1,s2): " << func::dpsi_ds2(s1, s2) << std::endl;
+  std::cout << "diff: " << std::abs(dpsi_ds2 - func::dpsi_ds2(s1, s2))
+            << std::endl;
+  Assert(std::abs(psi - func::psi(s1, s2)) < tol,
+         ExcMessage("No match for function value."));
+  Assert(std::abs(dpsi_ds1 - func::dpsi_ds1(s1, s2)) < tol,
+         ExcMessage("No match for first derivative."));
+  Assert(std::abs(dpsi_ds2 - func::dpsi_ds2(s1, s2)) < tol,
+         ExcMessage("No match for first derivative."));
+  if (AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2)
+    {
+      const ScalarNumberType d2psi_ds1_ds1 =
+        ad_helper.extract_hessian_component(D2psi, s1_dof, s1_dof);
+      const ScalarNumberType d2psi_ds2_ds1 =
+        ad_helper.extract_hessian_component(D2psi, s1_dof, s2_dof);
+      const ScalarNumberType d2psi_ds1_ds2 =
+        ad_helper.extract_hessian_component(D2psi, s2_dof, s1_dof);
+      const ScalarNumberType d2psi_ds2_ds2 =
+        ad_helper.extract_hessian_component(D2psi, s2_dof, s2_dof);
+      std::cout << "extracted D2psi (s1,s1): " << d2psi_ds1_ds1 << "\n"
+                << "extracted D2psi (s1,s2): " << d2psi_ds2_ds1 << "\n"
+                << "extracted D2psi (s2,s1): " << d2psi_ds1_ds2 << "\n"
+                << "extracted D2psi (s2,s2): " << d2psi_ds2_ds2 << "\n"
+                << std::endl;
+      std::cout << "d2psi_ds1_ds1:            " << d2psi_ds1_ds1 << std::endl;
+      std::cout << "func::d2psi_ds1_ds1(s1,s2): " << func::d2psi_ds1_ds1(s1, s2)
+                << std::endl;
+      std::cout << "diff: "
+                << std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1, s2))
+                << std::endl;
+      std::cout << "d2psi_ds2_ds1:            " << d2psi_ds2_ds1 << std::endl;
+      std::cout << "func::d2psi_ds2_ds1(s1,s2): " << func::d2psi_ds2_ds1(s1, s2)
+                << std::endl;
+      std::cout << "diff: "
+                << std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1, s2))
+                << std::endl;
+      std::cout << "d2psi_ds1_ds2:            " << d2psi_ds1_ds2 << std::endl;
+      std::cout << "func::d2psi_ds1_ds2(s1,s2): " << func::d2psi_ds1_ds2(s1, s2)
+                << std::endl;
+      std::cout << "diff: "
+                << std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1, s2))
+                << std::endl;
+      std::cout << "d2psi_ds2_ds2:            " << d2psi_ds2_ds2 << std::endl;
+      std::cout << "func::d2psi_ds2_ds2(s1,s2): " << func::d2psi_ds2_ds2(s1, s2)
+                << std::endl;
+      std::cout << "diff: "
+                << std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1, s2))
+                << std::endl;
+      Assert(std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+      Assert(std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+      Assert(std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+      Assert(std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+    }
+}
+
+
+
+template <int dim, typename number_t, enum AD::NumberTypes ad_type_code>
+void
+test_scalar_scalar_coupled(const unsigned int n_runs,
+                           const unsigned int retaping_modulo)
+{
+  for (unsigned int r = 0; r < n_runs; ++r)
+    {
+      const bool     overwrite_tape = (r % retaping_modulo == 0);
+      const number_t s1             = 3.1 * (1.0 + r / 10.0);
+      const number_t s2             = 5.9 * (1.0 + r / 10.0);
+      test_scalar_scalar_coupled_one_run<3, number_t, ad_type_code>(
+        s1, s2, overwrite_tape);
+    }
+}
+
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_scalar_scalar_coupled<2, double, AD::NumberTypes::adolc_taped>(12, 4);
+    test_scalar_scalar_coupled<3, double, AD::NumberTypes::adolc_taped>(4, 2);
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_retaping_01.output b/tests/adolc/helper_retaping_01.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_retaping_02.cc b/tests/adolc/helper_retaping_02.cc
new file mode 100644 (file)
index 0000000..4855401
--- /dev/null
@@ -0,0 +1,300 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Check that a tape can be recorded and reused. In particular, we check
+// that the state of the ADHelper can flipped from "record" to "reuse"
+// and vice versa at any point in a simulation. As opposed to
+// helper_retaping_01, the ADHelper is now a static object (i.e. it has
+// persistent lifetime, as might be the case if an  ADHelper is stored
+// as quadrature point member data).
+//
+// Based off of helper_scalar_coupled_2_components_01.h, this test the
+// evaluation of a coupled system (scalar + scalar components)
+// using a helper class.
+//
+// AD number type: ADOL-C taped
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/differentiation/ad.h>
+
+#include <deal.II/fe/fe_values_extractors.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+using namespace dealii;
+namespace AD = dealii::Differentiation::AD;
+
+// Function and its derivatives
+template <int dim, typename NumberType>
+struct FunctionsTestScalarScalarCoupled
+{
+  static NumberType
+  psi(const NumberType &s1, const NumberType &s2)
+  {
+    return 2.0 * std::pow(s1, 4) * std::pow(s2, 3);
+  };
+
+  static NumberType
+  dpsi_ds1(const NumberType &s1, const NumberType &s2)
+  {
+    return 8.0 * std::pow(s1, 3) * std::pow(s2, 3);
+  };
+
+  static NumberType
+  dpsi_ds2(const NumberType &s1, const NumberType &s2)
+  {
+    return 6.0 * std::pow(s1, 4) * std::pow(s2, 2);
+  };
+
+  static NumberType
+  d2psi_ds1_ds1(const NumberType &s1, const NumberType &s2)
+  {
+    return 24.0 * std::pow(s1, 2) * std::pow(s2, 3);
+  };
+
+  static NumberType
+  d2psi_ds2_ds1(const NumberType &s1, const NumberType &s2)
+  {
+    return 24.0 * std::pow(s1, 3) * std::pow(s2, 2);
+  };
+
+  static NumberType
+  d2psi_ds1_ds2(const NumberType &s1, const NumberType &s2)
+  {
+    return d2psi_ds2_ds1(s1, s2);
+  };
+
+  static NumberType
+  d2psi_ds2_ds2(const NumberType &s1, const NumberType &s2)
+  {
+    return 12.0 * std::pow(s1, 4) * std::pow(s2, 1);
+  };
+};
+
+template <int dim, typename number_t, enum AD::NumberTypes ad_type_code>
+void
+test_scalar_scalar_coupled_one_run(const number_t s1,
+                                   const number_t s2,
+                                   const bool     overwrite_tape)
+{
+  typedef AD::ScalarFunction<dim, ad_type_code, number_t> ADHelper;
+  typedef typename ADHelper::ad_type                      ADNumberType;
+  typedef typename ADHelper::scalar_type                  ScalarNumberType;
+
+  std::cout << "*** Test variables: Scalar + Scalar (coupled), "
+            << "dim = " << Utilities::to_string(dim) << ", "
+            << "Type code: " << static_cast<int>(ad_type_code) << ". "
+            << "Overwrite tape: " << std::boolalpha << overwrite_tape
+            << std::endl;
+
+  // Values computed from the AD energy function
+  ScalarNumberType             psi;
+  Vector<ScalarNumberType>     Dpsi;
+  FullMatrix<ScalarNumberType> D2psi;
+
+  // Function and its derivatives
+  typedef FunctionsTestScalarScalarCoupled<dim, ADNumberType> func_ad;
+
+  // Setup the variable components and choose a value at which to
+  // evaluate the tape
+  const FEValuesExtractors::Scalar s1_dof(0);
+  const FEValuesExtractors::Scalar s2_dof(1);
+  const unsigned int               n_AD_components = 2;
+  static ADHelper                  ad_helper(n_AD_components);
+  if (overwrite_tape)
+    {
+      ad_helper.reset();
+      ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the
+                                         // default values
+    }
+
+  // Configure tape
+  const int  tape_no = 1;
+  const bool is_recording =
+    ad_helper.start_recording_operations(tape_no /*material_id*/,
+                                         overwrite_tape,
+                                         true /*keep*/);
+  if (is_recording == true)
+    {
+      ad_helper.register_independent_variable(s1, s1_dof);
+      ad_helper.register_independent_variable(s2, s2_dof);
+
+      const ADNumberType s1_ad = ad_helper.get_sensitive_variables(s1_dof);
+      const ADNumberType s2_ad = ad_helper.get_sensitive_variables(s2_dof);
+
+      const ADNumberType psi(func_ad::psi(s1_ad, s2_ad));
+
+      ad_helper.register_dependent_variable(psi);
+      ad_helper.stop_recording_operations(false /*write_tapes_to_file*/);
+
+      std::cout << "Taped data..." << std::endl;
+      std::cout << "independent variable values: " << std::flush;
+      ad_helper.print_values(std::cout);
+      std::cout << "s1_ad: " << s1_ad << std::endl;
+      std::cout << "s2_ad: " << s2_ad << std::endl;
+      std::cout << "psi: " << psi << std::endl;
+      std::cout << std::endl;
+    }
+  else
+    {
+      std::cout
+        << "Using tape with different values for independent variables..."
+        << std::endl;
+      // Set a new evaluation point
+      ad_helper.activate_recorded_tape(tape_no);
+      ad_helper.set_independent_variable(s1, s1_dof);
+      ad_helper.set_independent_variable(s2, s2_dof);
+    }
+
+  std::cout << "independent variable values: " << std::flush;
+  ad_helper.print_values(std::cout);
+
+  // Compute the function value, gradient and hessian for the new evaluation
+  // point
+  psi = ad_helper.compute_value();
+  ad_helper.compute_gradient(Dpsi);
+  if (AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2)
+    {
+      ad_helper.compute_hessian(D2psi);
+    }
+
+  // Output the full stored function, gradient vector and hessian matrix
+  std::cout << "psi: " << psi << std::endl;
+  std::cout << "Dpsi: \n";
+  Dpsi.print(std::cout);
+  if (AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2)
+    {
+      std::cout << "D2psi: \n";
+      D2psi.print_formatted(std::cout, 3, true, 0, "0.0");
+    }
+
+  // Extract components of the solution
+  const ScalarNumberType dpsi_ds1 =
+    ad_helper.extract_gradient_component(Dpsi, s1_dof);
+  const ScalarNumberType dpsi_ds2 =
+    ad_helper.extract_gradient_component(Dpsi, s2_dof);
+  std::cout << "extracted Dpsi (s1): " << dpsi_ds1 << "\n"
+            << "extracted Dpsi (s2): " << dpsi_ds2 << "\n";
+
+  // Verify the result
+  typedef FunctionsTestScalarScalarCoupled<dim, ScalarNumberType> func;
+  static const ScalarNumberType                                   tol =
+    1e8 * std::numeric_limits<ScalarNumberType>::epsilon();
+  std::cout << "dpsi_ds1:            " << dpsi_ds1 << std::endl;
+  std::cout << "func::dpsi_ds1(s1,s2): " << func::dpsi_ds1(s1, s2) << std::endl;
+  std::cout << "diff: " << std::abs(dpsi_ds1 - func::dpsi_ds1(s1, s2))
+            << std::endl;
+  std::cout << "dpsi_ds2:            " << dpsi_ds2 << std::endl;
+  std::cout << "func::dpsi_ds2(s1,s2): " << func::dpsi_ds2(s1, s2) << std::endl;
+  std::cout << "diff: " << std::abs(dpsi_ds2 - func::dpsi_ds2(s1, s2))
+            << std::endl;
+  Assert(std::abs(psi - func::psi(s1, s2)) < tol,
+         ExcMessage("No match for function value."));
+  Assert(std::abs(dpsi_ds1 - func::dpsi_ds1(s1, s2)) < tol,
+         ExcMessage("No match for first derivative."));
+  Assert(std::abs(dpsi_ds2 - func::dpsi_ds2(s1, s2)) < tol,
+         ExcMessage("No match for first derivative."));
+  if (AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2)
+    {
+      const ScalarNumberType d2psi_ds1_ds1 =
+        ad_helper.extract_hessian_component(D2psi, s1_dof, s1_dof);
+      const ScalarNumberType d2psi_ds2_ds1 =
+        ad_helper.extract_hessian_component(D2psi, s1_dof, s2_dof);
+      const ScalarNumberType d2psi_ds1_ds2 =
+        ad_helper.extract_hessian_component(D2psi, s2_dof, s1_dof);
+      const ScalarNumberType d2psi_ds2_ds2 =
+        ad_helper.extract_hessian_component(D2psi, s2_dof, s2_dof);
+      std::cout << "extracted D2psi (s1,s1): " << d2psi_ds1_ds1 << "\n"
+                << "extracted D2psi (s1,s2): " << d2psi_ds2_ds1 << "\n"
+                << "extracted D2psi (s2,s1): " << d2psi_ds1_ds2 << "\n"
+                << "extracted D2psi (s2,s2): " << d2psi_ds2_ds2 << "\n"
+                << std::endl;
+      std::cout << "d2psi_ds1_ds1:            " << d2psi_ds1_ds1 << std::endl;
+      std::cout << "func::d2psi_ds1_ds1(s1,s2): " << func::d2psi_ds1_ds1(s1, s2)
+                << std::endl;
+      std::cout << "diff: "
+                << std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1, s2))
+                << std::endl;
+      std::cout << "d2psi_ds2_ds1:            " << d2psi_ds2_ds1 << std::endl;
+      std::cout << "func::d2psi_ds2_ds1(s1,s2): " << func::d2psi_ds2_ds1(s1, s2)
+                << std::endl;
+      std::cout << "diff: "
+                << std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1, s2))
+                << std::endl;
+      std::cout << "d2psi_ds1_ds2:            " << d2psi_ds1_ds2 << std::endl;
+      std::cout << "func::d2psi_ds1_ds2(s1,s2): " << func::d2psi_ds1_ds2(s1, s2)
+                << std::endl;
+      std::cout << "diff: "
+                << std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1, s2))
+                << std::endl;
+      std::cout << "d2psi_ds2_ds2:            " << d2psi_ds2_ds2 << std::endl;
+      std::cout << "func::d2psi_ds2_ds2(s1,s2): " << func::d2psi_ds2_ds2(s1, s2)
+                << std::endl;
+      std::cout << "diff: "
+                << std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1, s2))
+                << std::endl;
+      Assert(std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+      Assert(std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+      Assert(std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+      Assert(std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+    }
+}
+
+
+
+template <int dim, typename number_t, enum AD::NumberTypes ad_type_code>
+void
+test_scalar_scalar_coupled(const unsigned int n_runs,
+                           const unsigned int retaping_modulo)
+{
+  for (unsigned int r = 0; r < n_runs; ++r)
+    {
+      const bool     overwrite_tape = (r % retaping_modulo == 0);
+      const number_t s1             = 3.1 * (1.0 + r / 10.0);
+      const number_t s2             = 5.9 * (1.0 + r / 10.0);
+      test_scalar_scalar_coupled_one_run<3, number_t, ad_type_code>(
+        s1, s2, overwrite_tape);
+    }
+}
+
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_scalar_scalar_coupled<2, double, AD::NumberTypes::adolc_taped>(12, 4);
+    test_scalar_scalar_coupled<3, double, AD::NumberTypes::adolc_taped>(4, 2);
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_retaping_02.output b/tests/adolc/helper_retaping_02.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_retaping_03.cc b/tests/adolc/helper_retaping_03.cc
new file mode 100644 (file)
index 0000000..9a98225
--- /dev/null
@@ -0,0 +1,304 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Check that a tape can be recorded and reused. Here we specifically check
+// the status of the is_registered_tape() flag, and what happens after all
+// tapes are cleared. We explicitly do not allow tapes to be overwritten.
+// Note that this also checks if we can safely query the existence of a
+// tape that has not yet been traced (i.e. it is non-existent).
+//
+// Based off of helper_scalar_coupled_2_components_01.h, this test the
+// evaluation of a coupled system (scalar + scalar components)
+// using a helper class.
+//
+// AD number type: ADOL-C taped
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/differentiation/ad.h>
+
+#include <deal.II/fe/fe_values_extractors.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+using namespace dealii;
+namespace AD = dealii::Differentiation::AD;
+
+// Function and its derivatives
+template <int dim, typename NumberType>
+struct FunctionsTestScalarScalarCoupled
+{
+  static NumberType
+  psi(const NumberType &s1, const NumberType &s2)
+  {
+    return 2.0 * std::pow(s1, 4) * std::pow(s2, 3);
+  };
+
+  static NumberType
+  dpsi_ds1(const NumberType &s1, const NumberType &s2)
+  {
+    return 8.0 * std::pow(s1, 3) * std::pow(s2, 3);
+  };
+
+  static NumberType
+  dpsi_ds2(const NumberType &s1, const NumberType &s2)
+  {
+    return 6.0 * std::pow(s1, 4) * std::pow(s2, 2);
+  };
+
+  static NumberType
+  d2psi_ds1_ds1(const NumberType &s1, const NumberType &s2)
+  {
+    return 24.0 * std::pow(s1, 2) * std::pow(s2, 3);
+  };
+
+  static NumberType
+  d2psi_ds2_ds1(const NumberType &s1, const NumberType &s2)
+  {
+    return 24.0 * std::pow(s1, 3) * std::pow(s2, 2);
+  };
+
+  static NumberType
+  d2psi_ds1_ds2(const NumberType &s1, const NumberType &s2)
+  {
+    return d2psi_ds2_ds1(s1, s2);
+  };
+
+  static NumberType
+  d2psi_ds2_ds2(const NumberType &s1, const NumberType &s2)
+  {
+    return 12.0 * std::pow(s1, 4) * std::pow(s2, 1);
+  };
+};
+
+template <int dim, typename number_t, enum AD::NumberTypes ad_type_code>
+void
+test_scalar_scalar_coupled_one_run(const number_t s1,
+                                   const number_t s2,
+                                   const bool     clear_tapes)
+{
+  typedef AD::ScalarFunction<dim, ad_type_code, number_t> ADHelper;
+  typedef typename ADHelper::ad_type                      ADNumberType;
+  typedef typename ADHelper::scalar_type                  ScalarNumberType;
+
+  std::cout << "*** Test variables: Scalar + Scalar (coupled), "
+            << "dim = " << Utilities::to_string(dim) << ", "
+            << "Type code: " << static_cast<int>(ad_type_code) << ". "
+            << "Clear tapes: " << std::boolalpha << clear_tapes << std::endl;
+
+  // Values computed from the AD energy function
+  ScalarNumberType             psi;
+  Vector<ScalarNumberType>     Dpsi;
+  FullMatrix<ScalarNumberType> D2psi;
+
+  // Function and its derivatives
+  typedef FunctionsTestScalarScalarCoupled<dim, ADNumberType> func_ad;
+
+  // Setup the variable components and choose a value at which to
+  // evaluate the tape
+  const FEValuesExtractors::Scalar s1_dof(0);
+  const FEValuesExtractors::Scalar s2_dof(1);
+  const unsigned int               n_AD_components = 2;
+  ADHelper                         ad_helper(n_AD_components);
+  if (clear_tapes)
+    {
+      ad_helper.reset();
+      ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the
+                                         // default values
+    }
+
+  // Configure tape
+  const int tape_no = 1;
+  deallog << "is_registered_tape? " << std::boolalpha
+          << ad_helper.is_registered_tape(tape_no) << std::endl;
+
+  const bool is_recording =
+    ad_helper.start_recording_operations(tape_no /*material_id*/,
+                                         false /*overwrite_tape*/,
+                                         true /*keep*/);
+  if (is_recording == true)
+    {
+      deallog << "Recording new tape" << std::endl;
+      ad_helper.register_independent_variable(s1, s1_dof);
+      ad_helper.register_independent_variable(s2, s2_dof);
+
+      const ADNumberType s1_ad = ad_helper.get_sensitive_variables(s1_dof);
+      const ADNumberType s2_ad = ad_helper.get_sensitive_variables(s2_dof);
+
+      const ADNumberType psi(func_ad::psi(s1_ad, s2_ad));
+
+      ad_helper.register_dependent_variable(psi);
+      ad_helper.stop_recording_operations(false /*write_tapes_to_file*/);
+
+      std::cout << "Taped data..." << std::endl;
+      std::cout << "independent variable values: " << std::flush;
+      ad_helper.print_values(std::cout);
+      std::cout << "s1_ad: " << s1_ad << std::endl;
+      std::cout << "s2_ad: " << s2_ad << std::endl;
+      std::cout << "psi: " << psi << std::endl;
+      std::cout << std::endl;
+    }
+  else
+    {
+      deallog << "Reusing prerecorded tape" << std::endl;
+      std::cout
+        << "Using tape with different values for independent variables..."
+        << std::endl;
+      // Set a new evaluation point
+      ad_helper.activate_recorded_tape(tape_no);
+      ad_helper.set_independent_variable(s1, s1_dof);
+      ad_helper.set_independent_variable(s2, s2_dof);
+    }
+
+  std::cout << "independent variable values: " << std::flush;
+  ad_helper.print_values(std::cout);
+
+  // Compute the function value, gradient and hessian for the new evaluation
+  // point
+  psi = ad_helper.compute_value();
+  ad_helper.compute_gradient(Dpsi);
+  if (AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2)
+    {
+      ad_helper.compute_hessian(D2psi);
+    }
+
+  // Output the full stored function, gradient vector and hessian matrix
+  std::cout << "psi: " << psi << std::endl;
+  std::cout << "Dpsi: \n";
+  Dpsi.print(std::cout);
+  if (AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2)
+    {
+      std::cout << "D2psi: \n";
+      D2psi.print_formatted(std::cout, 3, true, 0, "0.0");
+    }
+
+  // Extract components of the solution
+  const ScalarNumberType dpsi_ds1 =
+    ad_helper.extract_gradient_component(Dpsi, s1_dof);
+  const ScalarNumberType dpsi_ds2 =
+    ad_helper.extract_gradient_component(Dpsi, s2_dof);
+  std::cout << "extracted Dpsi (s1): " << dpsi_ds1 << "\n"
+            << "extracted Dpsi (s2): " << dpsi_ds2 << "\n";
+
+  // Verify the result
+  typedef FunctionsTestScalarScalarCoupled<dim, ScalarNumberType> func;
+  static const ScalarNumberType                                   tol =
+    1e8 * std::numeric_limits<ScalarNumberType>::epsilon();
+  std::cout << "dpsi_ds1:            " << dpsi_ds1 << std::endl;
+  std::cout << "func::dpsi_ds1(s1,s2): " << func::dpsi_ds1(s1, s2) << std::endl;
+  std::cout << "diff: " << std::abs(dpsi_ds1 - func::dpsi_ds1(s1, s2))
+            << std::endl;
+  std::cout << "dpsi_ds2:            " << dpsi_ds2 << std::endl;
+  std::cout << "func::dpsi_ds2(s1,s2): " << func::dpsi_ds2(s1, s2) << std::endl;
+  std::cout << "diff: " << std::abs(dpsi_ds2 - func::dpsi_ds2(s1, s2))
+            << std::endl;
+  Assert(std::abs(psi - func::psi(s1, s2)) < tol,
+         ExcMessage("No match for function value."));
+  Assert(std::abs(dpsi_ds1 - func::dpsi_ds1(s1, s2)) < tol,
+         ExcMessage("No match for first derivative."));
+  Assert(std::abs(dpsi_ds2 - func::dpsi_ds2(s1, s2)) < tol,
+         ExcMessage("No match for first derivative."));
+  if (AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2)
+    {
+      const ScalarNumberType d2psi_ds1_ds1 =
+        ad_helper.extract_hessian_component(D2psi, s1_dof, s1_dof);
+      const ScalarNumberType d2psi_ds2_ds1 =
+        ad_helper.extract_hessian_component(D2psi, s1_dof, s2_dof);
+      const ScalarNumberType d2psi_ds1_ds2 =
+        ad_helper.extract_hessian_component(D2psi, s2_dof, s1_dof);
+      const ScalarNumberType d2psi_ds2_ds2 =
+        ad_helper.extract_hessian_component(D2psi, s2_dof, s2_dof);
+      std::cout << "extracted D2psi (s1,s1): " << d2psi_ds1_ds1 << "\n"
+                << "extracted D2psi (s1,s2): " << d2psi_ds2_ds1 << "\n"
+                << "extracted D2psi (s2,s1): " << d2psi_ds1_ds2 << "\n"
+                << "extracted D2psi (s2,s2): " << d2psi_ds2_ds2 << "\n"
+                << std::endl;
+      std::cout << "d2psi_ds1_ds1:            " << d2psi_ds1_ds1 << std::endl;
+      std::cout << "func::d2psi_ds1_ds1(s1,s2): " << func::d2psi_ds1_ds1(s1, s2)
+                << std::endl;
+      std::cout << "diff: "
+                << std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1, s2))
+                << std::endl;
+      std::cout << "d2psi_ds2_ds1:            " << d2psi_ds2_ds1 << std::endl;
+      std::cout << "func::d2psi_ds2_ds1(s1,s2): " << func::d2psi_ds2_ds1(s1, s2)
+                << std::endl;
+      std::cout << "diff: "
+                << std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1, s2))
+                << std::endl;
+      std::cout << "d2psi_ds1_ds2:            " << d2psi_ds1_ds2 << std::endl;
+      std::cout << "func::d2psi_ds1_ds2(s1,s2): " << func::d2psi_ds1_ds2(s1, s2)
+                << std::endl;
+      std::cout << "diff: "
+                << std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1, s2))
+                << std::endl;
+      std::cout << "d2psi_ds2_ds2:            " << d2psi_ds2_ds2 << std::endl;
+      std::cout << "func::d2psi_ds2_ds2(s1,s2): " << func::d2psi_ds2_ds2(s1, s2)
+                << std::endl;
+      std::cout << "diff: "
+                << std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1, s2))
+                << std::endl;
+      Assert(std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+      Assert(std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+      Assert(std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+      Assert(std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+    }
+}
+
+
+
+template <int dim, typename number_t, enum AD::NumberTypes ad_type_code>
+void
+test_scalar_scalar_coupled(const unsigned int n_runs,
+                           const unsigned int clearing_modulo)
+{
+  for (unsigned int r = 0; r < n_runs; ++r)
+    {
+      const bool     clear_tape = (r % clearing_modulo == 0);
+      const number_t s1         = 3.1 * (1.0 + r / 10.0);
+      const number_t s2         = 5.9 * (1.0 + r / 10.0);
+      test_scalar_scalar_coupled_one_run<3, number_t, ad_type_code>(s1,
+                                                                    s2,
+                                                                    clear_tape);
+    }
+}
+
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_scalar_scalar_coupled<2, double, AD::NumberTypes::adolc_taped>(12, 4);
+    test_scalar_scalar_coupled<3, double, AD::NumberTypes::adolc_taped>(4, 2);
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_retaping_03.output b/tests/adolc/helper_retaping_03.output
new file mode 100644 (file)
index 0000000..1dffec2
--- /dev/null
@@ -0,0 +1,35 @@
+
+DEAL:Double::is_registered_tape? false
+DEAL:Double::Recording new tape
+DEAL:Double::is_registered_tape? true
+DEAL:Double::Reusing prerecorded tape
+DEAL:Double::is_registered_tape? true
+DEAL:Double::Reusing prerecorded tape
+DEAL:Double::is_registered_tape? true
+DEAL:Double::Reusing prerecorded tape
+DEAL:Double::is_registered_tape? false
+DEAL:Double::Recording new tape
+DEAL:Double::is_registered_tape? true
+DEAL:Double::Reusing prerecorded tape
+DEAL:Double::is_registered_tape? true
+DEAL:Double::Reusing prerecorded tape
+DEAL:Double::is_registered_tape? true
+DEAL:Double::Reusing prerecorded tape
+DEAL:Double::is_registered_tape? false
+DEAL:Double::Recording new tape
+DEAL:Double::is_registered_tape? true
+DEAL:Double::Reusing prerecorded tape
+DEAL:Double::is_registered_tape? true
+DEAL:Double::Reusing prerecorded tape
+DEAL:Double::is_registered_tape? true
+DEAL:Double::Reusing prerecorded tape
+DEAL:Double::is_registered_tape? false
+DEAL:Double::Recording new tape
+DEAL:Double::is_registered_tape? true
+DEAL:Double::Reusing prerecorded tape
+DEAL:Double::is_registered_tape? false
+DEAL:Double::Recording new tape
+DEAL:Double::is_registered_tape? true
+DEAL:Double::Reusing prerecorded tape
+DEAL:Double::OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_01_1.cc b/tests/adolc/helper_scalar_coupled_2_components_01_1.cc
new file mode 100644 (file)
index 0000000..08cb486
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (scalar + scalar components)
+// using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_01.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_scalar_scalar_coupled<2, double, AD::NumberTypes::adolc_taped>();
+    test_scalar_scalar_coupled<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_01_1.output b/tests/adolc/helper_scalar_coupled_2_components_01_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_01_2.cc b/tests/adolc/helper_scalar_coupled_2_components_01_2.cc
new file mode 100644 (file)
index 0000000..e86a3c2
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (scalar + scalar components)
+// using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_01.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_scalar_scalar_coupled<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_scalar_scalar_coupled<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_01_2.output b/tests/adolc/helper_scalar_coupled_2_components_01_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_02_1.cc b/tests/adolc/helper_scalar_coupled_2_components_02_1.cc
new file mode 100644 (file)
index 0000000..a7fc503
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (vector + scalar components)
+// using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_02.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_vector_scalar_coupled<2, double, AD::NumberTypes::adolc_taped>();
+    test_vector_scalar_coupled<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_02_1.output b/tests/adolc/helper_scalar_coupled_2_components_02_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_02_2.cc b/tests/adolc/helper_scalar_coupled_2_components_02_2.cc
new file mode 100644 (file)
index 0000000..a5ddb71
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (vector + scalar components)
+// using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_02.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_vector_scalar_coupled<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_vector_scalar_coupled<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_02_2.output b/tests/adolc/helper_scalar_coupled_2_components_02_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_03_1.cc b/tests/adolc/helper_scalar_coupled_2_components_03_1.cc
new file mode 100644 (file)
index 0000000..905a93f
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (tensor + scalar components)
+// using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_03.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_tensor_scalar_coupled<2, double, AD::NumberTypes::adolc_taped>();
+    test_tensor_scalar_coupled<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_03_1.output b/tests/adolc/helper_scalar_coupled_2_components_03_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_03_2.cc b/tests/adolc/helper_scalar_coupled_2_components_03_2.cc
new file mode 100644 (file)
index 0000000..67b3d5a
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (tensor + scalar components)
+// using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_03.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_tensor_scalar_coupled<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_tensor_scalar_coupled<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_03_2.output b/tests/adolc/helper_scalar_coupled_2_components_03_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_04_1.cc b/tests/adolc/helper_scalar_coupled_2_components_04_1.cc
new file mode 100644 (file)
index 0000000..f16578a
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (vector + vector components)
+// using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_04.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_vector_vector_coupled<2, double, AD::NumberTypes::adolc_taped>();
+    test_vector_vector_coupled<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_04_1.output b/tests/adolc/helper_scalar_coupled_2_components_04_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_04_2.cc b/tests/adolc/helper_scalar_coupled_2_components_04_2.cc
new file mode 100644 (file)
index 0000000..295128d
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (vector + vector components)
+// using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_04.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_vector_vector_coupled<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_vector_vector_coupled<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_04_2.output b/tests/adolc/helper_scalar_coupled_2_components_04_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_05_1.cc b/tests/adolc/helper_scalar_coupled_2_components_05_1.cc
new file mode 100644 (file)
index 0000000..bb60dc5
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (tensor + vector components)
+// using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_05.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_tensor_vector_coupled<2, double, AD::NumberTypes::adolc_taped>();
+    test_tensor_vector_coupled<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_05_1.output b/tests/adolc/helper_scalar_coupled_2_components_05_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_05_2.cc b/tests/adolc/helper_scalar_coupled_2_components_05_2.cc
new file mode 100644 (file)
index 0000000..eb197e8
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (tensor + vector components)
+// using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_05.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_tensor_vector_coupled<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_tensor_vector_coupled<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_05_2.output b/tests/adolc/helper_scalar_coupled_2_components_05_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_06_1.cc b/tests/adolc/helper_scalar_coupled_2_components_06_1.cc
new file mode 100644 (file)
index 0000000..1c49818
--- /dev/null
@@ -0,0 +1,43 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (symmetric tensor + vector components)
+// using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_06.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_symmetric_tensor_vector_coupled<2,
+                                         double,
+                                         AD::NumberTypes::adolc_taped>();
+    test_symmetric_tensor_vector_coupled<3,
+                                         double,
+                                         AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_06_1.output b/tests/adolc/helper_scalar_coupled_2_components_06_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_06_2.cc b/tests/adolc/helper_scalar_coupled_2_components_06_2.cc
new file mode 100644 (file)
index 0000000..9252106
--- /dev/null
@@ -0,0 +1,43 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (symmetric tensor + vector components)
+// using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_06.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_symmetric_tensor_vector_coupled<2,
+                                         double,
+                                         AD::NumberTypes::adolc_tapeless>();
+    test_symmetric_tensor_vector_coupled<3,
+                                         double,
+                                         AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_06_2.output b/tests/adolc/helper_scalar_coupled_2_components_06_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_07_1.cc b/tests/adolc/helper_scalar_coupled_2_components_07_1.cc
new file mode 100644 (file)
index 0000000..c3aabfb
--- /dev/null
@@ -0,0 +1,43 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (symmetric tensor + tensor components)
+// using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_07.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_symmetric_tensor_tensor_coupled<2,
+                                         double,
+                                         AD::NumberTypes::adolc_taped>();
+    test_symmetric_tensor_tensor_coupled<3,
+                                         double,
+                                         AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_07_1.output b/tests/adolc/helper_scalar_coupled_2_components_07_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_07_2.cc b/tests/adolc/helper_scalar_coupled_2_components_07_2.cc
new file mode 100644 (file)
index 0000000..3b9dc68
--- /dev/null
@@ -0,0 +1,43 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (symmetric tensor + tensor components)
+// using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_07.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_symmetric_tensor_tensor_coupled<2,
+                                         double,
+                                         AD::NumberTypes::adolc_tapeless>();
+    test_symmetric_tensor_tensor_coupled<3,
+                                         double,
+                                         AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_07_2.output b/tests/adolc/helper_scalar_coupled_2_components_07_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_08_1.cc b/tests/adolc/helper_scalar_coupled_2_components_08_1.cc
new file mode 100644 (file)
index 0000000..126f69b
--- /dev/null
@@ -0,0 +1,45 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (symmetric tensor + symmetric tensor
+// components) using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_08.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_symmetric_tensor_symmetric_tensor_coupled<
+      2,
+      double,
+      AD::NumberTypes::adolc_taped>();
+    test_symmetric_tensor_symmetric_tensor_coupled<
+      3,
+      double,
+      AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_08_1.output b/tests/adolc/helper_scalar_coupled_2_components_08_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_2_components_08_2.cc b/tests/adolc/helper_scalar_coupled_2_components_08_2.cc
new file mode 100644 (file)
index 0000000..2bc2449
--- /dev/null
@@ -0,0 +1,45 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (symmetric tensor + symmetric tensor
+// components) using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_scalar_coupled_2_components_08.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_symmetric_tensor_symmetric_tensor_coupled<
+      2,
+      double,
+      AD::NumberTypes::adolc_tapeless>();
+    test_symmetric_tensor_symmetric_tensor_coupled<
+      3,
+      double,
+      AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_2_components_08_2.output b/tests/adolc/helper_scalar_coupled_2_components_08_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_1.cc b/tests/adolc/helper_scalar_coupled_3_components_01_1.cc
new file mode 100644 (file)
index 0000000..db7bb57
--- /dev/null
@@ -0,0 +1,44 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (tensor + vector + scalar components)
+// using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_scalar_coupled_3_components_01.h"
+#include "../tests.h"
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    mpi_log;
+
+  deallog.push("Double");
+  {
+    test_tensor_vector_scalar_coupled<2,
+                                      double,
+                                      AD::NumberTypes::adolc_taped>();
+    test_tensor_vector_scalar_coupled<3,
+                                      double,
+                                      AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_1.output b/tests/adolc/helper_scalar_coupled_3_components_01_1.output
new file mode 100644 (file)
index 0000000..aeccfbd
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL:0:Double:OK
+DEAL:0:OK
diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_1.with_mpi=true.mpirun=2.output b/tests/adolc/helper_scalar_coupled_3_components_01_1.with_mpi=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..4448125
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL:0:Double::OK
+DEAL:0:OK
+
+DEAL:1:Double::OK
+DEAL:1:OK
\ No newline at end of file
diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_2.cc b/tests/adolc/helper_scalar_coupled_3_components_01_2.cc
new file mode 100644 (file)
index 0000000..8a5c950
--- /dev/null
@@ -0,0 +1,44 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (tensor + vector + scalar components)
+// using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_scalar_coupled_3_components_01.h"
+#include "../tests.h"
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    mpi_log;
+
+  deallog.push("Double");
+  {
+    test_tensor_vector_scalar_coupled<2,
+                                      double,
+                                      AD::NumberTypes::adolc_tapeless>();
+    test_tensor_vector_scalar_coupled<3,
+                                      double,
+                                      AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_2.output b/tests/adolc/helper_scalar_coupled_3_components_01_2.output
new file mode 100644 (file)
index 0000000..aeccfbd
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL:0:Double:OK
+DEAL:0:OK
diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_2.with_mpi=true.mpirun=2.output b/tests/adolc/helper_scalar_coupled_3_components_01_2.with_mpi=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..4448125
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL:0:Double::OK
+DEAL:0:OK
+
+DEAL:1:Double::OK
+DEAL:1:OK
\ No newline at end of file
diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_tapeless_only_2.cc b/tests/adolc/helper_scalar_coupled_3_components_01_tapeless_only_2.cc
new file mode 100644 (file)
index 0000000..f13ceca
--- /dev/null
@@ -0,0 +1,48 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a coupled system (tensor + vector + scalar components)
+// using a helper class.
+// This test is based off of helper_scalar_coupled_3_components_01.h, and checks
+// that everything still works in tapeless only mode (i.e. when the
+// start_recording_operations and stop_recording_operations calls are
+// removed).
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_scalar_coupled_3_components_01_tapeless_only.h"
+#include "../tests.h"
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    mpi_log;
+
+  deallog.push("Double");
+  {
+    test_tensor_vector_scalar_coupled<2,
+                                      double,
+                                      AD::NumberTypes::adolc_tapeless>();
+    test_tensor_vector_scalar_coupled<3,
+                                      double,
+                                      AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_coupled_3_components_01_tapeless_only_2.output b/tests/adolc/helper_scalar_coupled_3_components_01_tapeless_only_2.output
new file mode 100644 (file)
index 0000000..aeccfbd
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL:0:Double:OK
+DEAL:0:OK
diff --git a/tests/adolc/helper_scalar_failure_01.cc b/tests/adolc/helper_scalar_failure_01.cc
new file mode 100644 (file)
index 0000000..09aa0b7
--- /dev/null
@@ -0,0 +1,121 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test that some cases that we expect to fail do in fact do so:
+// Overlapping indices of independent variables
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/differentiation/ad.h>
+
+#include <deal.II/fe/fe_values_extractors.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iomanip>
+
+#include "../tests.h"
+
+using namespace dealii;
+namespace AD = dealii::Differentiation::AD;
+
+template <int dim, typename number_t, enum AD::NumberTypes ad_type_code>
+bool
+test_overlapping_ind_var_components()
+{
+  typedef AD::ScalarFunction<dim, ad_type_code, number_t> ADHelper;
+  typedef typename ADHelper::ad_type                      ADNumberType;
+
+  std::cout << "*** Test variables: Scalar + Scalar (coupled), "
+            << (AD::ADNumberTraits<ADNumberType>::is_taped == true ? "Taped" :
+                                                                     "Tapeless")
+            << std::endl;
+
+  try
+    {
+      // NOTE: Both extractors set to same index
+      const FEValuesExtractors::Scalar s1_dof(0);
+      const FEValuesExtractors::Scalar s2_dof(0);
+      const unsigned int               n_AD_components = 2;
+      ADHelper                         ad_helper(n_AD_components);
+      ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the
+                                         // default values
+
+      double s1 = 3.1;
+      double s2 = 5.9;
+
+      // Configure tape
+      const int  tape_no = 1;
+      const bool is_recording =
+        ad_helper.start_recording_operations(tape_no /*material_id*/,
+                                             true /*overwrite_tape*/,
+                                             true /*keep*/);
+
+      ad_helper.register_independent_variable(s1, s1_dof);
+      ad_helper.register_independent_variable(s2, s2_dof);
+
+      const ADNumberType s1_ad = ad_helper.get_sensitive_variables(s1_dof);
+      const ADNumberType s2_ad = ad_helper.get_sensitive_variables(s2_dof);
+
+      const ADNumberType psi(s1_ad * s2_ad);
+
+      ad_helper.register_dependent_variable(psi);
+      ad_helper.stop_recording_operations(false /*write_tapes_to_file*/);
+
+      std::cout << "Unexpected success!" << std::endl;
+      return true;
+    }
+  catch (...)
+    {
+      std::cout << "Expected failure." << std::endl;
+      return false;
+    }
+}
+
+int
+main()
+{
+  initlog();
+
+  deal_II_exceptions::disable_abort_on_exception();
+#if DEBUG
+  // Asserts should be triggered
+  const bool expected_result = false;
+#else
+  // User beware: Asserts ignored
+  const bool expected_result = true;
+#endif
+
+  const unsigned int dim = 2;
+  AssertThrow(
+    (test_overlapping_ind_var_components<dim,
+                                         double,
+                                         AD::NumberTypes::adolc_taped>() ==
+     expected_result),
+    ExcInternalError());
+  AssertThrow(
+    (test_overlapping_ind_var_components<dim,
+                                         double,
+                                         AD::NumberTypes::adolc_tapeless>() ==
+     expected_result),
+    ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_failure_01.output b/tests/adolc/helper_scalar_failure_01.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_failure_02.cc b/tests/adolc/helper_scalar_failure_02.cc
new file mode 100644 (file)
index 0000000..59757f1
--- /dev/null
@@ -0,0 +1,123 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test that some cases that we expect to fail do in fact do so:
+// Overlapping indices of independent variables
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/differentiation/ad.h>
+
+#include <deal.II/fe/fe_values_extractors.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iomanip>
+
+#include "../tests.h"
+
+using namespace dealii;
+namespace AD = dealii::Differentiation::AD;
+
+template <int dim, typename number_t, enum AD::NumberTypes ad_type_code>
+bool
+test_overlapping_ind_var_components()
+{
+  typedef AD::ScalarFunction<dim, ad_type_code, number_t> ADHelper;
+  typedef typename ADHelper::ad_type                      ADNumberType;
+
+  std::cout << "*** Test variables: Vector + Scalar (coupled), "
+            << (AD::ADNumberTraits<ADNumberType>::is_taped == true ? "Taped" :
+                                                                     "Tapeless")
+            << std::endl;
+
+  try
+    {
+      // NOTE: Second extractor overlaps first ont.
+      const FEValuesExtractors::Vector v1_dof(0);
+      const FEValuesExtractors::Scalar s2_dof(1);
+      const unsigned int               n_AD_components =
+        Tensor<1, dim>::n_independent_components + 1;
+      ADHelper ad_helper(n_AD_components);
+      ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the
+                                         // default values
+
+      Tensor<1, dim> v1;
+      double         s2 = 5.9;
+
+      // Configure tape
+      const int  tape_no = 1;
+      const bool is_recording =
+        ad_helper.start_recording_operations(tape_no /*material_id*/,
+                                             true /*overwrite_tape*/,
+                                             true /*keep*/);
+
+      ad_helper.register_independent_variable(v1, v1_dof);
+      ad_helper.register_independent_variable(s2, s2_dof);
+
+      const Tensor<1, dim, ADNumberType> v1_ad =
+        ad_helper.get_sensitive_variables(v1_dof);
+      const ADNumberType s2_ad = ad_helper.get_sensitive_variables(s2_dof);
+
+      const ADNumberType psi((v1_ad * v1_ad) * s2_ad);
+
+      ad_helper.register_dependent_variable(psi);
+      ad_helper.stop_recording_operations(false /*write_tapes_to_file*/);
+
+      std::cout << "Unexpected success!" << std::endl;
+      return true;
+    }
+  catch (...)
+    {
+      std::cout << "Expected failure." << std::endl;
+      return false;
+    }
+}
+
+int
+main()
+{
+  initlog();
+
+  deal_II_exceptions::disable_abort_on_exception();
+#if DEBUG
+  // Asserts should be triggered
+  const bool expected_result = false;
+#else
+  // User beware: Asserts ignored
+  const bool expected_result = true;
+#endif
+
+  const unsigned int dim = 2;
+  AssertThrow(
+    (test_overlapping_ind_var_components<dim,
+                                         double,
+                                         AD::NumberTypes::adolc_taped>() ==
+     expected_result),
+    ExcInternalError());
+  AssertThrow(
+    (test_overlapping_ind_var_components<dim,
+                                         double,
+                                         AD::NumberTypes::adolc_tapeless>() ==
+     expected_result),
+    ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_failure_02.output b/tests/adolc/helper_scalar_failure_02.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_failure_03.cc b/tests/adolc/helper_scalar_failure_03.cc
new file mode 100644 (file)
index 0000000..110c4e0
--- /dev/null
@@ -0,0 +1,124 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test that some cases that we expect to fail do in fact do so:
+// Overlapping indices of independent variables
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/differentiation/ad.h>
+
+#include <deal.II/fe/fe_values_extractors.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iomanip>
+
+#include "../tests.h"
+
+using namespace dealii;
+namespace AD = dealii::Differentiation::AD;
+
+template <int dim, typename number_t, enum AD::NumberTypes ad_type_code>
+bool
+test_overlapping_ind_var_components()
+{
+  typedef AD::ScalarFunction<dim, ad_type_code, number_t> ADHelper;
+  typedef typename ADHelper::ad_type                      ADNumberType;
+
+  std::cout << "*** Test variables: Tensor + Vector (coupled), "
+            << (AD::ADNumberTraits<ADNumberType>::is_taped == true ? "Taped" :
+                                                                     "Tapeless")
+            << std::endl;
+
+  try
+    {
+      // NOTE: Second extractor overlaps first ont.
+      const FEValuesExtractors::Tensor<2> t1_dof(0);
+      const FEValuesExtractors::Vector    v2_dof(dim + 1);
+      const unsigned int                  n_AD_components =
+        Tensor<2, dim>::n_independent_components +
+        Tensor<1, dim>::n_independent_components;
+      ADHelper ad_helper(n_AD_components);
+      ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the
+                                         // default values
+
+      Tensor<2, dim> t1;
+      Tensor<1, dim> v2;
+
+      // Configure tape
+      const int  tape_no = 1;
+      const bool is_recording =
+        ad_helper.start_recording_operations(tape_no /*material_id*/,
+                                             true /*overwrite_tape*/,
+                                             true /*keep*/);
+
+      ad_helper.register_independent_variable(t1, t1_dof);
+      ad_helper.register_independent_variable(v2, v2_dof);
+
+      const Tensor<2, dim, ADNumberType> t1_ad =
+        ad_helper.get_sensitive_variables(t1_dof);
+      const Tensor<1, dim, ADNumberType> v2_ad =
+        ad_helper.get_sensitive_variables(v2_dof);
+
+      const ADNumberType psi((v2_ad * v2_ad) * trace(t1_ad));
+
+      ad_helper.register_dependent_variable(psi);
+      ad_helper.stop_recording_operations(false /*write_tapes_to_file*/);
+
+      std::cout << "Unexpected success!" << std::endl;
+      return true;
+    }
+  catch (...)
+    {
+      std::cout << "Expected failure." << std::endl;
+      return false;
+    }
+}
+
+int
+main()
+{
+  initlog();
+
+  deal_II_exceptions::disable_abort_on_exception();
+#if DEBUG
+  // Asserts should be triggered
+  const bool expected_result = false;
+#else
+  // User beware: Asserts ignored
+  const bool expected_result = true;
+#endif
+
+  const unsigned int dim = 2;
+  AssertThrow(
+    (test_overlapping_ind_var_components<dim,
+                                         double,
+                                         AD::NumberTypes::adolc_taped>() ==
+     expected_result),
+    ExcInternalError());
+  AssertThrow(
+    (test_overlapping_ind_var_components<dim,
+                                         double,
+                                         AD::NumberTypes::adolc_tapeless>() ==
+     expected_result),
+    ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_failure_03.output b/tests/adolc/helper_scalar_failure_03.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_failure_04.cc b/tests/adolc/helper_scalar_failure_04.cc
new file mode 100644 (file)
index 0000000..bd757b8
--- /dev/null
@@ -0,0 +1,123 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test that some cases that we expect to fail do in fact do so:
+// Incorrect setting of number of independent variables
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/differentiation/ad.h>
+
+#include <deal.II/fe/fe_values_extractors.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iomanip>
+
+#include "../tests.h"
+
+using namespace dealii;
+namespace AD = dealii::Differentiation::AD;
+
+template <int dim, typename number_t, enum AD::NumberTypes ad_type_code>
+bool
+test_incorrect_no_of_ind_vars()
+{
+  typedef AD::ScalarFunction<dim, ad_type_code, number_t> ADHelper;
+  typedef typename ADHelper::ad_type                      ADNumberType;
+
+  std::cout << "*** Test variables: Scalar + Scalar (coupled), "
+            << (AD::ADNumberTraits<ADNumberType>::is_taped == true ? "Taped" :
+                                                                     "Tapeless")
+            << std::endl;
+
+  try
+    {
+      const FEValuesExtractors::Scalar s1_dof(0);
+      const FEValuesExtractors::Scalar s2_dof(1);
+      // NOTE: Only 1 independent variable set instead of 2
+      const unsigned int n_AD_components = 1;
+      ADHelper           ad_helper(n_AD_components);
+      ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the
+                                         // default values
+
+      double s1 = 3.1;
+      double s2 = 5.9;
+
+      // Configure tape
+      const int  tape_no = 1;
+      const bool is_recording =
+        ad_helper.start_recording_operations(tape_no /*material_id*/,
+                                             true /*overwrite_tape*/,
+                                             true /*keep*/);
+
+      ad_helper.register_independent_variable(s1, s1_dof);
+      ad_helper.register_independent_variable(s2, s2_dof);
+
+      const ADNumberType s1_ad = ad_helper.get_sensitive_variables(s1_dof);
+      const ADNumberType s2_ad = ad_helper.get_sensitive_variables(
+        s2_dof); // This line might lead to memory corruption in tapeless mode.
+
+      const ADNumberType psi(s1_ad * s2_ad);
+
+      ad_helper.register_dependent_variable(psi);
+      ad_helper.stop_recording_operations(false /*write_tapes_to_file*/);
+
+      std::cout << "Unexpected success!" << std::endl;
+      return true;
+    }
+  catch (...)
+    {
+      std::cout << "Expected failure." << std::endl;
+      return false;
+    }
+}
+
+int
+main()
+{
+  initlog();
+
+  deal_II_exceptions::disable_abort_on_exception();
+#if DEBUG
+  // Asserts should be triggered
+  const bool expected_result_taped    = false;
+  const bool expected_result_tapeless = false;
+#else
+  // User beware: Asserts ignored
+  const bool expected_result_taped    = true;
+  const bool expected_result_tapeless = false;
+#endif
+
+  const unsigned int dim = 2;
+  AssertThrow((test_incorrect_no_of_ind_vars<dim,
+                                             double,
+                                             AD::NumberTypes::adolc_taped>() ==
+               expected_result_taped),
+              ExcInternalError());
+  AssertThrow(
+    (test_incorrect_no_of_ind_vars<dim,
+                                   double,
+                                   AD::NumberTypes::adolc_tapeless>() ==
+     expected_result_tapeless),
+    ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_failure_04.output b/tests/adolc/helper_scalar_failure_04.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_failure_05.cc b/tests/adolc/helper_scalar_failure_05.cc
new file mode 100644 (file)
index 0000000..63e41a6
--- /dev/null
@@ -0,0 +1,124 @@
+
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test that some cases that we expect to fail do in fact do so:
+// Registering more than one dependent variable
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/differentiation/ad.h>
+
+#include <deal.II/fe/fe_values_extractors.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iomanip>
+
+#include "../tests.h"
+
+using namespace dealii;
+namespace AD = dealii::Differentiation::AD;
+
+template <int dim, typename number_t, enum AD::NumberTypes ad_type_code>
+bool
+test_register_multiple_dep_vars()
+{
+  typedef AD::ScalarFunction<dim, ad_type_code, number_t> ADHelper;
+  typedef typename ADHelper::ad_type                      ADNumberType;
+
+  std::cout << "*** Test variables: Scalar + Scalar (coupled), "
+            << (AD::ADNumberTraits<ADNumberType>::is_taped == true ? "Taped" :
+                                                                     "Tapeless")
+            << std::endl;
+
+  try
+    {
+      const FEValuesExtractors::Scalar s1_dof(0);
+      const FEValuesExtractors::Scalar s2_dof(1);
+      const unsigned int               n_AD_components = 2;
+      ADHelper                         ad_helper(n_AD_components);
+      ad_helper.set_tape_buffer_sizes(); // Increase the buffer size from the
+                                         // default values
+
+      double s1 = 3.1;
+      double s2 = 5.9;
+
+      // Configure tape
+      const int  tape_no = 1;
+      const bool is_recording =
+        ad_helper.start_recording_operations(tape_no /*material_id*/,
+                                             true /*overwrite_tape*/,
+                                             true /*keep*/);
+
+      ad_helper.register_independent_variable(s1, s1_dof);
+      ad_helper.register_independent_variable(s2, s2_dof);
+
+      const ADNumberType s1_ad = ad_helper.get_sensitive_variables(s1_dof);
+      const ADNumberType s2_ad = ad_helper.get_sensitive_variables(s2_dof);
+
+      const ADNumberType psi(s1_ad * s2_ad);
+
+      // NOTE: Should only be calling this function once for this helper class
+      ad_helper.register_dependent_variable(psi);
+      ad_helper.register_dependent_variable(psi); // Wrong!
+
+      ad_helper.stop_recording_operations(false /*write_tapes_to_file*/);
+
+      std::cout << "Unexpected success!" << std::endl;
+      return true;
+    }
+  catch (...)
+    {
+      std::cout << "Expected failure." << std::endl;
+      return false;
+    }
+}
+
+int
+main()
+{
+  initlog();
+
+  deal_II_exceptions::disable_abort_on_exception();
+#if DEBUG
+  // Asserts should be triggered
+  const bool expected_result = false;
+#else
+  // User beware: Asserts ignored
+  const bool expected_result = true;
+#endif
+
+  const unsigned int dim = 2;
+  AssertThrow(
+    (test_register_multiple_dep_vars<dim,
+                                     double,
+                                     AD::NumberTypes::adolc_taped>() ==
+     expected_result),
+    ExcInternalError());
+  AssertThrow(
+    (test_register_multiple_dep_vars<dim,
+                                     double,
+                                     AD::NumberTypes::adolc_tapeless>() ==
+     expected_result),
+    ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_failure_05.output b/tests/adolc/helper_scalar_failure_05.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_single_component_01_1.cc b/tests/adolc/helper_scalar_single_component_01_1.cc
new file mode 100644 (file)
index 0000000..0ef35c8
--- /dev/null
@@ -0,0 +1,38 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a single component (scalar) system using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_scalar_single_component_01.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_scalar<2, double, AD::NumberTypes::adolc_taped>();
+    test_scalar<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_single_component_01_1.output b/tests/adolc/helper_scalar_single_component_01_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_single_component_01_2.cc b/tests/adolc/helper_scalar_single_component_01_2.cc
new file mode 100644 (file)
index 0000000..73db246
--- /dev/null
@@ -0,0 +1,38 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a single component (scalar) system using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_scalar_single_component_01.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_scalar<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_scalar<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_single_component_01_2.output b/tests/adolc/helper_scalar_single_component_01_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_single_component_02_1.cc b/tests/adolc/helper_scalar_single_component_02_1.cc
new file mode 100644 (file)
index 0000000..c25c0b5
--- /dev/null
@@ -0,0 +1,38 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a single component (vector) system using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_scalar_single_component_02.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_vector<2, double, AD::NumberTypes::adolc_taped>();
+    test_vector<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_single_component_02_1.output b/tests/adolc/helper_scalar_single_component_02_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_single_component_02_2.cc b/tests/adolc/helper_scalar_single_component_02_2.cc
new file mode 100644 (file)
index 0000000..ef10d02
--- /dev/null
@@ -0,0 +1,38 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a single component (vector) system using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_scalar_single_component_02.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_vector<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_vector<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_single_component_02_2.output b/tests/adolc/helper_scalar_single_component_02_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_single_component_03_1.cc b/tests/adolc/helper_scalar_single_component_03_1.cc
new file mode 100644 (file)
index 0000000..981d452
--- /dev/null
@@ -0,0 +1,38 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a single component (tensor) system using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_scalar_single_component_03.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_tensor<2, double, AD::NumberTypes::adolc_taped>();
+    test_tensor<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_single_component_03_1.output b/tests/adolc/helper_scalar_single_component_03_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_single_component_03_2.cc b/tests/adolc/helper_scalar_single_component_03_2.cc
new file mode 100644 (file)
index 0000000..3118095
--- /dev/null
@@ -0,0 +1,38 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a single component (tensor) system using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_scalar_single_component_03.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_tensor<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_tensor<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_single_component_03_2.output b/tests/adolc/helper_scalar_single_component_03_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_single_component_04_1.cc b/tests/adolc/helper_scalar_single_component_04_1.cc
new file mode 100644 (file)
index 0000000..999ca53
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a single component (symmetric tensor) system using a helper
+// class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_scalar_single_component_04.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_symmetric_tensor<2, double, AD::NumberTypes::adolc_taped>();
+    test_symmetric_tensor<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_single_component_04_1.output b/tests/adolc/helper_scalar_single_component_04_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_scalar_single_component_04_2.cc b/tests/adolc/helper_scalar_single_component_04_2.cc
new file mode 100644 (file)
index 0000000..feebc1c
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a single component (symmetric tensor) system using a helper
+// class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_scalar_single_component_04.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_symmetric_tensor<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_symmetric_tensor<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_scalar_single_component_04_2.output b/tests/adolc/helper_scalar_single_component_04_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_tape_index_01.cc b/tests/adolc/helper_tape_index_01.cc
new file mode 100644 (file)
index 0000000..6245a90
--- /dev/null
@@ -0,0 +1,338 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Verify that the maximum number of indices that can be used to record tapes is
+// indeed that stored within the HelperBase class
+//
+// AD number type: ADOL-C taped
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/differentiation/ad.h>
+
+#include <deal.II/fe/fe_values_extractors.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+using namespace dealii;
+namespace AD = dealii::Differentiation::AD;
+
+// Function and its derivatives
+template <int dim, typename NumberType>
+struct FunctionsTestScalarScalarCoupled
+{
+  static NumberType
+  psi(const NumberType &s1, const NumberType &s2)
+  {
+    return 2.0 * std::pow(s1, 4) * std::pow(s2, 3);
+  };
+
+  static NumberType
+  dpsi_ds1(const NumberType &s1, const NumberType &s2)
+  {
+    return 8.0 * std::pow(s1, 3) * std::pow(s2, 3);
+  };
+
+  static NumberType
+  dpsi_ds2(const NumberType &s1, const NumberType &s2)
+  {
+    return 6.0 * std::pow(s1, 4) * std::pow(s2, 2);
+  };
+
+  static NumberType
+  d2psi_ds1_ds1(const NumberType &s1, const NumberType &s2)
+  {
+    return 24.0 * std::pow(s1, 2) * std::pow(s2, 3);
+  };
+
+  static NumberType
+  d2psi_ds2_ds1(const NumberType &s1, const NumberType &s2)
+  {
+    return 24.0 * std::pow(s1, 3) * std::pow(s2, 2);
+  };
+
+  static NumberType
+  d2psi_ds1_ds2(const NumberType &s1, const NumberType &s2)
+  {
+    return d2psi_ds2_ds1(s1, s2);
+  };
+
+  static NumberType
+  d2psi_ds2_ds2(const NumberType &s1, const NumberType &s2)
+  {
+    return 12.0 * std::pow(s1, 4) * std::pow(s2, 1);
+  };
+};
+
+template <int dim,
+          typename number_t,
+          enum AD::NumberTypes ad_type_code,
+          typename tape_index_t>
+void
+test_scalar_scalar_coupled_one_run(const number_t     s1,
+                                   const number_t     s2,
+                                   const tape_index_t tape_no)
+{
+  typedef AD::ScalarFunction<dim, ad_type_code, number_t> ADHelper;
+  typedef typename ADHelper::ad_type                      ADNumberType;
+  typedef typename ADHelper::scalar_type                  ScalarNumberType;
+
+  std::cout
+    //      << "*** Test variables: Scalar + Scalar (coupled), "
+    //      << "dim = " << Utilities::to_string(dim) << ", "
+    //      << "Type code: " << static_cast<int>(ad_type_code) << ". "
+    << "Tape number: " << tape_no << std::endl;
+
+  // Values computed from the AD energy function
+  ScalarNumberType             psi;
+  Vector<ScalarNumberType>     Dpsi;
+  FullMatrix<ScalarNumberType> D2psi;
+
+  // Function and its derivatives
+  typedef FunctionsTestScalarScalarCoupled<dim, ADNumberType> func_ad;
+
+  // Setup the variable components and choose a value at which to
+  // evaluate the tape
+  const FEValuesExtractors::Scalar s1_dof(0);
+  const FEValuesExtractors::Scalar s2_dof(1);
+  const unsigned int               n_AD_components = 2;
+  ADHelper                         ad_helper(n_AD_components);
+  ad_helper.set_tape_buffer_sizes();     // Increase the buffer size from the
+                                         // default values
+  const bool write_tapes_to_file = true; // This flag will help check whether
+                                         // tape indices are universally valid
+
+  // Configure tape
+  const bool is_recording =
+    ad_helper.start_recording_operations(tape_no,
+                                         false /*overwrite_tape*/,
+                                         true /*keep*/);
+  if (is_recording == true)
+    {
+      ad_helper.register_independent_variable(s1, s1_dof);
+      ad_helper.register_independent_variable(s2, s2_dof);
+
+      const ADNumberType s1_ad = ad_helper.get_sensitive_variables(s1_dof);
+      const ADNumberType s2_ad = ad_helper.get_sensitive_variables(s2_dof);
+
+      const ADNumberType psi(func_ad::psi(s1_ad, s2_ad));
+
+      ad_helper.register_dependent_variable(psi);
+      ad_helper.stop_recording_operations(write_tapes_to_file);
+
+      //      std::cout << "Taped data..." << std::endl;
+      //      std::cout << "independent variable values: " << std::flush;
+      //      ad_helper.print_values(std::cout);
+      //      std::cout << "s1_ad: " << s1_ad << std::endl;
+      //      std::cout << "s2_ad: " << s2_ad << std::endl;
+      //      std::cout << "psi: " << psi << std::endl;
+      //      std::cout << std::endl;
+    }
+  else
+    {
+      //      std::cout << "Using tape with different values for independent
+      //      variables..." << std::endl;
+      // Set a new evaluation point
+      ad_helper.activate_recorded_tape(tape_no);
+      ad_helper.set_independent_variable(s1, s1_dof);
+      ad_helper.set_independent_variable(s2, s2_dof);
+    }
+
+  //  std::cout << "independent variable values: " << std::flush;
+  //  ad_helper.print_values(std::cout);
+
+  // Compute the function value, gradient and hessian for the new evaluation
+  // point
+  psi = ad_helper.compute_value();
+  ad_helper.compute_gradient(Dpsi);
+  if (AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2)
+    {
+      ad_helper.compute_hessian(D2psi);
+    }
+
+  // Output the full stored function, gradient vector and hessian matrix
+  //  std::cout << "psi: " << psi << std::endl;
+  //  std::cout << "Dpsi: \n";
+  //  Dpsi.print(std::cout);
+  //  if (AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2)
+  //    {
+  //      std::cout << "D2psi: \n";
+  //      D2psi.print_formatted(std::cout,3,true,0,"0.0");
+  //    }
+
+  // Extract components of the solution
+  const ScalarNumberType dpsi_ds1 =
+    ad_helper.extract_gradient_component(Dpsi, s1_dof);
+  const ScalarNumberType dpsi_ds2 =
+    ad_helper.extract_gradient_component(Dpsi, s2_dof);
+  //  std::cout
+  //      << "extracted Dpsi (s1): " << dpsi_ds1 << "\n"
+  //      << "extracted Dpsi (s2): " << dpsi_ds2 << "\n";
+
+  // Verify the result
+  typedef FunctionsTestScalarScalarCoupled<dim, ScalarNumberType> func;
+  static const ScalarNumberType                                   tol =
+    1e8 * std::numeric_limits<ScalarNumberType>::epsilon();
+  //  std::cout << "dpsi_ds1:            " << dpsi_ds1 << std::endl;
+  //  std::cout << "func::dpsi_ds1(s1,s2): " << func::dpsi_ds1(s1,s2) <<
+  //  std::endl; std::cout << "diff: " << std::abs(dpsi_ds1 -
+  //  func::dpsi_ds1(s1,s2)) << std::endl; std::cout << "dpsi_ds2:            "
+  //  << dpsi_ds2 << std::endl; std::cout << "func::dpsi_ds2(s1,s2): " <<
+  //  func::dpsi_ds2(s1,s2) << std::endl; std::cout << "diff: " <<
+  //  std::abs(dpsi_ds2 - func::dpsi_ds2(s1,s2)) << std::endl;
+  Assert(std::abs(psi - func::psi(s1, s2)) < tol,
+         ExcMessage("No match for function value."));
+  Assert(std::abs(dpsi_ds1 - func::dpsi_ds1(s1, s2)) < tol,
+         ExcMessage("No match for first derivative."));
+  Assert(std::abs(dpsi_ds2 - func::dpsi_ds2(s1, s2)) < tol,
+         ExcMessage("No match for first derivative."));
+  if (AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2)
+    {
+      const ScalarNumberType d2psi_ds1_ds1 =
+        ad_helper.extract_hessian_component(D2psi, s1_dof, s1_dof);
+      const ScalarNumberType d2psi_ds2_ds1 =
+        ad_helper.extract_hessian_component(D2psi, s1_dof, s2_dof);
+      const ScalarNumberType d2psi_ds1_ds2 =
+        ad_helper.extract_hessian_component(D2psi, s2_dof, s1_dof);
+      const ScalarNumberType d2psi_ds2_ds2 =
+        ad_helper.extract_hessian_component(D2psi, s2_dof, s2_dof);
+      //      std::cout
+      //          << "extracted D2psi (s1,s1): " << d2psi_ds1_ds1 << "\n"
+      //          << "extracted D2psi (s1,s2): " << d2psi_ds2_ds1 << "\n"
+      //          << "extracted D2psi (s2,s1): " << d2psi_ds1_ds2 << "\n"
+      //          << "extracted D2psi (s2,s2): " << d2psi_ds2_ds2 << "\n"
+      //          << std::endl;
+      //      std::cout << "d2psi_ds1_ds1:            " << d2psi_ds1_ds1 <<
+      //      std::endl; std::cout << "func::d2psi_ds1_ds1(s1,s2): " <<
+      //      func::d2psi_ds1_ds1(s1,s2) << std::endl; std::cout << "diff: " <<
+      //      std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1,s2)) << std::endl;
+      //      std::cout << "d2psi_ds2_ds1:            " << d2psi_ds2_ds1 <<
+      //      std::endl; std::cout << "func::d2psi_ds2_ds1(s1,s2): " <<
+      //      func::d2psi_ds2_ds1(s1,s2) << std::endl; std::cout << "diff: " <<
+      //      std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1,s2)) << std::endl;
+      //      std::cout << "d2psi_ds1_ds2:            " << d2psi_ds1_ds2 <<
+      //      std::endl; std::cout << "func::d2psi_ds1_ds2(s1,s2): " <<
+      //      func::d2psi_ds1_ds2(s1,s2) << std::endl; std::cout << "diff: " <<
+      //      std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1,s2)) << std::endl;
+      //      std::cout << "d2psi_ds2_ds2:            " << d2psi_ds2_ds2 <<
+      //      std::endl; std::cout << "func::d2psi_ds2_ds2(s1,s2): " <<
+      //      func::d2psi_ds2_ds2(s1,s2) << std::endl; std::cout << "diff: " <<
+      //      std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1,s2)) << std::endl;
+      Assert(std::abs(d2psi_ds1_ds1 - func::d2psi_ds1_ds1(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+      Assert(std::abs(d2psi_ds2_ds1 - func::d2psi_ds2_ds1(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+      Assert(std::abs(d2psi_ds1_ds2 - func::d2psi_ds1_ds2(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+      Assert(std::abs(d2psi_ds2_ds2 - func::d2psi_ds2_ds2(s1, s2)) < tol,
+             ExcMessage("No match for second derivative."));
+    }
+}
+
+
+
+template <int dim,
+          typename number_t,
+          enum AD::NumberTypes ad_type_code,
+          typename tape_index_t>
+void
+test_scalar_scalar_coupled(const tape_index_t tape_index_begin,
+                           const tape_index_t tape_index_end)
+{
+  for (tape_index_t tape_no = tape_index_begin; tape_no < tape_index_end;
+       ++tape_no)
+    {
+      const number_t s1 = 3.1 * (1.0 + tape_no / 10.0);
+      const number_t s2 = 5.9 * (1.0 + tape_no / 10.0);
+      test_scalar_scalar_coupled_one_run<3, number_t, ad_type_code>(s1,
+                                                                    s2,
+                                                                    tape_no);
+    }
+}
+
+
+int
+main()
+{
+  initlog();
+
+  using number_t                     = double;
+  const AD::NumberTypes ad_type_code = AD::NumberTypes::adolc_taped;
+  using ADNumberType = typename AD::HelperBase<ad_type_code, number_t>::ad_type;
+  using tape_index_t = typename AD::Types<ADNumberType>::tape_index;
+
+  const tape_index_t invalid_tape_index =
+    AD::Numbers<ADNumberType>::invalid_tape_index;
+  const tape_index_t max_tape_index = AD::Numbers<ADNumberType>::max_tape_index;
+
+  deallog.push("Valid range");
+  {
+    test_scalar_scalar_coupled<2, number_t, ad_type_code>(
+      static_cast<tape_index_t>(invalid_tape_index + 1), max_tape_index);
+    test_scalar_scalar_coupled<3, number_t, ad_type_code>(
+      static_cast<tape_index_t>(invalid_tape_index + 1), max_tape_index);
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  // We expect the following two checks to fail
+  deal_II_exceptions::disable_abort_on_exception();
+
+  deallog.push("Invalid index");
+  {
+    try
+      {
+        test_scalar_scalar_coupled<2, number_t, ad_type_code>(
+          invalid_tape_index,
+          static_cast<tape_index_t>(invalid_tape_index + 1));
+        test_scalar_scalar_coupled<3, number_t, ad_type_code>(
+          invalid_tape_index,
+          static_cast<tape_index_t>(invalid_tape_index + 1));
+        deallog << "Expected failure" << std::endl;
+      }
+    catch (...)
+      {
+        deallog << "OK" << std::endl;
+      }
+  }
+  deallog.pop();
+
+  deallog.push("Max index");
+  {
+    try
+      {
+        test_scalar_scalar_coupled<2, number_t, ad_type_code>(
+          max_tape_index, static_cast<tape_index_t>(max_tape_index + 1));
+        test_scalar_scalar_coupled<3, number_t, ad_type_code>(
+          max_tape_index, static_cast<tape_index_t>(max_tape_index + 1));
+        deallog << "Expected failure" << std::endl;
+      }
+    catch (...)
+      {
+        deallog << "OK" << std::endl;
+      }
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_tape_index_01.output b/tests/adolc/helper_tape_index_01.output
new file mode 100644 (file)
index 0000000..666b4bc
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL:Valid range::OK
+DEAL:Invalid index::OK
+DEAL:Max index::OK
+DEAL::OK
diff --git a/tests/adolc/helper_vector_1_indep_1_dep_vars_01_1.cc b/tests/adolc/helper_vector_1_indep_1_dep_vars_01_1.cc
new file mode 100644 (file)
index 0000000..687ea6b
--- /dev/null
@@ -0,0 +1,40 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a single component (scalar) system using a helper class
+// This is the equivalent of helper_scalar_single_component_01, but for
+// the other (vector) helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_vector_1_indep_1_dep_vars_01.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_taped>();
+    test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_vector_1_indep_1_dep_vars_01_1.output b/tests/adolc/helper_vector_1_indep_1_dep_vars_01_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_vector_1_indep_1_dep_vars_01_2.cc b/tests/adolc/helper_vector_1_indep_1_dep_vars_01_2.cc
new file mode 100644 (file)
index 0000000..525b2c7
--- /dev/null
@@ -0,0 +1,40 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a single component (scalar) system using a helper class
+// This is the equivalent of helper_scalar_single_component_01, but for
+// the other (vector) helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_vector_1_indep_1_dep_vars_01.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_vector_1_indep_1_dep_vars_01_2.output b/tests/adolc/helper_vector_1_indep_1_dep_vars_01_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_vector_1_indep_2_dep_vars_01_1.cc b/tests/adolc/helper_vector_1_indep_2_dep_vars_01_1.cc
new file mode 100644 (file)
index 0000000..b63f18e
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a vector of 2 dependent and 1 independent variables
+// using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_vector_1_indep_2_dep_vars_01.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_taped>();
+    test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_vector_1_indep_2_dep_vars_01_1.output b/tests/adolc/helper_vector_1_indep_2_dep_vars_01_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_vector_1_indep_2_dep_vars_01_2.cc b/tests/adolc/helper_vector_1_indep_2_dep_vars_01_2.cc
new file mode 100644 (file)
index 0000000..da14614
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a vector of 2 dependent and 1 independent variables
+// using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_vector_1_indep_2_dep_vars_01.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_vector_1_indep_2_dep_vars_01_2.output b/tests/adolc/helper_vector_1_indep_2_dep_vars_01_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_vector_2_indep_1_dep_vars_01_1.cc b/tests/adolc/helper_vector_2_indep_1_dep_vars_01_1.cc
new file mode 100644 (file)
index 0000000..a725273
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a vector of 1 dependent and 2 independent variables
+// using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_vector_2_indep_1_dep_vars_01.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_taped>();
+    test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_vector_2_indep_1_dep_vars_01_1.output b/tests/adolc/helper_vector_2_indep_1_dep_vars_01_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_vector_2_indep_1_dep_vars_01_2.cc b/tests/adolc/helper_vector_2_indep_1_dep_vars_01_2.cc
new file mode 100644 (file)
index 0000000..3371380
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a vector of 1 dependent and 2 independent variables
+// using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_vector_2_indep_1_dep_vars_01.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_vector_2_indep_1_dep_vars_01_2.output b/tests/adolc/helper_vector_2_indep_1_dep_vars_01_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_vector_2_indep_2_dep_vars_01_1.cc b/tests/adolc/helper_vector_2_indep_2_dep_vars_01_1.cc
new file mode 100644 (file)
index 0000000..264f6c2
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a vector of 2 dependent and 2 independent variables
+// using a helper class
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/helper_vector_2_indep_2_dep_vars_01.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_taped>();
+    test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_vector_2_indep_2_dep_vars_01_1.output b/tests/adolc/helper_vector_2_indep_2_dep_vars_01_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/helper_vector_2_indep_2_dep_vars_01_2.cc b/tests/adolc/helper_vector_2_indep_2_dep_vars_01_2.cc
new file mode 100644 (file)
index 0000000..969a2fa
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluation of a vector of 2 dependent and 2 independent variables
+// using a helper class
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/helper_vector_2_indep_2_dep_vars_01.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_AD_vector_jacobian<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_AD_vector_jacobian<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/helper_vector_2_indep_2_dep_vars_01_2.output b/tests/adolc/helper_vector_2_indep_2_dep_vars_01_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/step-44-helper_scalar_01_1.cc b/tests/adolc/step-44-helper_scalar_01_1.cc
new file mode 100644 (file)
index 0000000..b1de2d1
--- /dev/null
@@ -0,0 +1,68 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// This is a modified version of step-44, which tests the implementation of
+// QP-level auto-differentiation.
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/step-44-helper_scalar_01.h"
+#include "../tests.h"
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+
+  using namespace dealii;
+  namespace AD = dealii::Differentiation::AD;
+  using namespace Step44;
+  try
+    {
+      const unsigned int                               dim = 2;
+      Solid<dim, double, AD::NumberTypes::adolc_taped> solid(
+        SOURCE_DIR "/prm/parameters-step-44.prm");
+      solid.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  return 0;
+}
diff --git a/tests/adolc/step-44-helper_scalar_01_1.output b/tests/adolc/step-44-helper_scalar_01_1.output
new file mode 100644 (file)
index 0000000..8ba79ae
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::Timestep 1: 0.00 -0.000119
+DEAL::Timestep 2: 0.00 -0.000228
+DEAL::Timestep 3: 0.00 -0.000322
+DEAL::Timestep 4: 0.00 -0.000402
+DEAL::Timestep 5: 0.00 -0.000470
+DEAL::Timestep 6: 0.00 -0.000527
+DEAL::Timestep 7: 0.00 -0.000576
+DEAL::Timestep 8: 0.00 -0.000617
+DEAL::Timestep 9: 0.00 -0.000653
+DEAL::Timestep 10: 0.00 -0.000683
diff --git a/tests/adolc/step-44-helper_vector_01_1.cc b/tests/adolc/step-44-helper_vector_01_1.cc
new file mode 100644 (file)
index 0000000..2758838
--- /dev/null
@@ -0,0 +1,68 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// This is a modified version of step-44, which tests the implementation of
+// cell-level auto-differentiation (linearisation of a residual vector).
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/step-44-helper_vector_01.h"
+#include "../tests.h"
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+
+  using namespace dealii;
+  namespace AD = dealii::Differentiation::AD;
+  using namespace Step44;
+  try
+    {
+      const unsigned int                               dim = 2;
+      Solid<dim, double, AD::NumberTypes::adolc_taped> solid(
+        SOURCE_DIR "/prm/parameters-step-44.prm");
+      solid.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  return 0;
+}
diff --git a/tests/adolc/step-44-helper_vector_01_1.output b/tests/adolc/step-44-helper_vector_01_1.output
new file mode 100644 (file)
index 0000000..8ba79ae
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::Timestep 1: 0.00 -0.000119
+DEAL::Timestep 2: 0.00 -0.000228
+DEAL::Timestep 3: 0.00 -0.000322
+DEAL::Timestep 4: 0.00 -0.000402
+DEAL::Timestep 5: 0.00 -0.000470
+DEAL::Timestep 6: 0.00 -0.000527
+DEAL::Timestep 7: 0.00 -0.000576
+DEAL::Timestep 8: 0.00 -0.000617
+DEAL::Timestep 9: 0.00 -0.000653
+DEAL::Timestep 10: 0.00 -0.000683
diff --git a/tests/adolc/step-44-helper_vector_01_2.cc b/tests/adolc/step-44-helper_vector_01_2.cc
new file mode 100644 (file)
index 0000000..55f51f4
--- /dev/null
@@ -0,0 +1,68 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// This is a modified version of step-44, which tests the implementation of
+// cell-level auto-differentiation (linearisation of a residual vector).
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/step-44-helper_vector_01.h"
+#include "../tests.h"
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+
+  using namespace dealii;
+  namespace AD = dealii::Differentiation::AD;
+  using namespace Step44;
+  try
+    {
+      const unsigned int                                  dim = 2;
+      Solid<dim, double, AD::NumberTypes::adolc_tapeless> solid(
+        SOURCE_DIR "/prm/parameters-step-44.prm");
+      solid.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  return 0;
+}
diff --git a/tests/adolc/step-44-helper_vector_01_2.output b/tests/adolc/step-44-helper_vector_01_2.output
new file mode 100644 (file)
index 0000000..8ba79ae
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::Timestep 1: 0.00 -0.000119
+DEAL::Timestep 2: 0.00 -0.000228
+DEAL::Timestep 3: 0.00 -0.000322
+DEAL::Timestep 4: 0.00 -0.000402
+DEAL::Timestep 5: 0.00 -0.000470
+DEAL::Timestep 6: 0.00 -0.000527
+DEAL::Timestep 7: 0.00 -0.000576
+DEAL::Timestep 8: 0.00 -0.000617
+DEAL::Timestep 9: 0.00 -0.000653
+DEAL::Timestep 10: 0.00 -0.000683
diff --git a/tests/adolc/step-44.cc b/tests/adolc/step-44.cc
new file mode 100644 (file)
index 0000000..7e7c78f
--- /dev/null
@@ -0,0 +1,2128 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// This is a copy of step-44 (git rev 3f7e617) to use as a base-line for
+// results produced via different approaches to be compared to.
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/quadrature_point_data.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/base/work_stream.h>
+
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgp_monomial.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q_eulerian.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/packaged_operation.h>
+#include <deal.II/lac/precondition_selector.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_selector.h>
+#include <deal.II/lac/sparse_direct.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <functional>
+#include <iostream>
+
+#include "../tests.h"
+namespace Step44
+{
+  using namespace dealii;
+  namespace Parameters
+  {
+    struct FESystem
+    {
+      unsigned int poly_degree;
+      unsigned int quad_order;
+      static void
+      declare_parameters(ParameterHandler &prm);
+      void
+      parse_parameters(ParameterHandler &prm);
+    };
+    void
+    FESystem::declare_parameters(ParameterHandler &prm)
+    {
+      prm.enter_subsection("Finite element system");
+      {
+        prm.declare_entry("Polynomial degree",
+                          "2",
+                          Patterns::Integer(0),
+                          "Displacement system polynomial order");
+        prm.declare_entry("Quadrature order",
+                          "3",
+                          Patterns::Integer(0),
+                          "Gauss quadrature order");
+      }
+      prm.leave_subsection();
+    }
+    void
+    FESystem::parse_parameters(ParameterHandler &prm)
+    {
+      prm.enter_subsection("Finite element system");
+      {
+        poly_degree = prm.get_integer("Polynomial degree");
+        quad_order  = prm.get_integer("Quadrature order");
+      }
+      prm.leave_subsection();
+    }
+    struct Geometry
+    {
+      unsigned int global_refinement;
+      double       scale;
+      double       p_p0;
+      static void
+      declare_parameters(ParameterHandler &prm);
+      void
+      parse_parameters(ParameterHandler &prm);
+    };
+    void
+    Geometry::declare_parameters(ParameterHandler &prm)
+    {
+      prm.enter_subsection("Geometry");
+      {
+        prm.declare_entry("Global refinement",
+                          "2",
+                          Patterns::Integer(0),
+                          "Global refinement level");
+        prm.declare_entry("Grid scale",
+                          "1e-3",
+                          Patterns::Double(0.0),
+                          "Global grid scaling factor");
+        prm.declare_entry("Pressure ratio p/p0",
+                          "100",
+                          Patterns::Selection("20|40|60|80|100"),
+                          "Ratio of applied pressure to reference pressure");
+      }
+      prm.leave_subsection();
+    }
+    void
+    Geometry::parse_parameters(ParameterHandler &prm)
+    {
+      prm.enter_subsection("Geometry");
+      {
+        global_refinement = prm.get_integer("Global refinement");
+        scale             = prm.get_double("Grid scale");
+        p_p0              = prm.get_double("Pressure ratio p/p0");
+      }
+      prm.leave_subsection();
+    }
+    struct Materials
+    {
+      double nu;
+      double mu;
+      static void
+      declare_parameters(ParameterHandler &prm);
+      void
+      parse_parameters(ParameterHandler &prm);
+    };
+    void
+    Materials::declare_parameters(ParameterHandler &prm)
+    {
+      prm.enter_subsection("Material properties");
+      {
+        prm.declare_entry("Poisson's ratio",
+                          "0.4999",
+                          Patterns::Double(-1.0, 0.5),
+                          "Poisson's ratio");
+        prm.declare_entry("Shear modulus",
+                          "80.194e6",
+                          Patterns::Double(),
+                          "Shear modulus");
+      }
+      prm.leave_subsection();
+    }
+    void
+    Materials::parse_parameters(ParameterHandler &prm)
+    {
+      prm.enter_subsection("Material properties");
+      {
+        nu = prm.get_double("Poisson's ratio");
+        mu = prm.get_double("Shear modulus");
+      }
+      prm.leave_subsection();
+    }
+    struct LinearSolver
+    {
+      std::string type_lin;
+      double      tol_lin;
+      double      max_iterations_lin;
+      bool        use_static_condensation;
+      std::string preconditioner_type;
+      double      preconditioner_relaxation;
+      static void
+      declare_parameters(ParameterHandler &prm);
+      void
+      parse_parameters(ParameterHandler &prm);
+    };
+    void
+    LinearSolver::declare_parameters(ParameterHandler &prm)
+    {
+      prm.enter_subsection("Linear solver");
+      {
+        prm.declare_entry("Solver type",
+                          "CG",
+                          Patterns::Selection("CG|Direct"),
+                          "Type of solver used to solve the linear system");
+        prm.declare_entry("Residual",
+                          "1e-6",
+                          Patterns::Double(0.0),
+                          "Linear solver residual (scaled by residual norm)");
+        prm.declare_entry(
+          "Max iteration multiplier",
+          "1",
+          Patterns::Double(0.0),
+          "Linear solver iterations (multiples of the system matrix size)");
+        prm.declare_entry("Use static condensation",
+                          "true",
+                          Patterns::Bool(),
+                          "Solve the full block system or a reduced problem");
+        prm.declare_entry("Preconditioner type",
+                          "ssor",
+                          Patterns::Selection("jacobi|ssor"),
+                          "Type of preconditioner");
+        prm.declare_entry("Preconditioner relaxation",
+                          "0.65",
+                          Patterns::Double(0.0),
+                          "Preconditioner relaxation value");
+      }
+      prm.leave_subsection();
+    }
+    void
+    LinearSolver::parse_parameters(ParameterHandler &prm)
+    {
+      prm.enter_subsection("Linear solver");
+      {
+        type_lin                  = prm.get("Solver type");
+        tol_lin                   = prm.get_double("Residual");
+        max_iterations_lin        = prm.get_double("Max iteration multiplier");
+        use_static_condensation   = prm.get_bool("Use static condensation");
+        preconditioner_type       = prm.get("Preconditioner type");
+        preconditioner_relaxation = prm.get_double("Preconditioner relaxation");
+      }
+      prm.leave_subsection();
+    }
+    struct NonlinearSolver
+    {
+      unsigned int max_iterations_NR;
+      double       tol_f;
+      double       tol_u;
+      static void
+      declare_parameters(ParameterHandler &prm);
+      void
+      parse_parameters(ParameterHandler &prm);
+    };
+    void
+    NonlinearSolver::declare_parameters(ParameterHandler &prm)
+    {
+      prm.enter_subsection("Nonlinear solver");
+      {
+        prm.declare_entry("Max iterations Newton-Raphson",
+                          "10",
+                          Patterns::Integer(0),
+                          "Number of Newton-Raphson iterations allowed");
+        prm.declare_entry("Tolerance force",
+                          "1.0e-9",
+                          Patterns::Double(0.0),
+                          "Force residual tolerance");
+        prm.declare_entry("Tolerance displacement",
+                          "1.0e-6",
+                          Patterns::Double(0.0),
+                          "Displacement error tolerance");
+      }
+      prm.leave_subsection();
+    }
+    void
+    NonlinearSolver::parse_parameters(ParameterHandler &prm)
+    {
+      prm.enter_subsection("Nonlinear solver");
+      {
+        max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson");
+        tol_f             = prm.get_double("Tolerance force");
+        tol_u             = prm.get_double("Tolerance displacement");
+      }
+      prm.leave_subsection();
+    }
+    struct Time
+    {
+      double delta_t;
+      double end_time;
+      static void
+      declare_parameters(ParameterHandler &prm);
+      void
+      parse_parameters(ParameterHandler &prm);
+    };
+    void
+    Time::declare_parameters(ParameterHandler &prm)
+    {
+      prm.enter_subsection("Time");
+      {
+        prm.declare_entry("End time", "1", Patterns::Double(), "End time");
+        prm.declare_entry("Time step size",
+                          "0.1",
+                          Patterns::Double(),
+                          "Time step size");
+      }
+      prm.leave_subsection();
+    }
+    void
+    Time::parse_parameters(ParameterHandler &prm)
+    {
+      prm.enter_subsection("Time");
+      {
+        end_time = prm.get_double("End time");
+        delta_t  = prm.get_double("Time step size");
+      }
+      prm.leave_subsection();
+    }
+    struct AllParameters : public FESystem,
+                           public Geometry,
+                           public Materials,
+                           public LinearSolver,
+                           public NonlinearSolver,
+                           public Time
+    {
+      AllParameters(const std::string &input_file);
+      static void
+      declare_parameters(ParameterHandler &prm);
+      void
+      parse_parameters(ParameterHandler &prm);
+    };
+    AllParameters::AllParameters(const std::string &input_file)
+    {
+      ParameterHandler prm;
+      declare_parameters(prm);
+      prm.parse_input(input_file);
+      parse_parameters(prm);
+    }
+    void
+    AllParameters::declare_parameters(ParameterHandler &prm)
+    {
+      FESystem::declare_parameters(prm);
+      Geometry::declare_parameters(prm);
+      Materials::declare_parameters(prm);
+      LinearSolver::declare_parameters(prm);
+      NonlinearSolver::declare_parameters(prm);
+      Time::declare_parameters(prm);
+    }
+    void
+    AllParameters::parse_parameters(ParameterHandler &prm)
+    {
+      FESystem::parse_parameters(prm);
+      Geometry::parse_parameters(prm);
+      Materials::parse_parameters(prm);
+      LinearSolver::parse_parameters(prm);
+      NonlinearSolver::parse_parameters(prm);
+      Time::parse_parameters(prm);
+    }
+  } // namespace Parameters
+  template <int dim>
+  class StandardTensors
+  {
+  public:
+    static const SymmetricTensor<2, dim> I;
+    static const SymmetricTensor<4, dim> IxI;
+    static const SymmetricTensor<4, dim> II;
+    static const SymmetricTensor<4, dim> dev_P;
+  };
+  template <int dim>
+  const SymmetricTensor<2, dim>
+    StandardTensors<dim>::I = unit_symmetric_tensor<dim>();
+  template <int dim>
+  const SymmetricTensor<4, dim> StandardTensors<dim>::IxI = outer_product(I, I);
+  template <int dim>
+  const SymmetricTensor<4, dim>
+    StandardTensors<dim>::II = identity_tensor<dim>();
+  template <int dim>
+  const SymmetricTensor<4, dim>
+    StandardTensors<dim>::dev_P = deviator_tensor<dim>();
+  class Time
+  {
+  public:
+    Time(const double time_end, const double delta_t)
+      : timestep(0)
+      , time_current(0.0)
+      , time_end(time_end)
+      , delta_t(delta_t)
+    {}
+    virtual ~Time()
+    {}
+    double
+    current() const
+    {
+      return time_current;
+    }
+    double
+    end() const
+    {
+      return time_end;
+    }
+    double
+    get_delta_t() const
+    {
+      return delta_t;
+    }
+    unsigned int
+    get_timestep() const
+    {
+      return timestep;
+    }
+    void
+    increment()
+    {
+      time_current += delta_t;
+      ++timestep;
+    }
+
+  private:
+    unsigned int timestep;
+    double       time_current;
+    const double time_end;
+    const double delta_t;
+  };
+  template <int dim>
+  class Material_Compressible_Neo_Hook_Three_Field
+  {
+  public:
+    Material_Compressible_Neo_Hook_Three_Field(const double mu, const double nu)
+      : kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu)))
+      , c_1(mu / 2.0)
+      , det_F(1.0)
+      , p_tilde(0.0)
+      , J_tilde(1.0)
+      , b_bar(StandardTensors<dim>::I)
+    {
+      Assert(kappa > 0, ExcInternalError());
+    }
+    ~Material_Compressible_Neo_Hook_Three_Field()
+    {}
+    void
+    update_material_data(const Tensor<2, dim> &F,
+                         const double          p_tilde_in,
+                         const double          J_tilde_in)
+    {
+      det_F   = determinant(F);
+      b_bar   = std::pow(det_F, -2.0 / dim) * symmetrize(F * transpose(F));
+      p_tilde = p_tilde_in;
+      J_tilde = J_tilde_in;
+      Assert(det_F > 0, ExcInternalError());
+    }
+    SymmetricTensor<2, dim>
+    get_tau()
+    {
+      return get_tau_iso() + get_tau_vol();
+    }
+    SymmetricTensor<4, dim>
+    get_Jc() const
+    {
+      return get_Jc_vol() + get_Jc_iso();
+    }
+    double
+    get_dPsi_vol_dJ() const
+    {
+      return (kappa / 2.0) * (J_tilde - 1.0 / J_tilde);
+    }
+    double
+    get_d2Psi_vol_dJ2() const
+    {
+      return ((kappa / 2.0) * (1.0 + 1.0 / (J_tilde * J_tilde)));
+    }
+    double
+    get_det_F() const
+    {
+      return det_F;
+    }
+    double
+    get_p_tilde() const
+    {
+      return p_tilde;
+    }
+    double
+    get_J_tilde() const
+    {
+      return J_tilde;
+    }
+
+  protected:
+    const double            kappa;
+    const double            c_1;
+    double                  det_F;
+    double                  p_tilde;
+    double                  J_tilde;
+    SymmetricTensor<2, dim> b_bar;
+    SymmetricTensor<2, dim>
+    get_tau_vol() const
+    {
+      return p_tilde * det_F * StandardTensors<dim>::I;
+    }
+    SymmetricTensor<2, dim>
+    get_tau_iso() const
+    {
+      return StandardTensors<dim>::dev_P * get_tau_bar();
+    }
+    SymmetricTensor<2, dim>
+    get_tau_bar() const
+    {
+      return 2.0 * c_1 * b_bar;
+    }
+    SymmetricTensor<4, dim>
+    get_Jc_vol() const
+    {
+      return p_tilde * det_F *
+             (StandardTensors<dim>::IxI - (2.0 * StandardTensors<dim>::II));
+    }
+    SymmetricTensor<4, dim>
+    get_Jc_iso() const
+    {
+      const SymmetricTensor<2, dim> tau_bar = get_tau_bar();
+      const SymmetricTensor<2, dim> tau_iso = get_tau_iso();
+      const SymmetricTensor<4, dim> tau_iso_x_I =
+        outer_product(tau_iso, StandardTensors<dim>::I);
+      const SymmetricTensor<4, dim> I_x_tau_iso =
+        outer_product(StandardTensors<dim>::I, tau_iso);
+      const SymmetricTensor<4, dim> c_bar = get_c_bar();
+      return (2.0 / dim) * trace(tau_bar) * StandardTensors<dim>::dev_P -
+             (2.0 / dim) * (tau_iso_x_I + I_x_tau_iso) +
+             StandardTensors<dim>::dev_P * c_bar * StandardTensors<dim>::dev_P;
+    }
+    SymmetricTensor<4, dim>
+    get_c_bar() const
+    {
+      return SymmetricTensor<4, dim>();
+    }
+  };
+  template <int dim>
+  class PointHistory
+  {
+  public:
+    PointHistory()
+      : F_inv(StandardTensors<dim>::I)
+      , tau(SymmetricTensor<2, dim>())
+      , d2Psi_vol_dJ2(0.0)
+      , dPsi_vol_dJ(0.0)
+      , Jc(SymmetricTensor<4, dim>())
+    {}
+    virtual ~PointHistory()
+    {}
+    void
+    setup_lqp(const Parameters::AllParameters &parameters)
+    {
+      material.reset(
+        new Material_Compressible_Neo_Hook_Three_Field<dim>(parameters.mu,
+                                                            parameters.nu));
+      update_values(Tensor<2, dim>(), 0.0, 1.0);
+    }
+    void
+    update_values(const Tensor<2, dim> &Grad_u_n,
+                  const double          p_tilde,
+                  const double          J_tilde)
+    {
+      const Tensor<2, dim> F =
+        (Tensor<2, dim>(StandardTensors<dim>::I) + Grad_u_n);
+      material->update_material_data(F, p_tilde, J_tilde);
+      F_inv         = invert(F);
+      tau           = material->get_tau();
+      Jc            = material->get_Jc();
+      dPsi_vol_dJ   = material->get_dPsi_vol_dJ();
+      d2Psi_vol_dJ2 = material->get_d2Psi_vol_dJ2();
+    }
+    double
+    get_J_tilde() const
+    {
+      return material->get_J_tilde();
+    }
+    double
+    get_det_F() const
+    {
+      return material->get_det_F();
+    }
+    const Tensor<2, dim> &
+    get_F_inv() const
+    {
+      return F_inv;
+    }
+    double
+    get_p_tilde() const
+    {
+      return material->get_p_tilde();
+    }
+    const SymmetricTensor<2, dim> &
+    get_tau() const
+    {
+      return tau;
+    }
+    double
+    get_dPsi_vol_dJ() const
+    {
+      return dPsi_vol_dJ;
+    }
+    double
+    get_d2Psi_vol_dJ2() const
+    {
+      return d2Psi_vol_dJ2;
+    }
+    const SymmetricTensor<4, dim> &
+    get_Jc() const
+    {
+      return Jc;
+    }
+
+  private:
+    std::shared_ptr<Material_Compressible_Neo_Hook_Three_Field<dim>> material;
+    Tensor<2, dim>                                                   F_inv;
+    SymmetricTensor<2, dim>                                          tau;
+    double                  d2Psi_vol_dJ2;
+    double                  dPsi_vol_dJ;
+    SymmetricTensor<4, dim> Jc;
+  };
+  template <int dim>
+  class Solid
+  {
+  public:
+    Solid(const std::string &input_file);
+    virtual ~Solid();
+    void
+    run();
+
+  private:
+    struct PerTaskData_K;
+    struct ScratchData_K;
+    struct PerTaskData_RHS;
+    struct ScratchData_RHS;
+    struct PerTaskData_SC;
+    struct ScratchData_SC;
+    struct PerTaskData_UQPH;
+    struct ScratchData_UQPH;
+    void
+    make_grid();
+    void
+    system_setup();
+    void
+    determine_component_extractors();
+    void
+    assemble_system_tangent();
+    void
+    assemble_system_tangent_one_cell(
+      const typename DoFHandler<dim>::active_cell_iterator &cell,
+      ScratchData_K &                                       scratch,
+      PerTaskData_K &                                       data) const;
+    void
+    copy_local_to_global_K(const PerTaskData_K &data);
+    void
+    assemble_system_rhs();
+    void
+    assemble_system_rhs_one_cell(
+      const typename DoFHandler<dim>::active_cell_iterator &cell,
+      ScratchData_RHS &                                     scratch,
+      PerTaskData_RHS &                                     data) const;
+    void
+    copy_local_to_global_rhs(const PerTaskData_RHS &data);
+    void
+    assemble_sc();
+    void
+    assemble_sc_one_cell(
+      const typename DoFHandler<dim>::active_cell_iterator &cell,
+      ScratchData_SC &                                      scratch,
+      PerTaskData_SC &                                      data);
+    void
+    copy_local_to_global_sc(const PerTaskData_SC &data);
+    void
+    make_constraints(const int &it_nr);
+    void
+    setup_qph();
+    void
+    update_qph_incremental(const BlockVector<double> &solution_delta);
+    void
+    update_qph_incremental_one_cell(
+      const typename DoFHandler<dim>::active_cell_iterator &cell,
+      ScratchData_UQPH &                                    scratch,
+      PerTaskData_UQPH &                                    data);
+    void
+    copy_local_to_global_UQPH(const PerTaskData_UQPH & /*data*/)
+    {}
+    void
+    solve_nonlinear_timestep(BlockVector<double> &solution_delta);
+    std::pair<unsigned int, double>
+    solve_linear_system(BlockVector<double> &newton_update);
+    BlockVector<double>
+    get_total_solution(const BlockVector<double> &solution_delta) const;
+    void
+                              output_results() const;
+    Parameters::AllParameters parameters;
+    double                    vol_reference;
+    Triangulation<dim>        triangulation;
+    Time                      time;
+    mutable TimerOutput       timer;
+    CellDataStorage<typename Triangulation<dim>::cell_iterator,
+                    PointHistory<dim>>
+                                     quadrature_point_history;
+    const unsigned int               degree;
+    const FESystem<dim>              fe;
+    DoFHandler<dim>                  dof_handler_ref;
+    const unsigned int               dofs_per_cell;
+    const FEValuesExtractors::Vector u_fe;
+    const FEValuesExtractors::Scalar p_fe;
+    const FEValuesExtractors::Scalar J_fe;
+    static const unsigned int        n_blocks          = 3;
+    static const unsigned int        n_components      = dim + 2;
+    static const unsigned int        first_u_component = 0;
+    static const unsigned int        p_component       = dim;
+    static const unsigned int        J_component       = dim + 1;
+    enum
+    {
+      u_dof = 0,
+      p_dof = 1,
+      J_dof = 2
+    };
+    std::vector<types::global_dof_index> dofs_per_block;
+    std::vector<types::global_dof_index> element_indices_u;
+    std::vector<types::global_dof_index> element_indices_p;
+    std::vector<types::global_dof_index> element_indices_J;
+    const QGauss<dim>                    qf_cell;
+    const QGauss<dim - 1>                qf_face;
+    const unsigned int                   n_q_points;
+    const unsigned int                   n_q_points_f;
+    ConstraintMatrix                     constraints;
+    BlockSparsityPattern                 sparsity_pattern;
+    BlockSparseMatrix<double>            tangent_matrix;
+    BlockVector<double>                  system_rhs;
+    BlockVector<double>                  solution_n;
+    struct Errors
+    {
+      Errors()
+        : norm(1.0)
+        , u(1.0)
+        , p(1.0)
+        , J(1.0)
+      {}
+      void
+      reset()
+      {
+        norm = 1.0;
+        u    = 1.0;
+        p    = 1.0;
+        J    = 1.0;
+      }
+      void
+      normalise(const Errors &rhs)
+      {
+        if (rhs.norm != 0.0)
+          norm /= rhs.norm;
+        if (rhs.u != 0.0)
+          u /= rhs.u;
+        if (rhs.p != 0.0)
+          p /= rhs.p;
+        if (rhs.J != 0.0)
+          J /= rhs.J;
+      }
+      double norm, u, p, J;
+    };
+    Errors error_residual, error_residual_0, error_residual_norm, error_update,
+      error_update_0, error_update_norm;
+    void
+    get_error_residual(Errors &error_residual);
+    void
+    get_error_update(const BlockVector<double> &newton_update,
+                     Errors &                   error_update);
+    std::pair<double, double>
+    get_error_dilation() const;
+    double
+    compute_vol_current() const;
+    static void
+    print_conv_header();
+    void
+    print_conv_footer();
+  };
+  template <int dim>
+  Solid<dim>::Solid(const std::string &input_file)
+    : parameters(input_file)
+    , triangulation(Triangulation<dim>::maximum_smoothing)
+    , time(parameters.end_time, parameters.delta_t)
+    , timer(std::cout, TimerOutput::never, TimerOutput::wall_times)
+    , degree(parameters.poly_degree)
+    , fe(FE_Q<dim>(parameters.poly_degree),
+         dim, // displacement
+         FE_DGPMonomial<dim>(parameters.poly_degree - 1),
+         1, // pressure
+         FE_DGPMonomial<dim>(parameters.poly_degree - 1),
+         1)
+    , // dilatation
+    dof_handler_ref(triangulation)
+    , dofs_per_cell(fe.dofs_per_cell)
+    , u_fe(first_u_component)
+    , p_fe(p_component)
+    , J_fe(J_component)
+    , dofs_per_block(n_blocks)
+    , qf_cell(parameters.quad_order)
+    , qf_face(parameters.quad_order)
+    , n_q_points(qf_cell.size())
+    , n_q_points_f(qf_face.size())
+  {
+    Assert(dim == 2 || dim == 3,
+           ExcMessage("This problem only works in 2 or 3 space dimensions."));
+    determine_component_extractors();
+  }
+  template <int dim>
+  Solid<dim>::~Solid()
+  {
+    dof_handler_ref.clear();
+  }
+  template <int dim>
+  void
+  Solid<dim>::run()
+  {
+    make_grid();
+    system_setup();
+    {
+      ConstraintMatrix constraints;
+      constraints.close();
+      const ComponentSelectFunction<dim> J_mask(J_component, n_components);
+      VectorTools::project(dof_handler_ref,
+                           constraints,
+                           QGauss<dim>(degree + 2),
+                           J_mask,
+                           solution_n);
+    }
+    // output_results();
+    time.increment();
+    BlockVector<double> solution_delta(dofs_per_block);
+    while (time.current() < time.end())
+      {
+        solution_delta = 0.0;
+        solve_nonlinear_timestep(solution_delta);
+        solution_n += solution_delta;
+        // output_results();
+        // Output displacement at centre of traction surface
+        {
+          const Point<dim> soln_pt(
+            dim == 3 ? Point<dim>(0.0, 1.0, 0.0) * parameters.scale :
+                       Point<dim>(0.0, 1.0) * parameters.scale);
+          typename DoFHandler<dim>::active_cell_iterator
+            cell = dof_handler_ref.begin_active(),
+            endc = dof_handler_ref.end();
+          for (; cell != endc; ++cell)
+            for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
+                 ++v)
+              if (cell->vertex(v).distance(soln_pt) < 1e-6 * parameters.scale)
+                {
+                  Tensor<1, dim> soln;
+                  for (unsigned int d = 0; d < dim; ++d)
+                    soln[d] = solution_n(cell->vertex_dof_index(v, u_dof + d));
+                  deallog << "Timestep " << time.get_timestep() << ": " << soln
+                          << std::endl;
+                }
+        }
+        time.increment();
+      }
+  }
+  template <int dim>
+  struct Solid<dim>::PerTaskData_K
+  {
+    FullMatrix<double>                   cell_matrix;
+    std::vector<types::global_dof_index> local_dof_indices;
+    PerTaskData_K(const unsigned int dofs_per_cell)
+      : cell_matrix(dofs_per_cell, dofs_per_cell)
+      , local_dof_indices(dofs_per_cell)
+    {}
+    void
+    reset()
+    {
+      cell_matrix = 0.0;
+    }
+  };
+  template <int dim>
+  struct Solid<dim>::ScratchData_K
+  {
+    FEValues<dim>                                     fe_values_ref;
+    std::vector<std::vector<double>>                  Nx;
+    std::vector<std::vector<Tensor<2, dim>>>          grad_Nx;
+    std::vector<std::vector<SymmetricTensor<2, dim>>> symm_grad_Nx;
+    ScratchData_K(const FiniteElement<dim> &fe_cell,
+                  const QGauss<dim> &       qf_cell,
+                  const UpdateFlags         uf_cell)
+      : fe_values_ref(fe_cell, qf_cell, uf_cell)
+      , Nx(qf_cell.size(), std::vector<double>(fe_cell.dofs_per_cell))
+      , grad_Nx(qf_cell.size(),
+                std::vector<Tensor<2, dim>>(fe_cell.dofs_per_cell))
+      , symm_grad_Nx(qf_cell.size(),
+                     std::vector<SymmetricTensor<2, dim>>(
+                       fe_cell.dofs_per_cell))
+    {}
+    ScratchData_K(const ScratchData_K &rhs)
+      : fe_values_ref(rhs.fe_values_ref.get_fe(),
+                      rhs.fe_values_ref.get_quadrature(),
+                      rhs.fe_values_ref.get_update_flags())
+      , Nx(rhs.Nx)
+      , grad_Nx(rhs.grad_Nx)
+      , symm_grad_Nx(rhs.symm_grad_Nx)
+    {}
+    void
+    reset()
+    {
+      const unsigned int n_q_points      = Nx.size();
+      const unsigned int n_dofs_per_cell = Nx[0].size();
+      for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        {
+          Assert(Nx[q_point].size() == n_dofs_per_cell, ExcInternalError());
+          Assert(grad_Nx[q_point].size() == n_dofs_per_cell,
+                 ExcInternalError());
+          Assert(symm_grad_Nx[q_point].size() == n_dofs_per_cell,
+                 ExcInternalError());
+          for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
+            {
+              Nx[q_point][k]           = 0.0;
+              grad_Nx[q_point][k]      = 0.0;
+              symm_grad_Nx[q_point][k] = 0.0;
+            }
+        }
+    }
+  };
+  template <int dim>
+  struct Solid<dim>::PerTaskData_RHS
+  {
+    Vector<double>                       cell_rhs;
+    std::vector<types::global_dof_index> local_dof_indices;
+    PerTaskData_RHS(const unsigned int dofs_per_cell)
+      : cell_rhs(dofs_per_cell)
+      , local_dof_indices(dofs_per_cell)
+    {}
+    void
+    reset()
+    {
+      cell_rhs = 0.0;
+    }
+  };
+  template <int dim>
+  struct Solid<dim>::ScratchData_RHS
+  {
+    FEValues<dim>                                     fe_values_ref;
+    FEFaceValues<dim>                                 fe_face_values_ref;
+    std::vector<std::vector<double>>                  Nx;
+    std::vector<std::vector<SymmetricTensor<2, dim>>> symm_grad_Nx;
+    ScratchData_RHS(const FiniteElement<dim> &fe_cell,
+                    const QGauss<dim> &       qf_cell,
+                    const UpdateFlags         uf_cell,
+                    const QGauss<dim - 1> &   qf_face,
+                    const UpdateFlags         uf_face)
+      : fe_values_ref(fe_cell, qf_cell, uf_cell)
+      , fe_face_values_ref(fe_cell, qf_face, uf_face)
+      , Nx(qf_cell.size(), std::vector<double>(fe_cell.dofs_per_cell))
+      , symm_grad_Nx(qf_cell.size(),
+                     std::vector<SymmetricTensor<2, dim>>(
+                       fe_cell.dofs_per_cell))
+    {}
+    ScratchData_RHS(const ScratchData_RHS &rhs)
+      : fe_values_ref(rhs.fe_values_ref.get_fe(),
+                      rhs.fe_values_ref.get_quadrature(),
+                      rhs.fe_values_ref.get_update_flags())
+      , fe_face_values_ref(rhs.fe_face_values_ref.get_fe(),
+                           rhs.fe_face_values_ref.get_quadrature(),
+                           rhs.fe_face_values_ref.get_update_flags())
+      , Nx(rhs.Nx)
+      , symm_grad_Nx(rhs.symm_grad_Nx)
+    {}
+    void
+    reset()
+    {
+      const unsigned int n_q_points      = Nx.size();
+      const unsigned int n_dofs_per_cell = Nx[0].size();
+      for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        {
+          Assert(Nx[q_point].size() == n_dofs_per_cell, ExcInternalError());
+          Assert(symm_grad_Nx[q_point].size() == n_dofs_per_cell,
+                 ExcInternalError());
+          for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
+            {
+              Nx[q_point][k]           = 0.0;
+              symm_grad_Nx[q_point][k] = 0.0;
+            }
+        }
+    }
+  };
+  template <int dim>
+  struct Solid<dim>::PerTaskData_SC
+  {
+    FullMatrix<double>                   cell_matrix;
+    std::vector<types::global_dof_index> local_dof_indices;
+    FullMatrix<double>                   k_orig;
+    FullMatrix<double>                   k_pu;
+    FullMatrix<double>                   k_pJ;
+    FullMatrix<double>                   k_JJ;
+    FullMatrix<double>                   k_pJ_inv;
+    FullMatrix<double>                   k_bbar;
+    FullMatrix<double>                   A;
+    FullMatrix<double>                   B;
+    FullMatrix<double>                   C;
+    PerTaskData_SC(const unsigned int dofs_per_cell,
+                   const unsigned int n_u,
+                   const unsigned int n_p,
+                   const unsigned int n_J)
+      : cell_matrix(dofs_per_cell, dofs_per_cell)
+      , local_dof_indices(dofs_per_cell)
+      , k_orig(dofs_per_cell, dofs_per_cell)
+      , k_pu(n_p, n_u)
+      , k_pJ(n_p, n_J)
+      , k_JJ(n_J, n_J)
+      , k_pJ_inv(n_p, n_J)
+      , k_bbar(n_u, n_u)
+      , A(n_J, n_u)
+      , B(n_J, n_u)
+      , C(n_p, n_u)
+    {}
+    void
+    reset()
+    {}
+  };
+  template <int dim>
+  struct Solid<dim>::ScratchData_SC
+  {
+    void
+    reset()
+    {}
+  };
+  template <int dim>
+  struct Solid<dim>::PerTaskData_UQPH
+  {
+    void
+    reset()
+    {}
+  };
+  template <int dim>
+  struct Solid<dim>::ScratchData_UQPH
+  {
+    const BlockVector<double> & solution_total;
+    std::vector<Tensor<2, dim>> solution_grads_u_total;
+    std::vector<double>         solution_values_p_total;
+    std::vector<double>         solution_values_J_total;
+    FEValues<dim>               fe_values_ref;
+    ScratchData_UQPH(const FiniteElement<dim> & fe_cell,
+                     const QGauss<dim> &        qf_cell,
+                     const UpdateFlags          uf_cell,
+                     const BlockVector<double> &solution_total)
+      : solution_total(solution_total)
+      , solution_grads_u_total(qf_cell.size())
+      , solution_values_p_total(qf_cell.size())
+      , solution_values_J_total(qf_cell.size())
+      , fe_values_ref(fe_cell, qf_cell, uf_cell)
+    {}
+    ScratchData_UQPH(const ScratchData_UQPH &rhs)
+      : solution_total(rhs.solution_total)
+      , solution_grads_u_total(rhs.solution_grads_u_total)
+      , solution_values_p_total(rhs.solution_values_p_total)
+      , solution_values_J_total(rhs.solution_values_J_total)
+      , fe_values_ref(rhs.fe_values_ref.get_fe(),
+                      rhs.fe_values_ref.get_quadrature(),
+                      rhs.fe_values_ref.get_update_flags())
+    {}
+    void
+    reset()
+    {
+      const unsigned int n_q_points = solution_grads_u_total.size();
+      for (unsigned int q = 0; q < n_q_points; ++q)
+        {
+          solution_grads_u_total[q]  = 0.0;
+          solution_values_p_total[q] = 0.0;
+          solution_values_J_total[q] = 0.0;
+        }
+    }
+  };
+  template <int dim>
+  void
+  Solid<dim>::make_grid()
+  {
+    GridGenerator::hyper_rectangle(
+      triangulation,
+      (dim == 3 ? Point<dim>(0.0, 0.0, 0.0) : Point<dim>(0.0, 0.0)),
+      (dim == 3 ? Point<dim>(1.0, 1.0, 1.0) : Point<dim>(1.0, 1.0)),
+      true);
+    GridTools::scale(parameters.scale, triangulation);
+    triangulation.refine_global(std::max(1U, parameters.global_refinement));
+    vol_reference = GridTools::volume(triangulation);
+    std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl;
+    typename Triangulation<dim>::active_cell_iterator cell = triangulation
+                                                               .begin_active(),
+                                                      endc =
+                                                        triangulation.end();
+    for (; cell != endc; ++cell)
+      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
+           ++face)
+        {
+          if (cell->face(face)->at_boundary() == true &&
+              cell->face(face)->center()[1] == 1.0 * parameters.scale)
+            {
+              if (dim == 3)
+                {
+                  if (cell->face(face)->center()[0] < 0.5 * parameters.scale &&
+                      cell->face(face)->center()[2] < 0.5 * parameters.scale)
+                    cell->face(face)->set_boundary_id(6);
+                }
+              else
+                {
+                  if (cell->face(face)->center()[0] < 0.5 * parameters.scale)
+                    cell->face(face)->set_boundary_id(6);
+                }
+            }
+        }
+  }
+  template <int dim>
+  void
+  Solid<dim>::system_setup()
+  {
+    timer.enter_subsection("Setup system");
+    std::vector<unsigned int> block_component(n_components,
+                                              u_dof); // Displacement
+    block_component[p_component] = p_dof;             // Pressure
+    block_component[J_component] = J_dof;             // Dilatation
+    dof_handler_ref.distribute_dofs(fe);
+    DoFRenumbering::Cuthill_McKee(dof_handler_ref);
+    DoFRenumbering::component_wise(dof_handler_ref, block_component);
+    DoFTools::count_dofs_per_block(dof_handler_ref,
+                                   dofs_per_block,
+                                   block_component);
+    std::cout << "Triangulation:"
+              << "\n\t Number of active cells: "
+              << triangulation.n_active_cells()
+              << "\n\t Number of degrees of freedom: "
+              << dof_handler_ref.n_dofs() << std::endl;
+    tangent_matrix.clear();
+    {
+      const types::global_dof_index n_dofs_u = dofs_per_block[u_dof];
+      const types::global_dof_index n_dofs_p = dofs_per_block[p_dof];
+      const types::global_dof_index n_dofs_J = dofs_per_block[J_dof];
+      BlockDynamicSparsityPattern   dsp(n_blocks, n_blocks);
+      dsp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u);
+      dsp.block(u_dof, p_dof).reinit(n_dofs_u, n_dofs_p);
+      dsp.block(u_dof, J_dof).reinit(n_dofs_u, n_dofs_J);
+      dsp.block(p_dof, u_dof).reinit(n_dofs_p, n_dofs_u);
+      dsp.block(p_dof, p_dof).reinit(n_dofs_p, n_dofs_p);
+      dsp.block(p_dof, J_dof).reinit(n_dofs_p, n_dofs_J);
+      dsp.block(J_dof, u_dof).reinit(n_dofs_J, n_dofs_u);
+      dsp.block(J_dof, p_dof).reinit(n_dofs_J, n_dofs_p);
+      dsp.block(J_dof, J_dof).reinit(n_dofs_J, n_dofs_J);
+      dsp.collect_sizes();
+      Table<2, DoFTools::Coupling> coupling(n_components, n_components);
+      for (unsigned int ii = 0; ii < n_components; ++ii)
+        for (unsigned int jj = 0; jj < n_components; ++jj)
+          if (((ii < p_component) && (jj == J_component)) ||
+              ((ii == J_component) && (jj < p_component)) ||
+              ((ii == p_component) && (jj == p_component)))
+            coupling[ii][jj] = DoFTools::none;
+          else
+            coupling[ii][jj] = DoFTools::always;
+      DoFTools::make_sparsity_pattern(
+        dof_handler_ref, coupling, dsp, constraints, false);
+      sparsity_pattern.copy_from(dsp);
+    }
+    tangent_matrix.reinit(sparsity_pattern);
+    system_rhs.reinit(dofs_per_block);
+    system_rhs.collect_sizes();
+    solution_n.reinit(dofs_per_block);
+    solution_n.collect_sizes();
+    setup_qph();
+    timer.leave_subsection();
+  }
+  template <int dim>
+  void
+  Solid<dim>::determine_component_extractors()
+  {
+    element_indices_u.clear();
+    element_indices_p.clear();
+    element_indices_J.clear();
+    for (unsigned int k = 0; k < fe.dofs_per_cell; ++k)
+      {
+        const unsigned int k_group = fe.system_to_base_index(k).first.first;
+        if (k_group == u_dof)
+          element_indices_u.push_back(k);
+        else if (k_group == p_dof)
+          element_indices_p.push_back(k);
+        else if (k_group == J_dof)
+          element_indices_J.push_back(k);
+        else
+          {
+            Assert(k_group <= J_dof, ExcInternalError());
+          }
+      }
+  }
+  template <int dim>
+  void
+  Solid<dim>::setup_qph()
+  {
+    std::cout << "    Setting up quadrature point data..." << std::endl;
+    quadrature_point_history.initialize(triangulation.begin_active(),
+                                        triangulation.end(),
+                                        n_q_points);
+    for (typename Triangulation<dim>::active_cell_iterator cell =
+           triangulation.begin_active();
+         cell != triangulation.end();
+         ++cell)
+      {
+        const std::vector<std::shared_ptr<PointHistory<dim>>> lqph =
+          quadrature_point_history.get_data(cell);
+        Assert(lqph.size() == n_q_points, ExcInternalError());
+        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+          lqph[q_point]->setup_lqp(parameters);
+      }
+  }
+  template <int dim>
+  void
+  Solid<dim>::update_qph_incremental(const BlockVector<double> &solution_delta)
+  {
+    timer.enter_subsection("Update QPH data");
+    std::cout << " UQPH " << std::flush;
+    const BlockVector<double> solution_total(
+      get_total_solution(solution_delta));
+    const UpdateFlags uf_UQPH(update_values | update_gradients);
+    PerTaskData_UQPH  per_task_data_UQPH;
+    ScratchData_UQPH  scratch_data_UQPH(fe, qf_cell, uf_UQPH, solution_total);
+    WorkStream::run(dof_handler_ref.begin_active(),
+                    dof_handler_ref.end(),
+                    *this,
+                    &Solid::update_qph_incremental_one_cell,
+                    &Solid::copy_local_to_global_UQPH,
+                    scratch_data_UQPH,
+                    per_task_data_UQPH);
+    timer.leave_subsection();
+  }
+  template <int dim>
+  void
+  Solid<dim>::update_qph_incremental_one_cell(
+    const typename DoFHandler<dim>::active_cell_iterator &cell,
+    ScratchData_UQPH &                                    scratch,
+    PerTaskData_UQPH & /*data*/)
+  {
+    const std::vector<std::shared_ptr<PointHistory<dim>>> lqph =
+      quadrature_point_history.get_data(cell);
+    Assert(lqph.size() == n_q_points, ExcInternalError());
+    Assert(scratch.solution_grads_u_total.size() == n_q_points,
+           ExcInternalError());
+    Assert(scratch.solution_values_p_total.size() == n_q_points,
+           ExcInternalError());
+    Assert(scratch.solution_values_J_total.size() == n_q_points,
+           ExcInternalError());
+    scratch.reset();
+    scratch.fe_values_ref.reinit(cell);
+    scratch.fe_values_ref[u_fe].get_function_gradients(
+      scratch.solution_total, scratch.solution_grads_u_total);
+    scratch.fe_values_ref[p_fe].get_function_values(
+      scratch.solution_total, scratch.solution_values_p_total);
+    scratch.fe_values_ref[J_fe].get_function_values(
+      scratch.solution_total, scratch.solution_values_J_total);
+    for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+      lqph[q_point]->update_values(scratch.solution_grads_u_total[q_point],
+                                   scratch.solution_values_p_total[q_point],
+                                   scratch.solution_values_J_total[q_point]);
+  }
+  template <int dim>
+  void
+  Solid<dim>::solve_nonlinear_timestep(BlockVector<double> &solution_delta)
+  {
+    std::cout << std::endl
+              << "Timestep " << time.get_timestep() << " @ " << time.current()
+              << "s" << std::endl;
+    BlockVector<double> newton_update(dofs_per_block);
+    error_residual.reset();
+    error_residual_0.reset();
+    error_residual_norm.reset();
+    error_update.reset();
+    error_update_0.reset();
+    error_update_norm.reset();
+    print_conv_header();
+    unsigned int newton_iteration = 0;
+    for (; newton_iteration < parameters.max_iterations_NR; ++newton_iteration)
+      {
+        std::cout << " " << std::setw(2) << newton_iteration << " "
+                  << std::flush;
+        tangent_matrix = 0.0;
+        system_rhs     = 0.0;
+        assemble_system_rhs();
+        get_error_residual(error_residual);
+        if (newton_iteration == 0)
+          error_residual_0 = error_residual;
+        error_residual_norm = error_residual;
+        error_residual_norm.normalise(error_residual_0);
+        if (newton_iteration > 0 && error_update_norm.u <= parameters.tol_u &&
+            error_residual_norm.u <= parameters.tol_f)
+          {
+            std::cout << " CONVERGED! " << std::endl;
+            print_conv_footer();
+            break;
+          }
+        assemble_system_tangent();
+        make_constraints(newton_iteration);
+        constraints.condense(tangent_matrix, system_rhs);
+        const std::pair<unsigned int, double> lin_solver_output =
+          solve_linear_system(newton_update);
+        get_error_update(newton_update, error_update);
+        if (newton_iteration == 0)
+          error_update_0 = error_update;
+        error_update_norm = error_update;
+        error_update_norm.normalise(error_update_0);
+        solution_delta += newton_update;
+        update_qph_incremental(solution_delta);
+        std::cout << " | " << std::fixed << std::setprecision(3) << std::setw(7)
+                  << std::scientific << lin_solver_output.first << "  "
+                  << lin_solver_output.second << "  "
+                  << error_residual_norm.norm << "  " << error_residual_norm.u
+                  << "  " << error_residual_norm.p << "  "
+                  << error_residual_norm.J << "  " << error_update_norm.norm
+                  << "  " << error_update_norm.u << "  " << error_update_norm.p
+                  << "  " << error_update_norm.J << "  " << std::endl;
+      }
+    AssertThrow(newton_iteration <= parameters.max_iterations_NR,
+                ExcMessage("No convergence in nonlinear solver!"));
+  }
+  template <int dim>
+  void
+  Solid<dim>::print_conv_header()
+  {
+    static const unsigned int l_width = 155;
+    for (unsigned int i = 0; i < l_width; ++i)
+      std::cout << "_";
+    std::cout << std::endl;
+    std::cout << "                 SOLVER STEP                  "
+              << " |  LIN_IT   LIN_RES    RES_NORM    "
+              << " RES_U     RES_P      RES_J     NU_NORM     "
+              << " NU_U       NU_P       NU_J " << std::endl;
+    for (unsigned int i = 0; i < l_width; ++i)
+      std::cout << "_";
+    std::cout << std::endl;
+  }
+  template <int dim>
+  void
+  Solid<dim>::print_conv_footer()
+  {
+    static const unsigned int l_width = 155;
+    for (unsigned int i = 0; i < l_width; ++i)
+      std::cout << "_";
+    std::cout << std::endl;
+    const std::pair<double, double> error_dil = get_error_dilation();
+    std::cout << "Relative errors:" << std::endl
+              << "Displacement:\t" << error_update.u / error_update_0.u
+              << std::endl
+              << "Force: \t\t" << error_residual.u / error_residual_0.u
+              << std::endl
+              << "Dilatation:\t" << error_dil.first << std::endl
+              << "v / V_0:\t" << error_dil.second * vol_reference << " / "
+              << vol_reference << " = " << error_dil.second << std::endl;
+  }
+  template <int dim>
+  double
+  Solid<dim>::compute_vol_current() const
+  {
+    double        vol_current = 0.0;
+    FEValues<dim> fe_values_ref(fe, qf_cell, update_JxW_values);
+    for (typename Triangulation<dim>::active_cell_iterator cell =
+           triangulation.begin_active();
+         cell != triangulation.end();
+         ++cell)
+      {
+        fe_values_ref.reinit(cell);
+        const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
+          quadrature_point_history.get_data(cell);
+        Assert(lqph.size() == n_q_points, ExcInternalError());
+        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+          {
+            const double det_F_qp = lqph[q_point]->get_det_F();
+            const double JxW      = fe_values_ref.JxW(q_point);
+            vol_current += det_F_qp * JxW;
+          }
+      }
+    Assert(vol_current > 0.0, ExcInternalError());
+    return vol_current;
+  }
+  template <int dim>
+  std::pair<double, double>
+  Solid<dim>::get_error_dilation() const
+  {
+    double        dil_L2_error = 0.0;
+    FEValues<dim> fe_values_ref(fe, qf_cell, update_JxW_values);
+    for (typename Triangulation<dim>::active_cell_iterator cell =
+           triangulation.begin_active();
+         cell != triangulation.end();
+         ++cell)
+      {
+        fe_values_ref.reinit(cell);
+        const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
+          quadrature_point_history.get_data(cell);
+        Assert(lqph.size() == n_q_points, ExcInternalError());
+        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+          {
+            const double det_F_qp   = lqph[q_point]->get_det_F();
+            const double J_tilde_qp = lqph[q_point]->get_J_tilde();
+            const double the_error_qp_squared =
+              std::pow((det_F_qp - J_tilde_qp), 2);
+            const double JxW = fe_values_ref.JxW(q_point);
+            dil_L2_error += the_error_qp_squared * JxW;
+          }
+      }
+    return std::make_pair(std::sqrt(dil_L2_error),
+                          compute_vol_current() / vol_reference);
+  }
+  template <int dim>
+  void
+  Solid<dim>::get_error_residual(Errors &error_residual)
+  {
+    BlockVector<double> error_res(dofs_per_block);
+    for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
+      if (!constraints.is_constrained(i))
+        error_res(i) = system_rhs(i);
+    error_residual.norm = error_res.l2_norm();
+    error_residual.u    = error_res.block(u_dof).l2_norm();
+    error_residual.p    = error_res.block(p_dof).l2_norm();
+    error_residual.J    = error_res.block(J_dof).l2_norm();
+  }
+  template <int dim>
+  void
+  Solid<dim>::get_error_update(const BlockVector<double> &newton_update,
+                               Errors &                   error_update)
+  {
+    BlockVector<double> error_ud(dofs_per_block);
+    for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
+      if (!constraints.is_constrained(i))
+        error_ud(i) = newton_update(i);
+    error_update.norm = error_ud.l2_norm();
+    error_update.u    = error_ud.block(u_dof).l2_norm();
+    error_update.p    = error_ud.block(p_dof).l2_norm();
+    error_update.J    = error_ud.block(J_dof).l2_norm();
+  }
+  template <int dim>
+  BlockVector<double>
+  Solid<dim>::get_total_solution(
+    const BlockVector<double> &solution_delta) const
+  {
+    BlockVector<double> solution_total(solution_n);
+    solution_total += solution_delta;
+    return solution_total;
+  }
+  template <int dim>
+  void
+  Solid<dim>::assemble_system_tangent()
+  {
+    timer.enter_subsection("Assemble tangent matrix");
+    std::cout << " ASM_K " << std::flush;
+    tangent_matrix = 0.0;
+    const UpdateFlags uf_cell(update_values | update_gradients |
+                              update_JxW_values);
+    PerTaskData_K     per_task_data(dofs_per_cell);
+    ScratchData_K     scratch_data(fe, qf_cell, uf_cell);
+    WorkStream::run(dof_handler_ref.begin_active(),
+                    dof_handler_ref.end(),
+                    std::bind(&Solid<dim>::assemble_system_tangent_one_cell,
+                              this,
+                              std::placeholders::_1,
+                              std::placeholders::_2,
+                              std::placeholders::_3),
+                    std::bind(&Solid<dim>::copy_local_to_global_K,
+                              this,
+                              std::placeholders::_1),
+                    scratch_data,
+                    per_task_data);
+    timer.leave_subsection();
+  }
+  template <int dim>
+  void
+  Solid<dim>::copy_local_to_global_K(const PerTaskData_K &data)
+  {
+    for (unsigned int i = 0; i < dofs_per_cell; ++i)
+      for (unsigned int j = 0; j < dofs_per_cell; ++j)
+        tangent_matrix.add(data.local_dof_indices[i],
+                           data.local_dof_indices[j],
+                           data.cell_matrix(i, j));
+  }
+  template <int dim>
+  void
+  Solid<dim>::assemble_system_tangent_one_cell(
+    const typename DoFHandler<dim>::active_cell_iterator &cell,
+    ScratchData_K &                                       scratch,
+    PerTaskData_K &                                       data) const
+  {
+    data.reset();
+    scratch.reset();
+    scratch.fe_values_ref.reinit(cell);
+    cell->get_dof_indices(data.local_dof_indices);
+    const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
+      quadrature_point_history.get_data(cell);
+    Assert(lqph.size() == n_q_points, ExcInternalError());
+    for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+      {
+        const Tensor<2, dim> F_inv = lqph[q_point]->get_F_inv();
+        for (unsigned int k = 0; k < dofs_per_cell; ++k)
+          {
+            const unsigned int k_group = fe.system_to_base_index(k).first.first;
+            if (k_group == u_dof)
+              {
+                scratch.grad_Nx[q_point][k] =
+                  scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv;
+                scratch.symm_grad_Nx[q_point][k] =
+                  symmetrize(scratch.grad_Nx[q_point][k]);
+              }
+            else if (k_group == p_dof)
+              scratch.Nx[q_point][k] =
+                scratch.fe_values_ref[p_fe].value(k, q_point);
+            else if (k_group == J_dof)
+              scratch.Nx[q_point][k] =
+                scratch.fe_values_ref[J_fe].value(k, q_point);
+            else
+              Assert(k_group <= J_dof, ExcInternalError());
+          }
+      }
+    for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+      {
+        const Tensor<2, dim>          tau = lqph[q_point]->get_tau();
+        const SymmetricTensor<4, dim> Jc  = lqph[q_point]->get_Jc();
+        const double d2Psi_vol_dJ2        = lqph[q_point]->get_d2Psi_vol_dJ2();
+        const double det_F                = lqph[q_point]->get_det_F();
+        const std::vector<double> &                 N = scratch.Nx[q_point];
+        const std::vector<SymmetricTensor<2, dim>> &symm_grad_Nx =
+          scratch.symm_grad_Nx[q_point];
+        const std::vector<Tensor<2, dim>> &grad_Nx = scratch.grad_Nx[q_point];
+        const double JxW = scratch.fe_values_ref.JxW(q_point);
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          {
+            const unsigned int component_i =
+              fe.system_to_component_index(i).first;
+            const unsigned int i_group = fe.system_to_base_index(i).first.first;
+            for (unsigned int j = 0; j <= i; ++j)
+              {
+                const unsigned int component_j =
+                  fe.system_to_component_index(j).first;
+                const unsigned int j_group =
+                  fe.system_to_base_index(j).first.first;
+                if ((i_group == j_group) && (i_group == u_dof))
+                  {
+                    data.cell_matrix(i, j) += symm_grad_Nx[i] *
+                                              Jc // The material contribution:
+                                              * symm_grad_Nx[j] * JxW;
+                    if (component_i ==
+                        component_j) // geometrical stress contribution
+                      data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau *
+                                                grad_Nx[j][component_j] * JxW;
+                  }
+                else if ((i_group == p_dof) && (j_group == u_dof))
+                  {
+                    data.cell_matrix(i, j) +=
+                      N[i] * det_F *
+                      (symm_grad_Nx[j] * StandardTensors<dim>::I) * JxW;
+                  }
+                else if ((i_group == J_dof) && (j_group == p_dof))
+                  data.cell_matrix(i, j) -= N[i] * N[j] * JxW;
+                else if ((i_group == j_group) && (i_group == J_dof))
+                  data.cell_matrix(i, j) += N[i] * d2Psi_vol_dJ2 * N[j] * JxW;
+                else
+                  Assert((i_group <= J_dof) && (j_group <= J_dof),
+                         ExcInternalError());
+              }
+          }
+      }
+    for (unsigned int i = 0; i < dofs_per_cell; ++i)
+      for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
+        data.cell_matrix(i, j) = data.cell_matrix(j, i);
+  }
+  template <int dim>
+  void
+  Solid<dim>::assemble_system_rhs()
+  {
+    timer.enter_subsection("Assemble system right-hand side");
+    std::cout << " ASM_R " << std::flush;
+    system_rhs = 0.0;
+    const UpdateFlags uf_cell(update_values | update_gradients |
+                              update_JxW_values);
+    const UpdateFlags uf_face(update_values | update_normal_vectors |
+                              update_JxW_values);
+    PerTaskData_RHS   per_task_data(dofs_per_cell);
+    ScratchData_RHS   scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face);
+    WorkStream::run(dof_handler_ref.begin_active(),
+                    dof_handler_ref.end(),
+                    std::bind(&Solid<dim>::assemble_system_rhs_one_cell,
+                              this,
+                              std::placeholders::_1,
+                              std::placeholders::_2,
+                              std::placeholders::_3),
+                    std::bind(&Solid<dim>::copy_local_to_global_rhs,
+                              this,
+                              std::placeholders::_1),
+                    scratch_data,
+                    per_task_data);
+    timer.leave_subsection();
+  }
+  template <int dim>
+  void
+  Solid<dim>::copy_local_to_global_rhs(const PerTaskData_RHS &data)
+  {
+    for (unsigned int i = 0; i < dofs_per_cell; ++i)
+      system_rhs(data.local_dof_indices[i]) += data.cell_rhs(i);
+  }
+  template <int dim>
+  void
+  Solid<dim>::assemble_system_rhs_one_cell(
+    const typename DoFHandler<dim>::active_cell_iterator &cell,
+    ScratchData_RHS &                                     scratch,
+    PerTaskData_RHS &                                     data) const
+  {
+    data.reset();
+    scratch.reset();
+    scratch.fe_values_ref.reinit(cell);
+    cell->get_dof_indices(data.local_dof_indices);
+    const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
+      quadrature_point_history.get_data(cell);
+    Assert(lqph.size() == n_q_points, ExcInternalError());
+    for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+      {
+        const Tensor<2, dim> F_inv = lqph[q_point]->get_F_inv();
+        for (unsigned int k = 0; k < dofs_per_cell; ++k)
+          {
+            const unsigned int k_group = fe.system_to_base_index(k).first.first;
+            if (k_group == u_dof)
+              scratch.symm_grad_Nx[q_point][k] = symmetrize(
+                scratch.fe_values_ref[u_fe].gradient(k, q_point) * F_inv);
+            else if (k_group == p_dof)
+              scratch.Nx[q_point][k] =
+                scratch.fe_values_ref[p_fe].value(k, q_point);
+            else if (k_group == J_dof)
+              scratch.Nx[q_point][k] =
+                scratch.fe_values_ref[J_fe].value(k, q_point);
+            else
+              Assert(k_group <= J_dof, ExcInternalError());
+          }
+      }
+    for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+      {
+        const SymmetricTensor<2, dim> tau     = lqph[q_point]->get_tau();
+        const double                  det_F   = lqph[q_point]->get_det_F();
+        const double                  J_tilde = lqph[q_point]->get_J_tilde();
+        const double                  p_tilde = lqph[q_point]->get_p_tilde();
+        const double dPsi_vol_dJ = lqph[q_point]->get_dPsi_vol_dJ();
+        const std::vector<double> &                 N = scratch.Nx[q_point];
+        const std::vector<SymmetricTensor<2, dim>> &symm_grad_Nx =
+          scratch.symm_grad_Nx[q_point];
+        const double JxW = scratch.fe_values_ref.JxW(q_point);
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          {
+            const unsigned int i_group = fe.system_to_base_index(i).first.first;
+            if (i_group == u_dof)
+              data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
+            else if (i_group == p_dof)
+              data.cell_rhs(i) -= N[i] * (det_F - J_tilde) * JxW;
+            else if (i_group == J_dof)
+              data.cell_rhs(i) -= N[i] * (dPsi_vol_dJ - p_tilde) * JxW;
+            else
+              Assert(i_group <= J_dof, ExcInternalError());
+          }
+      }
+    for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
+         ++face)
+      if (cell->face(face)->at_boundary() == true &&
+          cell->face(face)->boundary_id() == 6)
+        {
+          scratch.fe_face_values_ref.reinit(cell, face);
+          for (unsigned int f_q_point = 0; f_q_point < n_q_points_f;
+               ++f_q_point)
+            {
+              const Tensor<1, dim> &N =
+                scratch.fe_face_values_ref.normal_vector(f_q_point);
+              static const double p0 =
+                -4.0 / (parameters.scale * parameters.scale);
+              const double         time_ramp = (time.current() / time.end());
+              const double         pressure  = p0 * parameters.p_p0 * time_ramp;
+              const Tensor<1, dim> traction  = pressure * N;
+              for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                {
+                  const unsigned int i_group =
+                    fe.system_to_base_index(i).first.first;
+                  if (i_group == u_dof)
+                    {
+                      const unsigned int component_i =
+                        fe.system_to_component_index(i).first;
+                      const double Ni =
+                        scratch.fe_face_values_ref.shape_value(i, f_q_point);
+                      const double JxW =
+                        scratch.fe_face_values_ref.JxW(f_q_point);
+                      data.cell_rhs(i) += (Ni * traction[component_i]) * JxW;
+                    }
+                }
+            }
+        }
+  }
+  template <int dim>
+  void
+  Solid<dim>::make_constraints(const int &it_nr)
+  {
+    std::cout << " CST " << std::flush;
+    if (it_nr > 1)
+      return;
+    constraints.clear();
+    const bool                       apply_dirichlet_bc = (it_nr == 0);
+    const FEValuesExtractors::Scalar x_displacement(0);
+    const FEValuesExtractors::Scalar y_displacement(1);
+    {
+      const int boundary_id = 0;
+      if (apply_dirichlet_bc == true)
+        VectorTools::interpolate_boundary_values(
+          dof_handler_ref,
+          boundary_id,
+          ZeroFunction<dim>(n_components),
+          constraints,
+          fe.component_mask(x_displacement));
+      else
+        VectorTools::interpolate_boundary_values(
+          dof_handler_ref,
+          boundary_id,
+          ZeroFunction<dim>(n_components),
+          constraints,
+          fe.component_mask(x_displacement));
+    }
+    {
+      const int boundary_id = 2;
+      if (apply_dirichlet_bc == true)
+        VectorTools::interpolate_boundary_values(
+          dof_handler_ref,
+          boundary_id,
+          ZeroFunction<dim>(n_components),
+          constraints,
+          fe.component_mask(y_displacement));
+      else
+        VectorTools::interpolate_boundary_values(
+          dof_handler_ref,
+          boundary_id,
+          ZeroFunction<dim>(n_components),
+          constraints,
+          fe.component_mask(y_displacement));
+    }
+    if (dim == 3)
+      {
+        const FEValuesExtractors::Scalar z_displacement(2);
+        {
+          const int boundary_id = 3;
+          if (apply_dirichlet_bc == true)
+            VectorTools::interpolate_boundary_values(
+              dof_handler_ref,
+              boundary_id,
+              ZeroFunction<dim>(n_components),
+              constraints,
+              (fe.component_mask(x_displacement) |
+               fe.component_mask(z_displacement)));
+          else
+            VectorTools::interpolate_boundary_values(
+              dof_handler_ref,
+              boundary_id,
+              ZeroFunction<dim>(n_components),
+              constraints,
+              (fe.component_mask(x_displacement) |
+               fe.component_mask(z_displacement)));
+        }
+        {
+          const int boundary_id = 4;
+          if (apply_dirichlet_bc == true)
+            VectorTools::interpolate_boundary_values(
+              dof_handler_ref,
+              boundary_id,
+              ZeroFunction<dim>(n_components),
+              constraints,
+              fe.component_mask(z_displacement));
+          else
+            VectorTools::interpolate_boundary_values(
+              dof_handler_ref,
+              boundary_id,
+              ZeroFunction<dim>(n_components),
+              constraints,
+              fe.component_mask(z_displacement));
+        }
+        {
+          const int boundary_id = 6;
+          if (apply_dirichlet_bc == true)
+            VectorTools::interpolate_boundary_values(
+              dof_handler_ref,
+              boundary_id,
+              ZeroFunction<dim>(n_components),
+              constraints,
+              (fe.component_mask(x_displacement) |
+               fe.component_mask(z_displacement)));
+          else
+            VectorTools::interpolate_boundary_values(
+              dof_handler_ref,
+              boundary_id,
+              ZeroFunction<dim>(n_components),
+              constraints,
+              (fe.component_mask(x_displacement) |
+               fe.component_mask(z_displacement)));
+        }
+      }
+    else
+      {
+        {
+          const int boundary_id = 3;
+          if (apply_dirichlet_bc == true)
+            VectorTools::interpolate_boundary_values(
+              dof_handler_ref,
+              boundary_id,
+              ZeroFunction<dim>(n_components),
+              constraints,
+              (fe.component_mask(x_displacement)));
+          else
+            VectorTools::interpolate_boundary_values(
+              dof_handler_ref,
+              boundary_id,
+              ZeroFunction<dim>(n_components),
+              constraints,
+              (fe.component_mask(x_displacement)));
+        }
+        {
+          const int boundary_id = 6;
+          if (apply_dirichlet_bc == true)
+            VectorTools::interpolate_boundary_values(
+              dof_handler_ref,
+              boundary_id,
+              ZeroFunction<dim>(n_components),
+              constraints,
+              (fe.component_mask(x_displacement)));
+          else
+            VectorTools::interpolate_boundary_values(
+              dof_handler_ref,
+              boundary_id,
+              ZeroFunction<dim>(n_components),
+              constraints,
+              (fe.component_mask(x_displacement)));
+        }
+      }
+    constraints.close();
+  }
+  template <int dim>
+  void
+  Solid<dim>::assemble_sc()
+  {
+    timer.enter_subsection("Perform static condensation");
+    std::cout << " ASM_SC " << std::flush;
+    PerTaskData_SC per_task_data(dofs_per_cell,
+                                 element_indices_u.size(),
+                                 element_indices_p.size(),
+                                 element_indices_J.size());
+    ScratchData_SC scratch_data;
+    WorkStream::run(dof_handler_ref.begin_active(),
+                    dof_handler_ref.end(),
+                    *this,
+                    &Solid::assemble_sc_one_cell,
+                    &Solid::copy_local_to_global_sc,
+                    scratch_data,
+                    per_task_data);
+    timer.leave_subsection();
+  }
+  template <int dim>
+  void
+  Solid<dim>::copy_local_to_global_sc(const PerTaskData_SC &data)
+  {
+    for (unsigned int i = 0; i < dofs_per_cell; ++i)
+      for (unsigned int j = 0; j < dofs_per_cell; ++j)
+        tangent_matrix.add(data.local_dof_indices[i],
+                           data.local_dof_indices[j],
+                           data.cell_matrix(i, j));
+  }
+  template <int dim>
+  void
+  Solid<dim>::assemble_sc_one_cell(
+    const typename DoFHandler<dim>::active_cell_iterator &cell,
+    ScratchData_SC &                                      scratch,
+    PerTaskData_SC &                                      data)
+  {
+    data.reset();
+    scratch.reset();
+    cell->get_dof_indices(data.local_dof_indices);
+    data.k_orig.extract_submatrix_from(tangent_matrix,
+                                       data.local_dof_indices,
+                                       data.local_dof_indices);
+    data.k_pu.extract_submatrix_from(data.k_orig,
+                                     element_indices_p,
+                                     element_indices_u);
+    data.k_pJ.extract_submatrix_from(data.k_orig,
+                                     element_indices_p,
+                                     element_indices_J);
+    data.k_JJ.extract_submatrix_from(data.k_orig,
+                                     element_indices_J,
+                                     element_indices_J);
+    data.k_pJ_inv.invert(data.k_pJ);
+    data.k_pJ_inv.mmult(data.A, data.k_pu);
+    data.k_JJ.mmult(data.B, data.A);
+    data.k_pJ_inv.Tmmult(data.C, data.B);
+    data.k_pu.Tmmult(data.k_bbar, data.C);
+    data.k_bbar.scatter_matrix_to(element_indices_u,
+                                  element_indices_u,
+                                  data.cell_matrix);
+    data.k_pJ_inv.add(-1.0, data.k_pJ);
+    data.k_pJ_inv.scatter_matrix_to(element_indices_p,
+                                    element_indices_J,
+                                    data.cell_matrix);
+  }
+  template <int dim>
+  std::pair<unsigned int, double>
+  Solid<dim>::solve_linear_system(BlockVector<double> &newton_update)
+  {
+    unsigned int lin_it  = 0;
+    double       lin_res = 0.0;
+    if (parameters.use_static_condensation == true)
+      {
+        BlockVector<double> A(dofs_per_block);
+        BlockVector<double> B(dofs_per_block);
+        {
+          assemble_sc();
+          tangent_matrix.block(p_dof, J_dof)
+            .vmult(A.block(J_dof), system_rhs.block(p_dof));
+          tangent_matrix.block(J_dof, J_dof)
+            .vmult(B.block(J_dof), A.block(J_dof));
+          A.block(J_dof) = system_rhs.block(J_dof);
+          A.block(J_dof) -= B.block(J_dof);
+          tangent_matrix.block(p_dof, J_dof)
+            .Tvmult(A.block(p_dof), A.block(J_dof));
+          tangent_matrix.block(u_dof, p_dof)
+            .vmult(A.block(u_dof), A.block(p_dof));
+          system_rhs.block(u_dof) -= A.block(u_dof);
+          timer.enter_subsection("Linear solver");
+          std::cout << " SLV " << std::flush;
+          if (parameters.type_lin == "CG")
+            {
+              const int solver_its = tangent_matrix.block(u_dof, u_dof).m() *
+                                     parameters.max_iterations_lin;
+              const double tol_sol =
+                parameters.tol_lin * system_rhs.block(u_dof).l2_norm();
+              SolverControl solver_control(solver_its, tol_sol, false, false);
+              GrowingVectorMemory<Vector<double>> GVM;
+              SolverCG<Vector<double>> solver_CG(solver_control, GVM);
+              PreconditionSelector<SparseMatrix<double>, Vector<double>>
+                preconditioner(parameters.preconditioner_type,
+                               parameters.preconditioner_relaxation);
+              preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof));
+              solver_CG.solve(tangent_matrix.block(u_dof, u_dof),
+                              newton_update.block(u_dof),
+                              system_rhs.block(u_dof),
+                              preconditioner);
+              lin_it  = solver_control.last_step();
+              lin_res = solver_control.last_value();
+            }
+          else if (parameters.type_lin == "Direct")
+            {
+              SparseDirectUMFPACK A_direct;
+              A_direct.initialize(tangent_matrix.block(u_dof, u_dof));
+              A_direct.vmult(newton_update.block(u_dof),
+                             system_rhs.block(u_dof));
+              lin_it  = 1;
+              lin_res = 0.0;
+            }
+          else
+            Assert(false, ExcMessage("Linear solver type not implemented"));
+          timer.leave_subsection();
+        }
+        constraints.distribute(newton_update);
+        timer.enter_subsection("Linear solver postprocessing");
+        std::cout << " PP " << std::flush;
+        {
+          tangent_matrix.block(p_dof, u_dof)
+            .vmult(A.block(p_dof), newton_update.block(u_dof));
+          A.block(p_dof) *= -1.0;
+          A.block(p_dof) += system_rhs.block(p_dof);
+          tangent_matrix.block(p_dof, J_dof)
+            .vmult(newton_update.block(J_dof), A.block(p_dof));
+        }
+        constraints.distribute(newton_update);
+        {
+          tangent_matrix.block(J_dof, J_dof)
+            .vmult(A.block(J_dof), newton_update.block(J_dof));
+          A.block(J_dof) *= -1.0;
+          A.block(J_dof) += system_rhs.block(J_dof);
+          tangent_matrix.block(p_dof, J_dof)
+            .Tvmult(newton_update.block(p_dof), A.block(J_dof));
+        }
+        constraints.distribute(newton_update);
+        timer.leave_subsection();
+      }
+    else
+      {
+        std::cout << " ------ " << std::flush;
+        timer.enter_subsection("Linear solver");
+        std::cout << " SLV " << std::flush;
+        if (parameters.type_lin == "CG")
+          {
+            const Vector<double> &f_u = system_rhs.block(u_dof);
+            const Vector<double> &f_p = system_rhs.block(p_dof);
+            const Vector<double> &f_J = system_rhs.block(J_dof);
+            Vector<double> &      d_u = newton_update.block(u_dof);
+            Vector<double> &      d_p = newton_update.block(p_dof);
+            Vector<double> &      d_J = newton_update.block(J_dof);
+            const auto            K_uu =
+              linear_operator(tangent_matrix.block(u_dof, u_dof));
+            const auto K_up =
+              linear_operator(tangent_matrix.block(u_dof, p_dof));
+            const auto K_pu =
+              linear_operator(tangent_matrix.block(p_dof, u_dof));
+            const auto K_Jp =
+              linear_operator(tangent_matrix.block(J_dof, p_dof));
+            const auto K_JJ =
+              linear_operator(tangent_matrix.block(J_dof, J_dof));
+            PreconditionSelector<SparseMatrix<double>, Vector<double>>
+              preconditioner_K_Jp_inv("jacobi");
+            preconditioner_K_Jp_inv.use_matrix(
+              tangent_matrix.block(J_dof, p_dof));
+            ReductionControl solver_control_K_Jp_inv(
+              tangent_matrix.block(J_dof, p_dof).m() *
+                parameters.max_iterations_lin,
+              1.0e-30,
+              parameters.tol_lin);
+            SolverSelector<Vector<double>> solver_K_Jp_inv;
+            solver_K_Jp_inv.select("cg");
+            solver_K_Jp_inv.set_control(solver_control_K_Jp_inv);
+            const auto K_Jp_inv =
+              inverse_operator(K_Jp, solver_K_Jp_inv, preconditioner_K_Jp_inv);
+            const auto K_pJ_inv     = transpose_operator(K_Jp_inv);
+            const auto K_pp_bar     = K_Jp_inv * K_JJ * K_pJ_inv;
+            const auto K_uu_bar_bar = K_up * K_pp_bar * K_pu;
+            const auto K_uu_con     = K_uu + K_uu_bar_bar;
+            PreconditionSelector<SparseMatrix<double>, Vector<double>>
+              preconditioner_K_con_inv(parameters.preconditioner_type,
+                                       parameters.preconditioner_relaxation);
+            preconditioner_K_con_inv.use_matrix(
+              tangent_matrix.block(u_dof, u_dof));
+            ReductionControl solver_control_K_con_inv(
+              tangent_matrix.block(u_dof, u_dof).m() *
+                parameters.max_iterations_lin,
+              1.0e-30,
+              parameters.tol_lin);
+            SolverSelector<Vector<double>> solver_K_con_inv;
+            solver_K_con_inv.select("cg");
+            solver_K_con_inv.set_control(solver_control_K_con_inv);
+            const auto K_uu_con_inv =
+              inverse_operator(K_uu_con,
+                               solver_K_con_inv,
+                               preconditioner_K_con_inv);
+            d_u =
+              K_uu_con_inv * (f_u - K_up * (K_Jp_inv * f_J - K_pp_bar * f_p));
+            timer.leave_subsection();
+            timer.enter_subsection("Linear solver postprocessing");
+            std::cout << " PP " << std::flush;
+            d_J     = K_pJ_inv * (f_p - K_pu * d_u);
+            d_p     = K_Jp_inv * (f_J - K_JJ * d_J);
+            lin_it  = solver_control_K_con_inv.last_step();
+            lin_res = solver_control_K_con_inv.last_value();
+          }
+        else if (parameters.type_lin == "Direct")
+          {
+            SparseDirectUMFPACK A_direct;
+            A_direct.initialize(tangent_matrix);
+            A_direct.vmult(newton_update, system_rhs);
+            lin_it  = 1;
+            lin_res = 0.0;
+            std::cout << " -- " << std::flush;
+          }
+        else
+          Assert(false, ExcMessage("Linear solver type not implemented"));
+        timer.leave_subsection();
+        constraints.distribute(newton_update);
+      }
+    return std::make_pair(lin_it, lin_res);
+  }
+  template <int dim>
+  void
+  Solid<dim>::output_results() const
+  {
+    DataOut<dim> data_out;
+    std::vector<DataComponentInterpretation::DataComponentInterpretation>
+      data_component_interpretation(
+        dim, DataComponentInterpretation::component_is_part_of_vector);
+    data_component_interpretation.push_back(
+      DataComponentInterpretation::component_is_scalar);
+    data_component_interpretation.push_back(
+      DataComponentInterpretation::component_is_scalar);
+    std::vector<std::string> solution_name(dim, "displacement");
+    solution_name.push_back("pressure");
+    solution_name.push_back("dilatation");
+    data_out.attach_dof_handler(dof_handler_ref);
+    data_out.add_data_vector(solution_n,
+                             solution_name,
+                             DataOut<dim>::type_dof_data,
+                             data_component_interpretation);
+    Vector<double> soln(solution_n.size());
+    for (unsigned int i = 0; i < soln.size(); ++i)
+      soln(i) = solution_n(i);
+    MappingQEulerian<dim> q_mapping(degree, dof_handler_ref, soln);
+    data_out.build_patches(q_mapping, degree);
+    std::ostringstream filename;
+    filename << "solution-" << dim << "d-" << time.get_timestep() << ".vtk";
+    std::ofstream output(filename.str().c_str());
+    data_out.write_vtk(output);
+  }
+} // namespace Step44
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+
+  using namespace dealii;
+  using namespace Step44;
+  try
+    {
+      const unsigned int dim = 3;
+      Solid<dim>         solid(SOURCE_DIR "/prm/parameters-step-44.prm");
+      solid.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  return 0;
+}
diff --git a/tests/adolc/step-44.output b/tests/adolc/step-44.output
new file mode 100644 (file)
index 0000000..9f2261a
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::Timestep 1: 0.00 -0.000129 0.00
+DEAL::Timestep 2: 0.00 -0.000256 0.00
+DEAL::Timestep 3: 0.00 -0.000376 0.00
+DEAL::Timestep 4: 0.00 -0.000485 0.00
+DEAL::Timestep 5: 0.00 -0.000580 0.00
+DEAL::Timestep 6: 0.00 -0.000660 0.00
+DEAL::Timestep 7: 0.00 -0.000727 0.00
+DEAL::Timestep 8: 0.00 -0.000783 0.00
+DEAL::Timestep 9: 0.00 -0.000828 0.00
+DEAL::Timestep 10: 0.00 -0.000866 0.00
diff --git a/tests/adolc/symmetric_tensor_functions_02_1.cc b/tests/adolc/symmetric_tensor_functions_02_1.cc
new file mode 100644 (file)
index 0000000..c3ba107
--- /dev/null
@@ -0,0 +1,40 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test to check that tensor functions both compile and produce the right
+// result when differentiated using the various auto-differentiable number
+// types: Tensor inverse
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/symmetric_tensor_functions_02.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_symmetric_tensor<2, double, AD::NumberTypes::adolc_taped>();
+    test_symmetric_tensor<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/symmetric_tensor_functions_02_1.output b/tests/adolc/symmetric_tensor_functions_02_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/symmetric_tensor_functions_02_2.cc b/tests/adolc/symmetric_tensor_functions_02_2.cc
new file mode 100644 (file)
index 0000000..259532b
--- /dev/null
@@ -0,0 +1,40 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test to check that tensor functions both compile and produce the right
+// result when differentiated using the various auto-differentiable number
+// types: Tensor inverse
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/symmetric_tensor_functions_02.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_symmetric_tensor<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_symmetric_tensor<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/symmetric_tensor_functions_02_2.output b/tests/adolc/symmetric_tensor_functions_02_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/symmetric_tensor_functions_03_1.cc b/tests/adolc/symmetric_tensor_functions_03_1.cc
new file mode 100644 (file)
index 0000000..03a52ad
--- /dev/null
@@ -0,0 +1,40 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test to check that tensor functions both compile and produce the right
+// result when differentiated using the various auto-differentiable number
+// types: Eigenvalues
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/symmetric_tensor_functions_03.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_symmetric_tensor<2, double, AD::NumberTypes::adolc_taped>();
+    test_symmetric_tensor<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/symmetric_tensor_functions_03_1.with_adolc_with_advanced_branching=false.output b/tests/adolc/symmetric_tensor_functions_03_1.with_adolc_with_advanced_branching=false.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/symmetric_tensor_functions_03_1.with_adolc_with_advanced_branching=true.expect=run.output b/tests/adolc/symmetric_tensor_functions_03_1.with_adolc_with_advanced_branching=true.expect=run.output
new file mode 100644 (file)
index 0000000..8b13789
--- /dev/null
@@ -0,0 +1 @@
+
diff --git a/tests/adolc/symmetric_tensor_functions_03_2.cc b/tests/adolc/symmetric_tensor_functions_03_2.cc
new file mode 100644 (file)
index 0000000..a0a74dd
--- /dev/null
@@ -0,0 +1,40 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test to check that tensor functions both compile and produce the right
+// result when differentiated using the various auto-differentiable number
+// types: Eigenvalues
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/symmetric_tensor_functions_03.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_symmetric_tensor<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_symmetric_tensor<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/symmetric_tensor_functions_03_2.output b/tests/adolc/symmetric_tensor_functions_03_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/symmetric_tensor_functions_04_1.cc b/tests/adolc/symmetric_tensor_functions_04_1.cc
new file mode 100644 (file)
index 0000000..735cf58
--- /dev/null
@@ -0,0 +1,40 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test to check that tensor functions both compile and produce the right
+// result when differentiated using the various auto-differentiable number
+// types: Eigenvalues and eignvectors
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/symmetric_tensor_functions_04.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_symmetric_tensor<2, double, AD::NumberTypes::adolc_taped>();
+    test_symmetric_tensor<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/symmetric_tensor_functions_04_1.with_adolc_with_advanced_branching=false.output b/tests/adolc/symmetric_tensor_functions_04_1.with_adolc_with_advanced_branching=false.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/symmetric_tensor_functions_04_1.with_adolc_with_advanced_branching=true.expect=run.output b/tests/adolc/symmetric_tensor_functions_04_1.with_adolc_with_advanced_branching=true.expect=run.output
new file mode 100644 (file)
index 0000000..8b13789
--- /dev/null
@@ -0,0 +1 @@
+
diff --git a/tests/adolc/symmetric_tensor_functions_04_2.cc b/tests/adolc/symmetric_tensor_functions_04_2.cc
new file mode 100644 (file)
index 0000000..1b03ca8
--- /dev/null
@@ -0,0 +1,40 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test to check that tensor functions both compile and produce the right
+// result when differentiated using the various auto-differentiable number
+// types: Eigenvalues and eignvectors
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/symmetric_tensor_functions_04.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_symmetric_tensor<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_symmetric_tensor<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/symmetric_tensor_functions_04_2.output b/tests/adolc/symmetric_tensor_functions_04_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/tensor_functions_02_1.cc b/tests/adolc/tensor_functions_02_1.cc
new file mode 100644 (file)
index 0000000..019aca9
--- /dev/null
@@ -0,0 +1,40 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test to check that tensor functions both compile and produce the right
+// result when differentiated using the various auto-differentiable number
+// types: Tensor inverse
+//
+// AD number type: ADOL-C taped
+
+#include "../ad_common_tests/tensor_functions_02.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_tensor<2, double, AD::NumberTypes::adolc_taped>();
+    test_tensor<3, double, AD::NumberTypes::adolc_taped>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/tensor_functions_02_1.output b/tests/adolc/tensor_functions_02_1.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/tensor_functions_02_2.cc b/tests/adolc/tensor_functions_02_2.cc
new file mode 100644 (file)
index 0000000..632453f
--- /dev/null
@@ -0,0 +1,40 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Test to check that tensor functions both compile and produce the right
+// result when differentiated using the various auto-differentiable number
+// types: Tensor inverse
+//
+// AD number type: ADOL-C tapeless
+
+#include "../ad_common_tests/tensor_functions_02.h"
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Double");
+  {
+    test_tensor<2, double, AD::NumberTypes::adolc_tapeless>();
+    test_tensor<3, double, AD::NumberTypes::adolc_tapeless>();
+    deallog << "OK" << std::endl;
+  }
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/tensor_functions_02_2.output b/tests/adolc/tensor_functions_02_2.output
new file mode 100644 (file)
index 0000000..3c6dd25
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Double:OK
+DEAL::OK
diff --git a/tests/adolc/tensor_ops_01.cc b/tests/adolc/tensor_ops_01.cc
new file mode 100644 (file)
index 0000000..11295a5
--- /dev/null
@@ -0,0 +1,236 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// compilation check constructor and operations acting on tensors of ADOL-C
+// number types
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/differentiation/ad.h>
+
+#include <fstream>
+#include <iomanip>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim,
+          typename number_t,
+          enum Differentiation::AD::NumberTypes ad_number_enum>
+void
+test_tensor()
+{
+  typedef typename Differentiation::AD::NumberTraits<number_t,
+                                                     ad_number_enum>::ad_type
+                                      ad_number_t;
+  typedef Tensor<2, dim, ad_number_t> AD_Tensor;
+  typedef Tensor<2, dim, number_t>    NonAD_Tensor;
+
+  // Constructors
+  NonAD_Tensor t1;
+  AD_Tensor    adt1, adt2;
+  AD_Tensor    adt3(adt1);
+
+  // Assignment
+  for (unsigned int i = 0; i < dim; ++i)
+    {
+      adt1[i][i] += 1.0;
+      for (unsigned int j = 0; j < dim; ++j)
+        {
+          t1[i][j] = dim * i + j;
+          adt1[i][j] += 2 * dim * i + j + 1;
+          adt2[i][j] += 3 * dim * i + j + 2;
+        }
+    }
+  AD_Tensor adt4 = adt1;
+
+  // Comparison
+  const bool res_1 = (adt1 == adt2);
+  const bool res_2 = (adt1 != adt2);
+
+  // Scalar operations
+  adt4 *= 2.0;
+  adt4 /= 4.0;
+  adt4[0][0] += 1.0;
+  adt4[0][0] *= 3.0;
+
+  // Vector space operations
+  adt4 = 2.0 * adt1;
+  adt4 += adt2;
+  adt4 -= adt2;
+  adt4 += t1;
+  adt4 -= t1;
+  const AD_Tensor adt5 = adt1 + adt2;
+  const AD_Tensor adt6 = adt1 - adt2;
+
+  // Contractions
+  const ad_number_t ad_res1 = scalar_product(adt1, adt2);
+  const ad_number_t ad_res2 = double_contract(adt1, adt2);
+  //  ad_res1 = double_contract(adt1,t1); // TODO: Not defined. Conflicting
+  //  number types ad_res1 = double_contract(t1,adt2); // TODO: Not defined.
+  //  Conflicting number types
+  AD_Tensor adt7 = adt1 * adt2;
+  adt7           = t1 * adt2;
+  adt7           = adt1 * t1;
+
+  // Outer product
+  const Tensor<1, dim, number_t>    v1;
+  const Tensor<1, dim, ad_number_t> av1;
+  const AD_Tensor                   adt8  = outer_product(v1, av1);
+  const AD_Tensor                   adt9  = outer_product(av1, v1);
+  const AD_Tensor                   adt10 = outer_product(av1, av1);
+
+  // Special operations
+  const ad_number_t det = determinant(adt1);
+  const ad_number_t tr  = trace(adt1);
+  //  const ad_number_t l1 = l1_norm(adt1);       // TODO: Defined with std::
+  //  math operator const ad_number_t linf = linfty_norm(adt1); // TODO: Defined
+  //  with std:: math operator
+  const AD_Tensor inv   = invert(adt1);
+  const AD_Tensor trans = transpose(adt1);
+}
+
+template <int dim,
+          typename number_t,
+          enum Differentiation::AD::NumberTypes ad_number_enum>
+void
+test_symmetric_tensor()
+{
+  typedef typename Differentiation::AD::NumberTraits<number_t,
+                                                     ad_number_enum>::ad_type
+                                               ad_number_t;
+  typedef SymmetricTensor<2, dim, ad_number_t> AD_STensor;
+  typedef SymmetricTensor<4, dim, ad_number_t> AD_STensor4;
+  typedef SymmetricTensor<2, dim, number_t>    NonAD_STensor;
+  typedef SymmetricTensor<4, dim, number_t>    NonAD_STensor4;
+
+  // Constructors
+  NonAD_STensor  t1;
+  NonAD_STensor4 t2;
+  AD_STensor     adt1, adt2;
+  AD_STensor     adt3(adt1);
+
+  // Assignment
+  for (unsigned int i = 0; i < dim; ++i)
+    {
+      for (unsigned int j = i; j < dim; ++j)
+        {
+          t1[i][j]   = dim * i + j;
+          adt1[i][j] = 2 * dim * i + j + 1;
+          adt2[i][j] += 3 * dim * i + j + 2;
+        }
+    }
+  AD_STensor adt4 = adt1;
+
+  // Comparison
+  bool res_1 = (adt1 == adt2);
+  bool res_2 = (adt1 != adt2);
+
+  // Scalar operations
+  adt4 *= 2.0;
+  adt4 /= 4.0;
+  adt4 = 2.0 * adt1;
+  adt4[0][0] += 1.0;
+  adt4[0][0] -= 0.5;
+  adt4[0][0] *= 3.0;
+  adt4[0][0] /= 1.5;
+
+  // Vector space operations
+  adt4 += adt2;
+  adt4 -= adt2;
+  adt4 += t1;
+  adt4 -= t1;
+  const AD_STensor adt5 = adt1 + adt2;
+  const AD_STensor adt6 = adt1 - adt2;
+
+  // Contractions
+  const ad_number_t ad_res1 = scalar_product(adt1, adt2);
+  const ad_number_t ad_res2 = adt1 * adt2;
+  //  const ad_number_t ad_res3 = scalar_product(t1,adt2);
+  //  const ad_number_t ad_res4 = scalar_product(adt1,t1);
+  // const ad_number_t ad_res5 = t1*adt2; // TODO: Overload issue
+  //  const ad_number_t ad_res6 = adt1*t1; // TODO: Overload issue
+
+  // Outer product
+  const Tensor<1, dim, number_t>    v1;
+  const Tensor<1, dim, ad_number_t> av1;
+  const AD_STensor                  adt8  = symmetrize(outer_product(v1, av1));
+  const AD_STensor                  adt9  = symmetrize(outer_product(av1, v1));
+  const AD_STensor                  adt10 = symmetrize(outer_product(av1, av1));
+
+  // Tensor outer product
+  const AD_STensor4 adt11 = outer_product(adt1, adt2);
+  // const AD_STensor4 adt12 = outer_product(t1,adt1); // TODO: Cannot deduce
+  // return type
+
+  // Special operations
+  const ad_number_t det  = determinant(adt1);
+  const ad_number_t tr   = trace(adt1);
+  const ad_number_t inv1 = first_invariant(adt1);
+  const ad_number_t inv2 = second_invariant(adt1);
+  const ad_number_t inv3 = third_invariant(adt1); // TODO: Return type incorrect
+  // const ad_number_t l1 = l1_norm(adt1); // TODO: Defined with std:: math
+  // operator const ad_number_t linf = linfty_norm(adt1); // TODO: Defined with
+  // std:: math operator
+  const AD_STensor inv   = invert(adt1);
+  const AD_STensor trans = transpose(adt1);
+  const AD_STensor dev   = deviator(adt1);
+
+  // Special tensors
+  const SymmetricTensor<2, dim, ad_number_t> ust =
+    unit_symmetric_tensor<dim, ad_number_t>();
+  const SymmetricTensor<4, dim, ad_number_t> dt =
+    deviator_tensor<dim, ad_number_t>();
+  //  SymmetricTensor<4,dim,ad_number_t> idt =
+  //  identity_tensor<4,dim,ad_number_t> (); // TODO: No number type
+  const SymmetricTensor<4, dim, ad_number_t> op = outer_product(adt1, adt2);
+  //  SymmetricTensor<4,dim,ad_number_t> inv_op = invert(op); // TODO: Not
+  //  implemented!
+}
+
+int
+main()
+{
+  initlog();
+
+  // It is assumed that whatever works for rank-2
+  // tensors will work for the other ranks
+
+  // --- Taped ---
+  // double type
+  test_tensor<3, double, Differentiation::AD::NumberTypes::adolc_taped>();
+  test_symmetric_tensor<3,
+                        double,
+                        Differentiation::AD::NumberTypes::adolc_taped>();
+  // complex double
+  //  test_tensor<3,std::complex<double>,Differentiation::AD::NumberTypes::adolc_taped>();
+  //  test_symmetric_tensor<3,std::complex<double>,Differentiation::AD::NumberTypes::adolc_taped>();
+
+  // --- Tapeless ---
+  // double type
+  test_tensor<3, double, Differentiation::AD::NumberTypes::adolc_tapeless>();
+  test_symmetric_tensor<3,
+                        double,
+                        Differentiation::AD::NumberTypes::adolc_tapeless>();
+  // complex double
+  //  test_tensor<3,std::complex<double>,Differentiation::AD::NumberTypes::adolc_tapeless>();
+  //  test_symmetric_tensor<3,std::complex<double>,Differentiation::AD::NumberTypes::adolc_tapeless>();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=false.output b/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=false.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.expect=run.debug.output b/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.expect=run.debug.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.release.output b/tests/adolc/tensor_ops_01.with_adolc_with_advanced_branching=true.release.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/adolc/tensor_ops_02.cc b/tests/adolc/tensor_ops_02.cc
new file mode 100644 (file)
index 0000000..020da8a
--- /dev/null
@@ -0,0 +1,105 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Check ability to cast from a tensor of ADOL-C numbers to floats
+// In particular, this checks that the specializations for
+// internal::NumberType<Number>::value() work as expected when
+// Number == double and the input to value() is an ADOL-C number.
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/differentiation/ad.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim,
+          typename number_t,
+          enum Differentiation::AD::NumberTypes ad_number_enum>
+void
+test_tensor()
+{
+  typedef typename Differentiation::AD::NumberTraits<number_t,
+                                                     ad_number_enum>::ad_type
+                                      ad_number_t;
+  typedef Tensor<2, dim, ad_number_t> AD_Tensor;
+  typedef Tensor<2, dim, number_t>    NonAD_Tensor;
+
+  AD_Tensor adt1;
+  for (unsigned int i = 0; i < NonAD_Tensor::n_independent_components; ++i)
+    adt1[NonAD_Tensor::unrolled_to_component_indices(i)] = i + 1;
+
+  const NonAD_Tensor t1 = NonAD_Tensor(adt1);
+
+  Assert(t1.norm() > 0.0, ExcMessage("Cast and copy unsuccessful"));
+}
+
+template <int dim,
+          typename number_t,
+          enum Differentiation::AD::NumberTypes ad_number_enum>
+void
+test_symmetric_tensor()
+{
+  typedef typename Differentiation::AD::NumberTraits<number_t,
+                                                     ad_number_enum>::ad_type
+                                               ad_number_t;
+  typedef SymmetricTensor<2, dim, ad_number_t> AD_STensor2;
+  typedef SymmetricTensor<4, dim, ad_number_t> AD_STensor4;
+  typedef SymmetricTensor<2, dim, number_t>    NonAD_STensor2;
+  typedef SymmetricTensor<4, dim, number_t>    NonAD_STensor4;
+
+  // Constructors
+  AD_STensor2 adt1;
+  AD_STensor4 adt2;
+
+  for (unsigned int i = 0; i < NonAD_STensor2::n_independent_components; ++i)
+    adt1[NonAD_STensor2::unrolled_to_component_indices(i)] = i + 1;
+
+  for (unsigned int i = 0; i < NonAD_STensor4::n_independent_components; ++i)
+    adt2.access_raw_entry(i) = i + 1;
+
+  const NonAD_STensor2 t1 = NonAD_STensor2(adt1);
+  const NonAD_STensor4 t2 = NonAD_STensor4(adt2);
+
+  Assert(t1.norm() > 0.0, ExcMessage("Cast and copy unsuccessful"));
+  Assert(t2.norm() > 0.0, ExcMessage("Cast and copy unsuccessful"));
+}
+
+int
+main()
+{
+  initlog();
+
+  // It is assumed that whatever works for rank-2
+  // tensors will work for the other ranks
+
+  // --- Taped ---
+  test_tensor<3, double, Differentiation::AD::NumberTypes::adolc_taped>();
+  test_symmetric_tensor<3,
+                        double,
+                        Differentiation::AD::NumberTypes::adolc_taped>();
+
+  // --- Tapeless ---
+  test_tensor<3, double, Differentiation::AD::NumberTypes::adolc_tapeless>();
+  test_symmetric_tensor<3,
+                        double,
+                        Differentiation::AD::NumberTypes::adolc_tapeless>();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/adolc/tensor_ops_02.output b/tests/adolc/tensor_ops_02.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK

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.