diff --git a/src/common/common.hh b/src/common/common.hh
index 0e970bd..ddac2f0 100644
--- a/src/common/common.hh
+++ b/src/common/common.hh
@@ -1,311 +1,319 @@
 /**
  * @file   common.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   01 May 2017
  *
  * @brief  Small definitions of commonly used types throughout µSpectre
  *
  * @section  LICENSE
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can 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 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with µSpectre; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * * Boston, MA 02111-1307, USA.
  *
  * Additional permission under GNU GPL version 3 section 7
  *
  * If you modify this Program, or any covered work, by linking or combining it
  * with proprietary FFT implementations or numerical libraries, containing parts
  * covered by the terms of those libraries' licenses, the licensors of this
  * Program grant you additional permission to convey the resulting work.
  */
 
 #include <array>
 #include <cmath>
 #include <complex>
 #include <iostream>
 #include <string>
 #include <type_traits>
 
 #ifndef SRC_COMMON_COMMON_HH_
 #define SRC_COMMON_COMMON_HH_
 
 namespace muSpectre {
 
   /**
    * Eigen uses signed integers for dimensions. For consistency,
    µSpectre uses them througout the code. needs to represent -1 for
    eigen
    */
   using Dim_t = int;
 
   constexpr Dim_t oneD{1};         //!< constant for a one-dimensional problem
   constexpr Dim_t twoD{2};         //!< constant for a two-dimensional problem
   constexpr Dim_t threeD{3};       //!< constant for a three-dimensional problem
   constexpr Dim_t firstOrder{1};   //!< constant for vectors
   constexpr Dim_t secondOrder{2};  //!< constant second-order tensors
   constexpr Dim_t fourthOrder{4};  //!< constant fourth-order tensors
 
   //@{
   //! @anchor scalars
   //! Scalar types used for mathematical calculations
   using Uint = unsigned int;
   using Int = int;
   using Real = double;
   using Complex = std::complex<Real>;
   //@}
 
   //! Ccoord_t are cell coordinates, i.e. integer coordinates
   template <Dim_t dim> using Ccoord_t = std::array<Dim_t, dim>;
   //! Real space coordinates
   template <Dim_t dim> using Rcoord_t = std::array<Real, dim>;
 
   /**
    * Allows inserting `muSpectre::Ccoord_t` and `muSpectre::Rcoord_t`
    * into `std::ostream`s
    */
   template <typename T, size_t dim>
   std::ostream & operator<<(std::ostream & os,
                             const std::array<T, dim> & index) {
     os << "(";
     for (size_t i = 0; i < dim - 1; ++i) {
       os << index[i] << ", ";
     }
     os << index.back() << ")";
     return os;
   }
 
   //! element-wise division
   template <size_t dim>
   Rcoord_t<dim> operator/(const Rcoord_t<dim> & a, const Rcoord_t<dim> & b) {
     Rcoord_t<dim> retval{a};
     for (size_t i = 0; i < dim; ++i) {
       retval[i] /= b[i];
     }
     return retval;
   }
 
   //! element-wise division
   template <size_t dim>
   Rcoord_t<dim> operator/(const Rcoord_t<dim> & a, const Ccoord_t<dim> & b) {
     Rcoord_t<dim> retval{a};
     for (size_t i = 0; i < dim; ++i) {
       retval[i] /= b[i];
     }
     return retval;
   }
 
   //! convenience definitions
   constexpr Real pi{3.1415926535897932384626433};
 
   //! compile-time potentiation required for field-size computations
   template <typename R, typename I> constexpr R ipow(R base, I exponent) {
     static_assert(std::is_integral<I>::value, "Type must be integer");
     R retval{1};
     for (I i = 0; i < exponent; ++i) {
       retval *= base;
     }
     return retval;
   }
 
   /**
    * Copyright banner to be printed to the terminal by executables
    * Arguments are the executable's name, year of writing and the name
    * + address of the copyright holder
    */
   void banner(std::string name, Uint year, std::string cpy_holder);
 
   /**
    * Planner flags for FFT (follows FFTW, hopefully this choice will
    * be compatible with alternative FFT implementations)
    * @enum muSpectre::FFT_PlanFlags
    */
   enum class FFT_PlanFlags {
     estimate,  //!< cheapest plan for slowest execution
     measure,   //!< more expensive plan for fast execution
     patient    //!< very expensive plan for fastest execution
   };
 
   //! continuum mechanics flags
   enum class Formulation {
     finite_strain,    //!< causes evaluation in PK1(F)
     small_strain,     //!< causes evaluation in   σ(ε)
     small_strain_sym  //!< symmetric storage as vector ε
   };
 
+  //! finite differences flags
+  enum class FiniteDiff {
+    forward,   //!< ∂f/∂x ≈ (f(x+Δx) - f(x))/Δx
+    backward,  //!< ∂f/∂x ≈ (f(x) - f(x-Δx))/Δx
+    centred    //!< ∂f/∂x ≈ (f(x+Δx) - f(x-Δx))/2Δx
+  };
+
+
   /**
    * compile time computation of voigt vector
    */
   template <bool sym = true> constexpr Dim_t vsize(Dim_t dim) {
     if (sym) {
       return (dim * (dim - 1) / 2 + dim);
     } else {
       return dim * dim;
     }
   }
 
   //! compute the number of degrees of freedom to store for the strain
   //! tenor given dimension dim
   constexpr Dim_t dof_for_formulation(const Formulation form, const Dim_t dim) {
     switch (form) {
     case Formulation::small_strain_sym: {
       return vsize(dim);
       break;
     }
     default:
       return ipow(dim, 2);
       break;
     }
   }
 
   //! inserts `muSpectre::Formulation`s into `std::ostream`s
   std::ostream & operator<<(std::ostream & os, Formulation f);
 
   /* ---------------------------------------------------------------------- */
   //! Material laws can declare which type of stress measure they provide,
   //! and µSpectre will handle conversions
   enum class StressMeasure {
     Cauchy,     //!< Cauchy stress σ
     PK1,        //!< First Piola-Kirchhoff stress
     PK2,        //!< Second Piola-Kirchhoff stress
     Kirchhoff,  //!< Kirchhoff stress τ
     Biot,       //!< Biot stress
     Mandel,     //!< Mandel stress
     no_stress_  //!< only for triggering static_asserts
   };
   //! inserts `muSpectre::StressMeasure`s into `std::ostream`s
   std::ostream & operator<<(std::ostream & os, StressMeasure s);
 
   /* ---------------------------------------------------------------------- */
   //! Material laws can declare which type of strain measure they require and
   //! µSpectre will provide it
   enum class StrainMeasure {
     Gradient,       //!< placement gradient (δy/δx)
     Infinitesimal,  //!< small strain tensor .5(∇u + ∇uᵀ)
     GreenLagrange,  //!< Green-Lagrange strain .5(Fᵀ·F - I)
     Biot,           //!< Biot strain
     Log,            //!< logarithmic strain
     Almansi,        //!< Almansi strain
     RCauchyGreen,   //!< Right Cauchy-Green tensor
     LCauchyGreen,   //!< Left Cauchy-Green tensor
     no_strain_      //!< only for triggering static_assert
   };
   //! inserts `muSpectre::StrainMeasure`s into `std::ostream`s
   std::ostream & operator<<(std::ostream & os, StrainMeasure s);
 
   /* ---------------------------------------------------------------------- */
   /**
    * all isotropic elastic moduli to identify conversions, such as E
    * = µ(3λ + 2µ)/(λ+µ). For the full description, see
    * https://en.wikipedia.org/wiki/Lam%C3%A9_parameters
    * Not all the conversions are implemented, so please add as needed
    */
   enum class ElasticModulus {
     Bulk,          //!< Bulk modulus K
     K = Bulk,      //!< alias for ``ElasticModulus::Bulk``
     Young,         //!< Young's modulus E
     E = Young,     //!< alias for ``ElasticModulus::Young``
     lambda,        //!< Lamé's first parameter λ
     Shear,         //!< Shear modulus G or µ
     G = Shear,     //!< alias for ``ElasticModulus::Shear``
     mu = Shear,    //!< alias for ``ElasticModulus::Shear``
     Poisson,       //!< Poisson's ratio ν
     nu = Poisson,  //!< alias for ``ElasticModulus::Poisson``
     Pwave,         //!< P-wave modulus M
     M = Pwave,     //!< alias for ``ElasticModulus::Pwave``
     no_modulus_
   };  //!< only for triggering static_asserts
 
   /**
    * define comparison in order to exploit that moduli can be
    * expressed in terms of any two other moduli in any order (e.g. K
    * = K(E, ν) = K(ν, E)
    */
   constexpr inline bool operator<(ElasticModulus A, ElasticModulus B) {
     return static_cast<int>(A) < static_cast<int>(B);
   }
   /* ---------------------------------------------------------------------- */
   /** Compile-time function to g strain measure stored by muSpectre
       depending on the formulation
    **/
   constexpr StrainMeasure get_stored_strain_type(Formulation form) {
     switch (form) {
     case Formulation::finite_strain: {
       return StrainMeasure::Gradient;
       break;
     }
     case Formulation::small_strain: {
       return StrainMeasure::Infinitesimal;
       break;
     }
     default:
       return StrainMeasure::no_strain_;
       break;
     }
   }
 
   /** Compile-time function to g stress measure stored by muSpectre
       depending on the formulation
    **/
   constexpr StressMeasure get_stored_stress_type(Formulation form) {
     switch (form) {
     case Formulation::finite_strain: {
       return StressMeasure::PK1;
       break;
     }
     case Formulation::small_strain: {
       return StressMeasure::Cauchy;
       break;
     }
     default:
       return StressMeasure::no_stress_;
       break;
     }
   }
 
   /* ---------------------------------------------------------------------- */
   /** Compile-time functions to get the stress and strain measures
       after they may have been modified by choosing a formulation.
 
       For instance, a law that expecs a Green-Lagrange strain as input
       will get the infinitesimal strain tensor instead in a small
       strain computation
    **/
   constexpr StrainMeasure get_formulation_strain_type(Formulation form,
                                                       StrainMeasure expected) {
     switch (form) {
     case Formulation::finite_strain: {
       return expected;
       break;
     }
     case Formulation::small_strain: {
       return get_stored_strain_type(form);
       break;
     }
     default:
       return StrainMeasure::no_strain_;
       break;
     }
   }
 
 }  // namespace muSpectre
 
 #ifndef EXPLICITLY_TURNED_ON_CXX17
 #include "common/utilities.hh"
 #endif
 
 #endif  // SRC_COMMON_COMMON_HH_
diff --git a/src/materials/material_evaluator.hh b/src/materials/material_evaluator.hh
index beec9d5..615253a 100644
--- a/src/materials/material_evaluator.hh
+++ b/src/materials/material_evaluator.hh
@@ -1,184 +1,267 @@
 /**
  * @file   material_evaluator.hh
  *
  * @author Till Junge <till.junge@altermail.ch>
  *
  * @date   12 Dec 2018
  *
  * @brief  Helper to evaluate material laws
  *
  * Copyright © 2018 Till Junge
  *
  * µSpectre is free software; you can 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 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with µSpectre; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  *
  * Additional permission under GNU GPL version 3 section 7
  *
  * If you modify this Program, or any covered work, by linking or combining it
  * with proprietary FFT implementations or numerical libraries, containing parts
  * covered by the terms of those libraries' licenses, the licensors of this
  * Program grant you additional permission to convey the resulting work.
  */
 
 #ifndef SRC_MATERIALS_MATERIAL_EVALUATOR_HH_
 #define SRC_MATERIALS_MATERIAL_EVALUATOR_HH_
 
 #include "common/common.hh"
 #include "common/T4_map_proxy.hh"
 #include "common/ccoord_operations.hh"
 #include "materials/materials_toolbox.hh"
 #include "common/field.hh"
 
 namespace muSpectre {
 
   template <class Material>
   class MaterialEvaluator {
    public:
     static constexpr Dim_t DimM{Material::mdim()};
 
     using traits = typename Material::traits;
     using T2_t = Eigen::Matrix<Real, DimM, DimM>;
     using T4_t = T4Mat<Real, DimM>;
 
     using T2_map = Eigen::Map<T2_t>;
     using T4_map = T4MatMap<Real, DimM>;
 
     using T2_const_map = Eigen::Map<const T2_t>;
     using T4_const_map = T4MatMap<Real, DimM, true>;
 
     using FieldColl_t = GlobalFieldCollection<DimM>;
     using T2Field_t = TensorField<FieldColl_t, Real, secondOrder, DimM>;
     using T4Field_t = TensorField<FieldColl_t, Real, fourthOrder, DimM>;
 
     //! Default constructor
     MaterialEvaluator() = delete;
 
     /**
      * The constructor takes all arguments after the name of the underlying
      * Material's constructor. E.g., if the underlying material is a
      * `muSpectre::MaterialLinearElastic1<threeD>`, these would be Young's
      * modulus and Poisson's ratio.
      */
     template <class... ConstructorArgs>
     explicit MaterialEvaluator(ConstructorArgs... args)
         : material("name", std::forward<ConstructorArgs>(args)...),
           strain{make_field<T2Field_t>("gradient", this->collection)},
           stress{make_field<T2Field_t>("stress", this->collection)},
           tangent{make_field<T4Field_t>("tangent", this->collection)} {
       this->collection.initialise(CcoordOps::get_cube<DimM>(1), {0});
     }
 
     //! Copy constructor
     MaterialEvaluator(const MaterialEvaluator & other) = delete;
 
     //! Move constructor
     MaterialEvaluator(MaterialEvaluator && other) = default;
 
     //! Destructor
     virtual ~MaterialEvaluator() = default;
 
     //! Copy assignment operator
     MaterialEvaluator & operator=(const MaterialEvaluator & other) = delete;
 
     //! Move assignment operator
     MaterialEvaluator & operator=(MaterialEvaluator && other) = default;
 
     /**
      * Some materials have per-pixel info, which then must be given here
      * (i.e., those are the arguments to Material::add_pixel())
      */
     template <class... ArgTypes>
     void initialise(ArgTypes &&... args) {
       this->material.add_pixel({0}, std::forward<ArgTypes>(args)...);
       this->is_initialised = true;
     }
 
     /**
      * for materials with state variables. See `muSpectre::MaterialBase` for
      * details
      */
     void save_history_variables() { this->material.save_history_variables(); }
 
     /**
      * Evaluates the underlying materials constitutive law and returns the
      * stress P or σ as a function of the placement gradient F or small strain
      * tensor ε depending on the formulation
      * (`muSpectre::Formulation::small_strain` for σ(ε),
      * `muSpectre::Formulation::finite_strain` for P(F))
      */
     inline T2_const_map evaluate_stress(const Eigen::Ref<const T2_t> & grad,
                                         const Formulation & form);
 
     /**
      * Evaluates the underlying materials constitutive law and returns the the
      * stress P or σ and the tangent moduli K as a function of the placement
      * gradient F or small strain tensor ε depending on the formulation
      * (`muSpectre::Formulation::small_strain` for σ(ε),
      * `muSpectre::Formulation::finite_strain` for P(F))
      */
     inline std::tuple<T2_const_map, T4_const_map>
     evaluate_stress_tangent(const Eigen::Ref<const T2_t> & grad,
                             const Formulation & form);
 
+    /**
+     * estimate the tangent using finite difference
+     */
+    inline T4_t
+    estimate_tangent(const Eigen::Ref<const T2_t> & grad,
+                     const Formulation & form, const Real step,
+                     const FiniteDiff diff_type = FiniteDiff::centred);
+
    protected:
     /**
      * throws a runtime error if the material's per-pixel data has not been set.
      */
     void check_init() const;
 
     Material material;
     FieldColl_t collection{};
     T2Field_t & strain;
     T2Field_t & stress;
     T4Field_t & tangent;
     bool is_initialised{false};
 
    private:
   };
 
   /* ---------------------------------------------------------------------- */
   template <class Material>
   auto MaterialEvaluator<Material>::evaluate_stress(
       const Eigen::Ref<const T2_t> & grad, const Formulation & form)
       -> T2_const_map {
     this->check_init();
     this->strain.get_map()[0] = grad;
     this->material.compute_stresses(this->strain, this->stress, form);
     return T2_const_map(this->stress.get_map()[0].data());
   }
 
   /* ---------------------------------------------------------------------- */
   template <class Material>
   auto MaterialEvaluator<Material>::evaluate_stress_tangent(
       const Eigen::Ref<const T2_t> & grad, const Formulation & form)
       -> std::tuple<T2_const_map, T4_const_map> {
     this->check_init();
     this->strain.get_map()[0] = grad;
     this->material.compute_stresses_tangent(this->strain, this->stress,
                                             this->tangent, form);
     return std::make_tuple(T2_const_map(this->stress.get_map()[0].data()),
                            T4_const_map(this->tangent.get_map()[0].data()));
   }
 
   /* ---------------------------------------------------------------------- */
   template <class Material>
   void MaterialEvaluator<Material>::check_init() const {
     if (not this->is_initialised) {
       throw std::runtime_error("Not initialised");
     }
   }
+
+  /* ---------------------------------------------------------------------- */
+  template <class Material>
+  auto MaterialEvaluator<Material>::estimate_tangent(
+      const Eigen::Ref<const T2_t> & grad, const Formulation & form,
+      const Real delta, const FiniteDiff diff_type) -> T4_t {
+    using T2_vec = Eigen::Map<Eigen::Matrix<Real, DimM * DimM, 1>>;
+
+    T4_t tangent{T4_t::Zero()};
+
+    const T2_t stress{this->evaluate_stress(grad, form)};
+
+    auto fun{[&form, this](const auto & grad_) {
+      return this->evaluate_stress(grad_, form);
+    }};
+
+    auto symmetrise{
+        [](auto & eps) { eps = .5 * (eps + eps.transpose().eval()); }};
+
+    static_assert(Int(decltype(tangent.col(0))::SizeAtCompileTime) ==
+                      Int(T2_t::SizeAtCompileTime),
+                  "wrong column size");
+
+    for (Dim_t i{}; i < DimM * DimM; ++i) {
+      T2_t strain2{grad};
+      T2_vec strain2_vec{strain2.data()};
+      switch (diff_type) {
+      case FiniteDiff::forward: {
+        strain2_vec(i) += delta;
+        if (form == Formulation::small_strain) {
+          symmetrise(strain2);
+        }
+
+        T2_t del_f_del{(fun(strain2) - stress) / delta};
+
+        tangent.col(i) = T2_vec(del_f_del.data());
+        break;
+      }
+      case FiniteDiff::backward: {
+        strain2_vec(i) -= delta;
+        if (form == Formulation::small_strain) {
+          symmetrise(strain2);
+        }
+        T2_t del_f_del{(stress - fun(strain2)) / delta};
+
+        tangent.col(i) = T2_vec(del_f_del.data());
+        break;
+      }
+      case FiniteDiff::centred: {
+        T2_t strain1{grad};
+        T2_vec strain1_vec{strain1.data()};
+        strain1_vec(i) += delta;
+        strain2_vec(i) -= delta;
+        if (form == Formulation::small_strain) {
+          symmetrise(strain1);
+          symmetrise(strain2);
+        }
+        T2_t del_f_del{(fun(strain1).eval() - fun(strain2).eval()) /
+                       (2 * delta)};
+
+        tangent.col(i) = T2_vec(del_f_del.data());
+        break;
+      }
+      default: {
+        throw std::runtime_error("Unknown FiniteDiff value");
+        break;
+      }
+      }
+      static_assert(Int(decltype(tangent.col(i))::SizeAtCompileTime) ==
+                        Int(T2_t::SizeAtCompileTime),
+                    "wrong column size");
+    }
+    return tangent;
+  }
+
 }  // namespace muSpectre
 
 #endif  // SRC_MATERIALS_MATERIAL_EVALUATOR_HH_
diff --git a/tests/test_material_evaluator.cc b/tests/test_material_evaluator.cc
index 30b4bfa..6e7f397 100644
--- a/tests/test_material_evaluator.cc
+++ b/tests/test_material_evaluator.cc
@@ -1,189 +1,270 @@
 /**
  * @file   test_material_evaluator.cc
  *
  * @author Till Junge <till.junge@altermail.ch>
  *
  * @date   13 Jan 2019
  *
  * @brief  tests for the material evaluator mechanism
  *
  * Copyright © 2019 Till Junge
  *
  * µSpectre is free software; you can 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 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with µSpectre; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  *
  * Additional permission under GNU GPL version 3 section 7
  *
  * If you modify this Program, or any covered work, by linking or combining it
  * with proprietary FFT implementations or numerical libraries, containing parts
  * covered by the terms of those libraries' licenses, the licensors of this
  * Program grant you additional permission to convey the resulting work.
  */
 
 #include "tests.hh"
 #include "common/T4_map_proxy.hh"
 #include "materials/material_linear_elastic2.hh"
 #include "materials/material_evaluator.hh"
 
 #include "Eigen/Dense"
 
 namespace muSpectre {
 
   BOOST_AUTO_TEST_SUITE(material_evaluator_tests);
+
+  /* ---------------------------------------------------------------------- */
   BOOST_AUTO_TEST_CASE(without_per_pixel_data) {
     using Mat_t = MaterialLinearElastic1<twoD, twoD>;
     using Eval_t = MaterialEvaluator<Mat_t>;
 
     constexpr Real Young{210e9};
     constexpr Real Poisson{.33};
 
     Eval_t evaluator{Young, Poisson};
 
     using T2_t = Eigen::Matrix<Real, twoD, twoD>;
     using T4_t = T4Mat<Real, twoD>;
     const T2_t F{(T2_t::Random() - (T2_t::Ones() * .5)) * 1e-4 +
                  T2_t::Identity()};
     const T2_t eps{
         .5 * ((F - T2_t::Identity()) + (F - T2_t::Identity()).transpose())};
 
     BOOST_CHECK_THROW(evaluator.evaluate_stress(eps, Formulation::small_strain),
                       std::runtime_error);
 
     evaluator.initialise();
 
     const T2_t sigma{evaluator.evaluate_stress(eps, Formulation::small_strain)};
     const T2_t P{evaluator.evaluate_stress(F, Formulation::finite_strain)};
 
     auto J{F.determinant()};
     auto P_reconstruct{J * sigma * F.inverse().transpose()};
 
     auto error_comp{[](const auto & a, const auto & b) {
       return (a - b).norm() / (a + b).norm();
     }};
     auto error{error_comp(P, P_reconstruct)};
 
     constexpr Real small_strain_tol{1e-3};
-    if ( not (error <= small_strain_tol)) {
+    if (not(error <= small_strain_tol)) {
       std::cout << "F =" << std::endl << F << std::endl;
       std::cout << "ε =" << std::endl << eps << std::endl;
       std::cout << "P =" << std::endl << P << std::endl;
       std::cout << "σ =" << std::endl << sigma << std::endl;
       std::cout << "P_reconstructed =" << std::endl
                 << P_reconstruct << std::endl;
     }
 
     BOOST_CHECK_LE(error, small_strain_tol);
 
     T2_t sigma2, P2;
     T4_t C, K;
 
     std::tie(sigma2, C) =
       evaluator.evaluate_stress_tangent(eps, Formulation::small_strain);
     std::tie(P2, K) =
       evaluator.evaluate_stress_tangent(F, Formulation::finite_strain);
 
     error = error_comp(sigma2, sigma);
     BOOST_CHECK_LE(error, tol);
     error = error_comp(P2, P);
     BOOST_CHECK_LE(error, tol);
 
     error = error_comp(C, K);
-    if ( not (error <= small_strain_tol)) {
+    if (not(error <= small_strain_tol)) {
       std::cout << "F =" << std::endl << F << std::endl;
       std::cout << "ε =" << std::endl << eps << std::endl;
       std::cout << "P =" << std::endl << P << std::endl;
       std::cout << "σ =" << std::endl << sigma << std::endl;
       std::cout << "K =" << std::endl << K << std::endl;
       std::cout << "C =" << std::endl << C << std::endl;
     }
     BOOST_CHECK_LE(error, small_strain_tol);
   }
 
+  /* ---------------------------------------------------------------------- */
   BOOST_AUTO_TEST_CASE(with_per_pixel_data) {
     using Mat_t = MaterialLinearElastic2<twoD, twoD>;
     using Eval_t = MaterialEvaluator<Mat_t>;
 
     constexpr Real Young{210e9};
     constexpr Real Poisson{.33};
 
     Eval_t evaluator{Young, Poisson};
 
     using T2_t = Eigen::Matrix<Real, twoD, twoD>;
     using T4_t = T4Mat<Real, twoD>;
     const T2_t F{(T2_t::Random() - (T2_t::Ones() * .5)) * 1e-4 +
                  T2_t::Identity()};
     const T2_t eps{
         .5 * ((F - T2_t::Identity()) + (F - T2_t::Identity()).transpose())};
 
     BOOST_CHECK_THROW(evaluator.evaluate_stress(eps, Formulation::small_strain),
                       std::runtime_error);
 
     T2_t eigen_strain{[](auto x) {
       return 1e-4 * (x + x.transpose());
     }(T2_t::Random() - T2_t::Ones() * .5)};
 
     evaluator.initialise(eigen_strain);
 
     const T2_t sigma{evaluator.evaluate_stress(eps, Formulation::small_strain)};
     const T2_t P{evaluator.evaluate_stress(F, Formulation::finite_strain)};
 
     auto J{F.determinant()};
     auto P_reconstruct{J * sigma * F.inverse().transpose()};
 
     auto error_comp{[](const auto & a, const auto & b) {
       return (a - b).norm() / (a + b).norm();
     }};
     auto error{error_comp(P, P_reconstruct)};
 
     constexpr Real small_strain_tol{1e-3};
-    if ( not (error <= small_strain_tol)) {
+    if (not(error <= small_strain_tol)) {
       std::cout << "F =" << std::endl << F << std::endl;
       std::cout << "ε =" << std::endl << eps << std::endl;
       std::cout << "P =" << std::endl << P << std::endl;
       std::cout << "σ =" << std::endl << sigma << std::endl;
       std::cout << "P_reconstructed =" << std::endl
                 << P_reconstruct << std::endl;
     }
 
     BOOST_CHECK_LE(error, small_strain_tol);
 
     T2_t sigma2, P2;
     T4_t C, K;
 
     std::tie(sigma2, C) =
       evaluator.evaluate_stress_tangent(eps, Formulation::small_strain);
     std::tie(P2, K) =
       evaluator.evaluate_stress_tangent(F, Formulation::finite_strain);
 
     error = error_comp(sigma2, sigma);
     BOOST_CHECK_LE(error, tol);
     error = error_comp(P2, P);
     BOOST_CHECK_LE(error, tol);
 
     error = error_comp(C, K);
-    if ( not (error <= small_strain_tol)) {
+    if (not(error <= small_strain_tol)) {
       std::cout << "F =" << std::endl << F << std::endl;
       std::cout << "ε =" << std::endl << eps << std::endl;
       std::cout << "P =" << std::endl << P << std::endl;
       std::cout << "σ =" << std::endl << sigma << std::endl;
       std::cout << "K =" << std::endl << K << std::endl;
       std::cout << "C =" << std::endl << C << std::endl;
     }
     BOOST_CHECK_LE(error, small_strain_tol);
   }
 
+  /* ---------------------------------------------------------------------- */
+  BOOST_AUTO_TEST_CASE(tangent_estimation) {
+    using Mat_t = MaterialLinearElastic1<twoD, twoD>;
+    using Eval_t = MaterialEvaluator<Mat_t>;
+
+    constexpr Real Young{210e9};
+    constexpr Real Poisson{.33};
+
+    Eval_t evaluator{Young, Poisson};
+
+    using T2_t = Eigen::Matrix<Real, twoD, twoD>;
+    using T4_t = T4Mat<Real, twoD>;
+    const T2_t F{(T2_t::Random() - (T2_t::Ones() * .5)) * 1e-4 +
+                 T2_t::Identity()};
+    const T2_t eps{
+        .5 * ((F - T2_t::Identity()) + (F - T2_t::Identity()).transpose())};
+
+    BOOST_CHECK_THROW(evaluator.evaluate_stress(eps, Formulation::small_strain),
+                      std::runtime_error);
+
+    evaluator.initialise();
+
+    T2_t sigma, P;
+    T4_t C, K;
+
+    std::tie(sigma, C) =
+      evaluator.evaluate_stress_tangent(eps, Formulation::small_strain);
+    std::tie(P, K) =
+      evaluator.evaluate_stress_tangent(F, Formulation::finite_strain);
+
+    constexpr Real linear_step{1.};
+    constexpr Real nonlin_step{1.e-6};
+    T4_t C_estim{evaluator.estimate_tangent(eps, Formulation::small_strain,
+                                            linear_step)};
+    T4_t K_estim{evaluator.estimate_tangent(F, Formulation::finite_strain,
+                                            nonlin_step)};
+
+    auto error_comp{[](const auto & a, const auto & b) {
+      return (a - b).norm() / (a + b).norm();
+    }};
+
+    constexpr Real finite_diff_tol{1e-9};
+    Real error {error_comp(K, K_estim)};
+    if (not(error <= finite_diff_tol)) {
+      std::cout << "K =" << std::endl << K << std::endl;
+      std::cout << "K_estim =" << std::endl << K_estim << std::endl;
+    }
+    BOOST_CHECK_LE(error, finite_diff_tol);
+
+    error = error_comp(C, C_estim);
+    if (not(error <= tol)) {
+      std::cout << "centred difference:" << std::endl;
+      std::cout << "C =" << std::endl << C << std::endl;
+      std::cout << "C_estim =" << std::endl << C_estim << std::endl;
+    }
+    BOOST_CHECK_LE(error, tol);
+
+    C_estim = evaluator.estimate_tangent(eps, Formulation::small_strain,
+                                         linear_step, FiniteDiff::forward);
+    error = error_comp(C, C_estim);
+    if (not(error <= tol)) {
+      std::cout << "forward difference:" << std::endl;
+      std::cout << "C =" << std::endl << C << std::endl;
+      std::cout << "C_estim =" << std::endl << C_estim << std::endl;
+    }
+    BOOST_CHECK_LE(error, tol);
+
+    C_estim = evaluator.estimate_tangent(eps, Formulation::small_strain,
+                                         linear_step, FiniteDiff::backward);
+    error = error_comp(C, C_estim);
+    if (not(error <= tol)) {
+      std::cout << "backward difference:" << std::endl;
+      std::cout << "C =" << std::endl << C << std::endl;
+      std::cout << "C_estim =" << std::endl << C_estim << std::endl;
+    }
+    BOOST_CHECK_LE(error, tol);
+}
+
   BOOST_AUTO_TEST_SUITE_END();
 
 }  // namespace muSpectre