From 3261926bfaccc6fb7f8c70ecb41965a79e5068b1 Mon Sep 17 00:00:00 2001
From: Thomas Helfer <thomas.helfer@cea.fr>
Date: Tue, 25 May 2021 15:12:07 +0200
Subject: [PATCH] Synchronize with TFEL-3.4.1 release and later

---
 .../MGIS/Behaviour/MaterialDataManager.h      |  55 +++++
 bindings/c/src/Integrate.cxx                  |   1 -
 bindings/c/src/MaterialDataManager.cxx        | 139 ++++++++++-
 .../MGIS/FEniCS/FEniCSTestingUtilities.hxx    |   8 +-
 bindings/fortran/src/mgis_behaviour.f95       | 182 ++++++++++++++
 .../include/MGIS/Julia/JuliaUtilities.hxx     |  16 +-
 bindings/julia/src/BehaviourDataView.cxx      |   8 +
 bindings/julia/src/Integrate.cxx              |   1 -
 bindings/python/src/Integrate.cxx             | 105 ++++++---
 bindings/python/src/MaterialDataManager.cxx   |  42 +++-
 docs/web/mgis-template.html                   |  48 ++--
 include/MGIS/Behaviour/Behaviour.hxx          |   4 +-
 include/MGIS/Behaviour/BehaviourData.hxx      | 201 ++++++++--------
 include/MGIS/Behaviour/BehaviourDataView.hxx  |  15 +-
 include/MGIS/Behaviour/Integrate.hxx          | 223 +++++++++++++-----
 .../MGIS/Behaviour/MaterialDataManager.hxx    | 147 +++++++++++-
 include/MGIS/Config-c.h                       |  21 +-
 include/MGIS/LibrariesManager.hxx             |   9 +
 include/MGIS/ThreadedTaskResult.hxx           |  19 +-
 include/MGIS/ThreadedTaskResult.ixx           | 133 ++---------
 src/Behaviour.cxx                             |   9 +
 src/BehaviourData.cxx                         |   8 +-
 src/Integrate.cxx                             | 200 ++++++++++------
 src/LibrariesManager.cxx                      | 130 +++++-----
 src/MaterialDataManager.cxx                   | 141 +++++++++--
 src/ThreadedTaskResult.cxx                    |   2 +-
 tests/IntegrateTest.cxx                       |   2 +-
 27 files changed, 1360 insertions(+), 509 deletions(-)

diff --git a/bindings/c/include/MGIS/Behaviour/MaterialDataManager.h b/bindings/c/include/MGIS/Behaviour/MaterialDataManager.h
index 4aa86c500..90b5c9536 100644
--- a/bindings/c/include/MGIS/Behaviour/MaterialDataManager.h
+++ b/bindings/c/include/MGIS/Behaviour/MaterialDataManager.h
@@ -57,6 +57,15 @@ MGIS_C_EXPORT mgis_status mgis_bv_create_material_data_manager_initializer(
 MGIS_C_EXPORT
 mgis_status mgis_bv_material_data_manager_initializer_bind_tangent_operator(
     mgis_bv_MaterialDataManagerInitializer*, mgis_real* const, mgis_size_type);
+/*!
+ * \brief bind the speed of sound to the given array
+ * \param[in,out] d: initializer
+ * \param[in] p: pointer to a memory area meant to store the speed of sound values
+ * \param[in] s: size of the memory area
+ */
+MGIS_C_EXPORT
+mgis_status mgis_bv_material_data_manager_initializer_bind_speed_of_sound(
+    mgis_bv_MaterialDataManagerInitializer*, mgis_real* const, mgis_size_type);
 /*!
  * \brief set the state at the beginning of the time step
  * \param[out] s: pointer to a pointer to the state initializer
@@ -102,6 +111,52 @@ MGIS_C_EXPORT mgis_status mgis_bv_create_material_data_manager_with_initializer(
     const mgis_bv_Behaviour* const,
     const mgis_size_type,
     const mgis_bv_MaterialDataManagerInitializer* const);
+/*!
+ * \param[in]  d: material data manager
+ * \param[in]  b: boolean value
+ */
+MGIS_C_EXPORT mgis_status mgis_bv_material_data_manager_set_thread_safe(
+    mgis_bv_MaterialDataManager*, const mgis_size_type);
+/*!
+ * \param[in]  d: material data manager
+ */
+MGIS_C_EXPORT mgis_status
+mgis_bv_material_data_manager_allocate_array_of_tangent_operator_blocks(
+    mgis_bv_MaterialDataManager* const);
+/*!
+ * \param[in] d: material data manager
+ * \param[in] p: pointer to a memory buffer
+ * \param[in] s: size of the a memory buffer
+ */
+MGIS_C_EXPORT mgis_status
+mgis_bv_material_data_manager_use_external_array_of_tangent_operator_blocks(
+    mgis_bv_MaterialDataManager* const, mgis_real* const, const mgis_size_type);
+/*!
+ * \param[in]  d: material data manager
+ */
+MGIS_C_EXPORT mgis_status
+mgis_bv_material_data_manager_release_array_of_tangent_operator_blocks(
+    mgis_bv_MaterialDataManager*const);
+/*!
+ * \param[in]  d: material data manager
+ */
+MGIS_C_EXPORT mgis_status
+mgis_bv_material_data_manager_allocate_array_of_speed_of_sounds(
+    mgis_bv_MaterialDataManager*const);
+/*!
+ * \param[in] d: material data manager
+ * \param[in] p: pointer to a memory buffer
+ * \param[in] s: size of the a memory buffer
+ */
+MGIS_C_EXPORT mgis_status
+mgis_bv_material_data_manager_use_external_array_of_speed_of_sounds(
+    mgis_bv_MaterialDataManager* const, mgis_real* const, const mgis_size_type);
+/*!
+ * \param[in]  d: material data manager
+ */
+MGIS_C_EXPORT mgis_status
+mgis_bv_material_data_manager_release_array_of_speed_of_sounds(
+    mgis_bv_MaterialDataManager*const);
 /*!
  * \brief set the state at the beginning of the time step
  * \param[out] s: pointer to a pointer to the state
diff --git a/bindings/c/src/Integrate.cxx b/bindings/c/src/Integrate.cxx
index ab5f18ea4..6dbb3e4ad 100644
--- a/bindings/c/src/Integrate.cxx
+++ b/bindings/c/src/Integrate.cxx
@@ -60,7 +60,6 @@ mgis_status mgis_bv_integrate_2(int* const r,
                                 const mgis_bv_Behaviour* const b) {
   auto v = mgis::behaviour::make_view(*d);
   auto s = mgis_bv_integrate(r, &v, b);
-  d->rdt = v.rdt;
   return s;
 }  // end of mgis_bv_integrate
 
diff --git a/bindings/c/src/MaterialDataManager.cxx b/bindings/c/src/MaterialDataManager.cxx
index 0a2748fed..c576608e0 100644
--- a/bindings/c/src/MaterialDataManager.cxx
+++ b/bindings/c/src/MaterialDataManager.cxx
@@ -52,7 +52,29 @@ mgis_status mgis_bv_material_data_manager_initializer_bind_tangent_operator(
     return mgis_handle_cxx_exception();
   }
   return mgis_report_success();
-}  // end of mgis_bv_create_material_data_manager_initializer
+}  // end of mgis_bv_material_data_manager_initializer_bind_tangent_operator
+
+mgis_status mgis_bv_material_data_manager_initializer_bind_speed_of_sound(
+    mgis_bv_MaterialDataManagerInitializer* d,
+    mgis_real* const p,
+    mgis_size_type s) {
+  if (d == nullptr) {
+    return mgis_report_failure(
+        "mgis_bv_material_data_manager_initializer_bind_speed_of_sound: "
+        "null argument");
+  }
+  if (p == nullptr) {
+    return mgis_report_failure(
+        "mgis_bv_material_data_manager_initializer_bind_speed_of_sound: "
+        "invalid tangent operator");
+  }
+  try {
+    d->speed_of_sound = mgis::span<mgis::real>(p, s);
+  } catch (...) {
+    return mgis_handle_cxx_exception();
+  }
+  return mgis_report_success();
+}  // end of mgis_bv_create_material_data_manager_initializer_bind_speed_of_sound
 
 mgis_status mgis_bv_material_data_manager_initializer_get_state_0_initializer(
     mgis_bv_MaterialStateManagerInitializer** s,
@@ -135,6 +157,121 @@ mgis_status mgis_bv_create_material_data_manager_with_initializer(
   return mgis_report_success();
 }  // end of mgis_bv_create_material_data_manager_with_initializer
 
+mgis_status mgis_bv_material_data_manager_set_thread_safe(
+    mgis_bv_MaterialDataManager* const d, const mgis_size_type b) {
+  if (d == nullptr) {
+    return mgis_report_failure(
+        "mgis_bv_material_data_manager_set_thread_safe: "
+        "null behaviour");
+  }
+  try {
+    d->setThreadSafe(static_cast<bool>(b));
+  } catch (...) {
+    return mgis_handle_cxx_exception();
+  }
+  return mgis_report_success();
+}  // end of mgis_bv_material_data_manager_set_thread_safe
+
+mgis_status
+mgis_bv_material_data_manager_allocate_array_of_tangent_operator_blocks(
+    mgis_bv_MaterialDataManager* const d) {
+  if (d == nullptr) {
+    return mgis_report_failure(
+        "mgis_bv_material_data_manager_set_thread_safe: "
+        "null behaviour");
+  }
+  try {
+    d->allocateArrayOfTangentOperatorBlocks();
+  } catch (...) {
+    return mgis_handle_cxx_exception();
+  }
+  return mgis_report_success();
+}  // end of
+   // mgis_bv_material_data_manager_allocate_array_of_tangent_operator_blocks
+
+mgis_status
+mgis_bv_material_data_manager_use_external_array_of_tangent_operator_blocks(
+    mgis_bv_MaterialDataManager* const d,
+    mgis_real* const p,
+    const mgis_size_type n) {
+  if (d == nullptr) {
+    return mgis_report_failure(
+        "mgis_bv_material_data_manager_set_thread_safe: "
+        "null behaviour");
+  }
+  try {
+    d->useExternalArrayOfTangentOperatorBlocks(mgis::span<mgis::real>(p, n));
+  } catch (...) {
+    return mgis_handle_cxx_exception();
+  }
+  return mgis_report_success();
+}  // end of
+   // mgis_bv_material_data_manager_use_external_array_of_tangent_operator_blocks
+
+mgis_status
+mgis_bv_material_data_manager_release_array_of_tangent_operator_blocks(
+    mgis_bv_MaterialDataManager* const d) {
+  if (d == nullptr) {
+    return mgis_report_failure(
+        "mgis_bv_material_data_manager_set_thread_safe: "
+        "null behaviour");
+  }
+  try {
+    d->releaseArrayOfTangentOperatorBlocks();
+  } catch (...) {
+    return mgis_handle_cxx_exception();
+  }
+  return mgis_report_success();
+}  // end of
+   // mgis_bv_material_data_manager_release_array_of_tangent_operator_blocks
+
+mgis_status
+mgis_bv_material_data_manager_allocate_array_of_speed_of_sounds(
+    mgis_bv_MaterialDataManager* const d) {
+  if (d == nullptr) {
+    return mgis_report_failure(
+        "mgis_bv_material_data_manager_set_thread_safe: "
+        "null behaviour");
+  }
+  try {
+    d->allocateArrayOfSpeedOfSounds();
+  } catch (...) {
+    return mgis_handle_cxx_exception();
+  }
+  return mgis_report_success();
+}  // end of mgis_bv_material_data_manager_allocate_array_of_speed_of_sounds
+
+mgis_status mgis_bv_material_data_manager_use_external_array_of_speed_of_sounds(
+    mgis_bv_MaterialDataManager* const d,
+    mgis_real* const p,
+    const mgis_size_type n) {
+  if (d == nullptr) {
+    return mgis_report_failure(
+        "mgis_bv_material_data_manager_set_thread_safe: "
+        "null behaviour");
+  }
+  try {
+    d->useExternalArrayOfSpeedOfSounds(mgis::span<mgis::real>(p, n));
+  } catch (...) {
+    return mgis_handle_cxx_exception();
+  }
+  return mgis_report_success();
+}  // end of mgis_bv_material_data_manager_use_external_array_of_speed_of_sounds
+
+mgis_status mgis_bv_material_data_manager_release_array_of_speed_of_sounds(
+    mgis_bv_MaterialDataManager* const d) {
+  if (d == nullptr) {
+    return mgis_report_failure(
+        "mgis_bv_material_data_manager_set_thread_safe: "
+        "null behaviour");
+  }
+  try {
+    d->releaseArrayOfSpeedOfSounds();
+  } catch (...) {
+    return mgis_handle_cxx_exception();
+  }
+  return mgis_report_success();
+}  // end of mgis_bv_material_data_manager_release_array_of_speed_of_sounds
 mgis_status mgis_bv_material_data_manager_get_state_0(
     mgis_bv_MaterialStateManager** s, mgis_bv_MaterialDataManager* const d) {
   if (d == nullptr) {
diff --git a/bindings/fenics/tests/include/MGIS/FEniCS/FEniCSTestingUtilities.hxx b/bindings/fenics/tests/include/MGIS/FEniCS/FEniCSTestingUtilities.hxx
index ade9a6eb0..e0c4a5a96 100644
--- a/bindings/fenics/tests/include/MGIS/FEniCS/FEniCSTestingUtilities.hxx
+++ b/bindings/fenics/tests/include/MGIS/FEniCS/FEniCSTestingUtilities.hxx
@@ -52,8 +52,8 @@ namespace mgis {
      * - sy2: defined by the points satisfying y==1
      */
     MGIS_VISIBILITY_EXPORT
-        std::map<std::string, std::shared_ptr<dolfin::SubDomain>>
-        getUnitSquareBoundaries();
+    std::map<std::string, std::shared_ptr<dolfin::SubDomain>>
+    getUnitSquareBoundaries();
     /*!
      * \return all the boundaries associated with a unit cube
      * The following boundaries are defined:
@@ -65,8 +65,8 @@ namespace mgis {
      * - sz2: defined by the points satisfying z==1
      */
     MGIS_VISIBILITY_EXPORT
-        std::map<std::string, std::shared_ptr<dolfin::SubDomain>>
-        getUnitCubeBoundaries();
+    std::map<std::string, std::shared_ptr<dolfin::SubDomain>>
+    getUnitCubeBoundaries();
 
   }  // end of namespace fenics
 
diff --git a/bindings/fortran/src/mgis_behaviour.f95 b/bindings/fortran/src/mgis_behaviour.f95
index 2f852734b..149b3ee6d 100644
--- a/bindings/fortran/src/mgis_behaviour.f95
+++ b/bindings/fortran/src/mgis_behaviour.f95
@@ -2386,6 +2386,35 @@ end function material_data_manager_initializer_bind_tangent_operator_wrapper
     s = material_data_manager_initializer_bind_tangent_operator_wrapper(d%ptr,K_ptr,nc)
   end function material_data_manager_initializer_bind_tangent_operator
   !
+  function material_data_manager_initializer_bind_speed_of_sound(d,p,n) result(s)
+    use, intrinsic :: iso_c_binding, only: c_size_t, c_loc
+    use mgis_fortran_utilities
+    use mgis, only: mgis_status, report_failure
+    implicit none
+    interface
+       function material_data_manager_initializer_bind_speed_of_sound_wrapper(d,p,s) &
+            bind(c,name = 'mgis_bv_material_data_manager_initializer_bind_speed_of_sound') &
+            result(r)
+         use, intrinsic :: iso_c_binding, only: c_ptr, c_size_t
+         use mgis, only: mgis_status
+         implicit none
+         type(c_ptr), intent(in), value :: d
+         type(c_ptr), intent(in), value :: p
+         integer(kind=c_size_t), intent(in), value :: s
+         type(mgis_status) :: r
+       end function material_data_manager_initializer_bind_speed_of_sound_wrapper
+    end interface
+    type(MaterialDataManagerInitializer), intent(in) :: d
+    real(kind=8), dimension(:,:), target, intent(out) :: p
+    integer, intent(in) :: n
+    type(mgis_status) :: s
+    type(c_ptr) p_ptr
+    integer(kind=c_size_t) nc
+    nc=n
+    p_ptr = c_loc(p)
+    s = material_data_manager_initializer_bind_speed_of_sound_wrapper(d%ptr,p_ptr,nc)
+  end function material_data_manager_initializer_bind_speed_of_sound
+  !
   function free_material_data_manager_initializer(ptr) result(r)
     use, intrinsic :: iso_c_binding, only: c_associated
     use mgis
@@ -2469,6 +2498,159 @@ end function create_material_data_manager_with_initializer_wrapper
     s = create_material_data_manager_with_initializer_wrapper(d%ptr, b%ptr, nc, i%ptr)
   end function create_material_data_manager_with_initializer
   !
+  function material_data_manager_set_thread_safe(d, b) result(r)
+    use, intrinsic :: iso_c_binding, only: c_size_t
+    use mgis, only: mgis_status
+    implicit none
+    interface
+       function material_data_manager_set_thread_safe_wrapper(d, b) &
+            bind(c,name = 'mgis_bv_material_data_manager_set_thread_safe') result(r)
+         use, intrinsic :: iso_c_binding, only: c_ptr, c_size_t
+         use mgis, only: mgis_status
+         implicit none
+         type(c_ptr), intent(in), value :: d
+         integer(kind=c_size_t), intent(in), value :: b
+         type(mgis_status) :: r
+       end function material_data_manager_set_thread_safe_wrapper
+    end interface
+    type(MaterialDataManager), intent(in) :: d
+    integer(kind=c_size_t), intent(in) :: b
+    type(mgis_status) :: r
+    r = material_data_manager_set_thread_safe_wrapper(d%ptr, b)
+  end function material_data_manager_set_thread_safe
+  !
+  function material_data_manager_allocate_array_of_tangent_operator_blocks(d) result(r)
+    use mgis, only: mgis_status
+    implicit none
+    interface
+       function mdm_allocate_array_of_tangent_operator_blocks_wrapper(d) &
+            bind(c,name = 'mgis_bv_material_data_manager_allocate_array_of_tangent_operator_blocks') &
+            result(r)
+         use, intrinsic :: iso_c_binding, only: c_ptr
+         use mgis, only: mgis_status
+         implicit none
+         type(c_ptr), intent(in), value :: d
+         type(mgis_status) :: r
+       end function mdm_allocate_array_of_tangent_operator_blocks_wrapper
+    end interface
+    type(MaterialDataManager), intent(in) :: d
+    type(mgis_status) :: r
+    r = mdm_allocate_array_of_tangent_operator_blocks_wrapper(d%ptr)
+  end function material_data_manager_allocate_array_of_tangent_operator_blocks
+  !
+  function material_data_manager_use_array_of_tangent_operator_blocks(d, p , n) result(r)
+    use, intrinsic :: iso_c_binding, only: c_size_t, c_loc
+    use mgis, only: mgis_status
+    implicit none
+    interface
+       function mdm_use_array_of_tangent_operator_blocks_wrapper(d, p, n) &
+            bind(c,name = 'mgis_bv_material_data_manager_use_external_array_of_tangent_operator_blocks') &
+            result(r)
+         use, intrinsic :: iso_c_binding, only: c_ptr, c_size_t
+         use mgis, only: mgis_status
+         implicit none
+         type(c_ptr), intent(in), value :: d
+         type(c_ptr), intent(in), value :: p
+         integer(kind=c_size_t), intent(in), value :: n
+         type(mgis_status) :: r
+       end function mdm_use_array_of_tangent_operator_blocks_wrapper
+    end interface
+    type(MaterialDataManager), intent(in) :: d
+    real(kind=8), dimension(:,:), target, intent(out) :: p
+    integer, intent(in) :: n
+    type(mgis_status) :: r
+    type(c_ptr) p_ptr
+    integer(kind=c_size_t) nc
+    nc=n
+    p_ptr = c_loc(p)
+    r = mdm_use_array_of_tangent_operator_blocks_wrapper(d%ptr, p_ptr, nc)
+  end function material_data_manager_use_array_of_tangent_operator_blocks
+  !
+  function material_data_manager_release_array_of_tangent_operator_blocks(d) result(r)
+    use mgis, only: mgis_status
+    implicit none
+    interface
+       function mdm_release_array_of_tangent_operator_blocks_wrapper(d) &
+            bind(c,name = 'mgis_bv_material_data_manager_release_array_of_tangent_operator_blocks') &
+            result(r)
+         use, intrinsic :: iso_c_binding, only: c_ptr
+         use mgis, only: mgis_status
+         implicit none
+         type(c_ptr), intent(in), value :: d
+         type(mgis_status) :: r
+       end function mdm_release_array_of_tangent_operator_blocks_wrapper
+    end interface
+    type(MaterialDataManager), intent(in) :: d
+    type(mgis_status) :: r
+    r = mdm_release_array_of_tangent_operator_blocks_wrapper(d%ptr)
+  end function material_data_manager_release_array_of_tangent_operator_blocks
+  !
+  function material_data_manager_allocate_array_of_speed_of_sounds(d) result(r)
+    use mgis, only: mgis_status
+    implicit none
+    interface
+       function mdm_allocate_array_of_speed_of_sounds_wrapper(d) &
+            bind(c,name = 'mgis_bv_material_data_manager_allocate_array_of_speed_of_sounds') &
+            result(r)
+         use, intrinsic :: iso_c_binding, only: c_ptr
+         use mgis, only: mgis_status
+         implicit none
+         type(c_ptr), intent(in), value :: d
+         type(mgis_status) :: r
+       end function mdm_allocate_array_of_speed_of_sounds_wrapper
+    end interface
+    type(MaterialDataManager), intent(in) :: d
+    type(mgis_status) :: r
+    r = mdm_allocate_array_of_speed_of_sounds_wrapper(d%ptr)
+  end function material_data_manager_allocate_array_of_speed_of_sounds
+  !
+  function material_data_manager_use_array_of_speed_of_sounds(d, p , n) result(r)
+    use, intrinsic :: iso_c_binding, only: c_size_t, c_loc
+    use mgis, only: mgis_status
+    implicit none
+    interface
+       function mdm_use_array_of_speed_of_sounds_wrapper(d, p, n) &
+            bind(c,name = 'mgis_bv_material_data_manager_use_external_array_of_speed_of_sounds') &
+            result(r)
+         use, intrinsic :: iso_c_binding, only: c_ptr, c_size_t
+         use mgis, only: mgis_status
+         implicit none
+         type(c_ptr), intent(in), value :: d
+         type(c_ptr), intent(in), value :: p
+         integer(kind=c_size_t), intent(in), value :: n
+         type(mgis_status) :: r
+       end function mdm_use_array_of_speed_of_sounds_wrapper
+    end interface
+    type(MaterialDataManager), intent(in) :: d
+    real(kind=8), dimension(:,:), target, intent(out) :: p
+    integer, intent(in) :: n
+    type(mgis_status) :: r
+    type(c_ptr) p_ptr
+    integer(kind=c_size_t) nc
+    nc=n
+    p_ptr = c_loc(p)
+    r = mdm_use_array_of_speed_of_sounds_wrapper(d%ptr, p_ptr, nc)
+  end function material_data_manager_use_array_of_speed_of_sounds
+  !
+  function material_data_manager_release_array_of_speed_of_sounds(d) result(r)
+    use mgis, only: mgis_status
+    implicit none
+    interface
+       function mdm_release_array_of_speed_of_sounds_wrapper(d) &
+            bind(c,name = 'mgis_bv_material_data_manager_release_array_of_speed_of_sounds') &
+            result(r)
+         use, intrinsic :: iso_c_binding, only: c_ptr
+         use mgis, only: mgis_status
+         implicit none
+         type(c_ptr), intent(in), value :: d
+         type(mgis_status) :: r
+       end function mdm_release_array_of_speed_of_sounds_wrapper
+    end interface
+    type(MaterialDataManager), intent(in) :: d
+    type(mgis_status) :: r
+    r = mdm_release_array_of_speed_of_sounds_wrapper(d%ptr)
+  end function material_data_manager_release_array_of_speed_of_sounds
+  !
   function material_data_manager_get_state_0(s0,d) result(s)
     use mgis, only: mgis_status
     implicit none
diff --git a/bindings/julia/include/MGIS/Julia/JuliaUtilities.hxx b/bindings/julia/include/MGIS/Julia/JuliaUtilities.hxx
index c2336e71f..06689a13b 100644
--- a/bindings/julia/include/MGIS/Julia/JuliaUtilities.hxx
+++ b/bindings/julia/include/MGIS/Julia/JuliaUtilities.hxx
@@ -29,30 +29,30 @@ namespace jlcxx {
 
   //!
   template <>
-  struct IsBits<mgis::behaviour::Variable::Type> : std::true_type {};
+  struct IsMirroredType<mgis::behaviour::Variable::Type> : std::true_type {};
   //!
   template <>
-  struct IsBits<mgis::behaviour::Hypothesis> : std::true_type {};
+  struct IsMirroredType<mgis::behaviour::Hypothesis> : std::true_type {};
   //!
   template <>
-  struct IsBits<mgis::behaviour::Behaviour::Symmetry> : std::true_type {};
+  struct IsMirroredType<mgis::behaviour::Behaviour::Symmetry> : std::true_type {};
   //!
   template <>
-  struct IsBits<mgis::behaviour::Behaviour::BehaviourType> : std::true_type {};
+  struct IsMirroredType<mgis::behaviour::Behaviour::BehaviourType> : std::true_type {};
   //!
   template <>
-  struct IsBits<mgis::behaviour::Behaviour::Kinematic> : std::true_type {};
+  struct IsMirroredType<mgis::behaviour::Behaviour::Kinematic> : std::true_type {};
   //!
   template <>
-  struct IsBits<mgis::behaviour::FiniteStrainBehaviourOptions::StressMeasure>
+  struct IsMirroredType<mgis::behaviour::FiniteStrainBehaviourOptions::StressMeasure>
       : std::true_type {};
   //!
   template <>
-  struct IsBits<mgis::behaviour::FiniteStrainBehaviourOptions::TangentOperator>
+  struct IsMirroredType<mgis::behaviour::FiniteStrainBehaviourOptions::TangentOperator>
       : std::true_type {};
   //! \brief
   template <>
-  struct IsBits<mgis::behaviour::IntegrationType> : std::true_type {};
+  struct IsMirroredType<mgis::behaviour::IntegrationType> : std::true_type {};
 }  // end of namespace jlcxx
 
 namespace mgis {
diff --git a/bindings/julia/src/BehaviourDataView.cxx b/bindings/julia/src/BehaviourDataView.cxx
index 5a28bec9c..d02d37a6c 100644
--- a/bindings/julia/src/BehaviourDataView.cxx
+++ b/bindings/julia/src/BehaviourDataView.cxx
@@ -16,6 +16,14 @@
 #include "MGIS/Behaviour/BehaviourDataView.hxx"
 #include "MGIS/Julia/JuliaUtilities.hxx"
 
+namespace jlcxx {
+
+  //! \brief
+  template <>
+  struct IsMirroredType<mgis::behaviour::BehaviourDataView> : std::false_type {
+  };
+
+}  // end of namespace jlcxx
 void declareBehaviourDataView();
 
 void declareBehaviourDataView(jlcxx::Module& m) {
diff --git a/bindings/julia/src/Integrate.cxx b/bindings/julia/src/Integrate.cxx
index a2d98de6f..037a79f72 100644
--- a/bindings/julia/src/Integrate.cxx
+++ b/bindings/julia/src/Integrate.cxx
@@ -23,7 +23,6 @@ static int integrateBehaviourData1(mgis::behaviour::BehaviourData& d,
                                    const mgis::behaviour::Behaviour& b) {
   auto v = mgis::behaviour::make_view(d);
   const auto s = mgis::behaviour::integrate(v, b);
-  d.rdt = v.rdt;
   return s;
 }  // end of integrateBehaviourData
 
diff --git a/bindings/python/src/Integrate.cxx b/bindings/python/src/Integrate.cxx
index 8d7bda106..a263180cf 100644
--- a/bindings/python/src/Integrate.cxx
+++ b/bindings/python/src/Integrate.cxx
@@ -14,10 +14,12 @@
 
 #include <boost/python/def.hpp>
 #include <boost/python/enum.hpp>
+#include <boost/python/class.hpp>
 #include "MGIS/ThreadPool.hxx"
 #include "MGIS/Behaviour/BehaviourData.hxx"
 #include "MGIS/Behaviour/MaterialDataManager.hxx"
 #include "MGIS/Behaviour/Integrate.hxx"
+#include "MGIS/Python/VectorConverter.hxx"
 
 void declareIntegrate();
 
@@ -25,61 +27,104 @@ static int integrateBehaviourData1(mgis::behaviour::BehaviourData& d,
                                    const mgis::behaviour::Behaviour& b) {
   auto v = mgis::behaviour::make_view(d);
   const auto s = mgis::behaviour::integrate(v, b);
-  d.rdt = v.rdt;
   return s;
 }  // end of integrateBehaviourData
 
 void declareIntegrate() {
-  boost::python::enum_<mgis::behaviour::IntegrationType>("IntegrationType")
+  using namespace mgis::behaviour;
+
+  boost::python::enum_<IntegrationType>("IntegrationType")
       .value("PREDICTION_TANGENT_OPERATOR",
-             mgis::behaviour::IntegrationType::PREDICTION_TANGENT_OPERATOR)
+             IntegrationType::PREDICTION_TANGENT_OPERATOR)
       .value("PredictionWithTangentOperator",
-             mgis::behaviour::IntegrationType::PREDICTION_TANGENT_OPERATOR)
+             IntegrationType::PREDICTION_TANGENT_OPERATOR)
       .value("PREDICTION_SECANT_OPERATOR",
-             mgis::behaviour::IntegrationType::PREDICTION_SECANT_OPERATOR)
+             IntegrationType::PREDICTION_SECANT_OPERATOR)
       .value("PredictionWithSecantOperator",
-             mgis::behaviour::IntegrationType::PREDICTION_SECANT_OPERATOR)
+             IntegrationType::PREDICTION_SECANT_OPERATOR)
       .value("PREDICTION_ELASTIC_OPERATOR",
-             mgis::behaviour::IntegrationType::PREDICTION_ELASTIC_OPERATOR)
+             IntegrationType::PREDICTION_ELASTIC_OPERATOR)
       .value("PredictionWithElasticOperator",
-             mgis::behaviour::IntegrationType::PREDICTION_ELASTIC_OPERATOR)
+             IntegrationType::PREDICTION_ELASTIC_OPERATOR)
       .value("INTEGRATION_NO_TANGENT_OPERATOR",
-             mgis::behaviour::IntegrationType::INTEGRATION_NO_TANGENT_OPERATOR)
+             IntegrationType::INTEGRATION_NO_TANGENT_OPERATOR)
       .value("IntegrationWithoutTangentOperator",
-             mgis::behaviour::IntegrationType::INTEGRATION_NO_TANGENT_OPERATOR)
+             IntegrationType::INTEGRATION_NO_TANGENT_OPERATOR)
       .value("INTEGRATION_ELASTIC_OPERATOR",
-             mgis::behaviour::IntegrationType::INTEGRATION_ELASTIC_OPERATOR)
+             IntegrationType::INTEGRATION_ELASTIC_OPERATOR)
       .value("IntegrationWithElasticOperator",
-             mgis::behaviour::IntegrationType::INTEGRATION_ELASTIC_OPERATOR)
+             IntegrationType::INTEGRATION_ELASTIC_OPERATOR)
       .value("INTEGRATION_SECANT_OPERATOR",
-             mgis::behaviour::IntegrationType::INTEGRATION_SECANT_OPERATOR)
+             IntegrationType::INTEGRATION_SECANT_OPERATOR)
       .value("IntegrationWithSecantOperator",
-             mgis::behaviour::IntegrationType::INTEGRATION_SECANT_OPERATOR)
+             IntegrationType::INTEGRATION_SECANT_OPERATOR)
       .value("INTEGRATION_TANGENT_OPERATOR",
-             mgis::behaviour::IntegrationType::INTEGRATION_TANGENT_OPERATOR)
+             IntegrationType::INTEGRATION_TANGENT_OPERATOR)
       .value("IntegrationWithTangentOperator",
-             mgis::behaviour::IntegrationType::INTEGRATION_TANGENT_OPERATOR)
+             IntegrationType::INTEGRATION_TANGENT_OPERATOR)
       .value("INTEGRATION_CONSISTENT_TANGENT_OPERATOR",
-             mgis::behaviour::IntegrationType::
-                 INTEGRATION_CONSISTENT_TANGENT_OPERATOR)
+             IntegrationType::INTEGRATION_CONSISTENT_TANGENT_OPERATOR)
       .value("IntegrationWithConsistentTangentOperator",
-             mgis::behaviour::IntegrationType::
-                 INTEGRATION_CONSISTENT_TANGENT_OPERATOR);
+             IntegrationType::INTEGRATION_CONSISTENT_TANGENT_OPERATOR);
+
+  boost::python::class_<BehaviourIntegrationOptions>(
+      "BehaviourIntegrationOptions")
+      .add_property("integration_type",
+                    &BehaviourIntegrationOptions::integration_type)
+      .add_property("compute_speed_of_sound",
+                    &BehaviourIntegrationOptions::compute_speed_of_sound);
+
+  boost::python::class_<BehaviourIntegrationResult>(
+      "BehaviourIntegrationResult")
+      .add_property("exit_status", &BehaviourIntegrationResult::exit_status,
+                    "The returned value has the following meaning:\n"
+                    "- -1: integration failed for at least one Gauss point\n"
+                    "-  0: all integrations succeeded but results are\n "
+                    "      unreliable for at least one Gauss point\n"
+                    "-  1: integration succeeded and results are reliable.")
+      .add_property("time_step_increase_factor",
+                    &BehaviourIntegrationResult::time_step_increase_factor)
+      .add_property("n", &BehaviourIntegrationResult::n,
+                    "number of the integration point that failed\n"
+                    "or number of the last integration point \n"
+                    "that reported unreliable results")
+      .add_property("error_message",
+                    &BehaviourIntegrationResult::error_message);
+
+  // wrapping std::vector<BehaviourIntegrationResult>
+  mgis::python::initializeVectorConverter<
+      std::vector<BehaviourIntegrationResult>>();
+
+  boost::python::class_<MultiThreadedBehaviourIntegrationResult>(
+      "MultiThreadedBehaviourIntegrationResult")
+      .add_property("exit_status",
+                    &MultiThreadedBehaviourIntegrationResult::exit_status,
+                    "The returned value has the following meaning:\n"
+                    "- -1: integration failed for at least one Gauss point\n"
+                    "-  0: all integrations succeeded but results are\n "
+                    "      unreliable for at least one Gauss point\n"
+                    "-  1: integration succeeded and results are reliable.")
+      .add_property("results",
+                    &MultiThreadedBehaviourIntegrationResult::results);
 
-  int (*integrate_ptr1)(mgis::behaviour::BehaviourDataView&,
-                        const mgis::behaviour::Behaviour&) =
-      mgis::behaviour::integrate;
-  int (*integrate_ptr2)(mgis::behaviour::MaterialDataManager&,
-                        const mgis::behaviour::IntegrationType,
+  int (*integrate_ptr1)(BehaviourDataView&, const Behaviour&) = integrate;
+  int (*integrate_ptr2)(MaterialDataManager&, const IntegrationType,
                         const mgis::real, const mgis::size_type,
-                        const mgis::size_type) = mgis::behaviour::integrate;
-  int (*integrate_ptr3)(mgis::ThreadPool&,
-                        mgis::behaviour::MaterialDataManager&,
-                        const mgis::behaviour::IntegrationType,
-                        const mgis::real) = mgis::behaviour::integrate;
+                        const mgis::size_type) = integrate;
+  int (*integrate_ptr3)(mgis::ThreadPool&, MaterialDataManager&,
+                        const IntegrationType, const mgis::real) = integrate;
+  BehaviourIntegrationResult (*integrate_ptr4)(
+      MaterialDataManager&, const BehaviourIntegrationOptions&,
+      const mgis::real, const mgis::size_type, const mgis::size_type) =
+      integrate;
+  MultiThreadedBehaviourIntegrationResult (*integrate_ptr5)(
+      mgis::ThreadPool&, MaterialDataManager&,
+      const BehaviourIntegrationOptions&, const mgis::real) = integrate;
 
   boost::python::def("integrate", &integrateBehaviourData1);
   boost::python::def("integrate", integrate_ptr1);
   boost::python::def("integrate", integrate_ptr2);
   boost::python::def("integrate", integrate_ptr3);
+  boost::python::def("integrate", integrate_ptr4);
+  boost::python::def("integrate", integrate_ptr5);
 }  // end of declareIntegrate
diff --git a/bindings/python/src/MaterialDataManager.cxx b/bindings/python/src/MaterialDataManager.cxx
index b5ed45bdb..cc8248a4a 100644
--- a/bindings/python/src/MaterialDataManager.cxx
+++ b/bindings/python/src/MaterialDataManager.cxx
@@ -26,6 +26,22 @@ static void MaterialDataManagerInitializer_bindTangentOperator(
   i.K = mgis::python::mgis_convert_to_span(K);
 }  // end of MaterialDataManagerInitializer_bindTangentOperator
 
+static void MaterialDataManagerInitializer_bindSpeedOfSound(
+    mgis::behaviour::MaterialDataManagerInitializer& i,
+    boost::python::object vs) {
+  i.speed_of_sound = mgis::python::mgis_convert_to_span(vs);
+}  // end of MaterialDataManagerInitializer_bindSpeedOfSound
+
+static void MaterialDataManager_useExternalArrayOfTangentOperatorBlocks(
+    mgis::behaviour::MaterialDataManager& m, boost::python::object K) {
+  m.useExternalArrayOfTangentOperatorBlocks(
+      mgis::python::mgis_convert_to_span(K));
+}  // end of MaterialDataManager_useExternalArrayOfTangentOperatorBlocks
+
+static void MaterialDataManager_useExternalArrayOfSpeedOfSounds(
+    mgis::behaviour::MaterialDataManager& m, boost::python::object vs) {
+  m.useExternalArrayOfSpeedOfSounds(mgis::python::mgis_convert_to_span(vs));
+}  // end of MaterialDataManager_useExternalArrayOfSpeedOfSounds
 static boost::python::object MaterialDataManager_getK(
     mgis::behaviour::MaterialDataManager& d) {
   if (d.b.to_blocks.size() == 1u) {
@@ -54,13 +70,37 @@ void declareMaterialDataManager() {
       .add_property("s1", &MaterialDataManagerInitializer::s1)
       .def("bindTangentOperator",
            &MaterialDataManagerInitializer_bindTangentOperator,
-           "use the given array to store the tangent operator");
+           "use the given array to store the tangent operator blocks")
+      .def("bindSpeedOfSound", &MaterialDataManagerInitializer_bindSpeedOfSound,
+           "use the given array to store the speed of sounds");
   // exporting the MaterialDataManager class
   boost::python::class_<MaterialDataManager, boost::noncopyable>(
       "MaterialDataManager",
       boost::python::init<const Behaviour&, const size_type>())
       .def(boost::python::init<const Behaviour&, const size_type,
                                const MaterialDataManagerInitializer&>())
+      .def("setThreadSafe", &MaterialDataManager::setThreadSafe,
+           "specify if various operations performed by the MaterialDataManager "
+           "(memory allocations) shall be protected by a mutex")
+      .def("allocateArrayOfTangentOperatorBlocks",
+           &MaterialDataManager::allocateArrayOfTangentOperatorBlocks,
+           "allocate the tangent operator blocks")
+      .def("useExternalArrayOfTangentOperatorBlocks",
+           &MaterialDataManager_useExternalArrayOfTangentOperatorBlocks,
+           "use and externally defined array to store the tangent operator "
+           "blocks")
+      .def("releaseArrayOfTangentOperatorBlocks",
+           &MaterialDataManager::releaseArrayOfTangentOperatorBlocks,
+           "release the arrays of tangent operator blocks")
+      .def("allocateArrayOfSpeedOfSounds",
+           &MaterialDataManager::allocateArrayOfSpeedOfSounds,
+           "allocate the array of speed of sounds")
+      .def("useExternalArrayOfSpeedOfSounds",
+           &MaterialDataManager_useExternalArrayOfSpeedOfSounds,
+           "use and externally defined array to store the spped of sounds")
+      .def("releaseArrayOfSpeedOfSounds",
+           &MaterialDataManager::releaseArrayOfSpeedOfSounds,
+           "release the array of speed of sounds")
       .def_readonly("n", &MaterialDataManager::n)
       .def_readonly("number_of_integration_points", &MaterialDataManager::n)
       .add_property("s0", &MaterialDataManager::s0)
diff --git a/docs/web/mgis-template.html b/docs/web/mgis-template.html
index 12851f5f2..580ca3f8f 100644
--- a/docs/web/mgis-template.html
+++ b/docs/web/mgis-template.html
@@ -102,28 +102,28 @@
       </ul>
     </header>
     <br></br>
-    $if(title)$
-    <div id="$idprefix$header">
-      <h1 class="title">$title$</h1>
-      $if(subtitle)$
-      <h1 class="subtitle">$subtitle$</h1>
-      $endif$
-      $for(author)$
-      <h2 class="author">$author$</h2>
-      $endfor$
-      $if(date)$
-      <h3 class="date">$date$</h3>
-      $endif$
-    </div>
-    $endif$
-    $if(toc)$
-    <div id="$idprefix$TOC">
-      $toc$
-    </div>
-    $endif$
-    $body$
-    $for(include-after)$
-    $include-after$
-    $endfor$
-  </body>
+$if(title)$
+<div id="$idprefix$header">
+<h1 class="title">$title$</h1>
+$if(subtitle)$
+<h1 class="subtitle">$subtitle$</h1>
+$endif$
+$for(author)$
+<h2 class="author">$author$</h2>
+$endfor$
+$if(date)$
+<h3 class="date">$date$</h3>
+$endif$
+</div>
+$endif$
+$if(toc)$
+<div id="$idprefix$TOC">
+$toc$
+</div>
+$endif$
+$body$
+$for(include-after)$
+$include-after$
+$endfor$
+</body>
 </html>
diff --git a/include/MGIS/Behaviour/Behaviour.hxx b/include/MGIS/Behaviour/Behaviour.hxx
index 5b04966e3..e2b01f1e6 100644
--- a/include/MGIS/Behaviour/Behaviour.hxx
+++ b/include/MGIS/Behaviour/Behaviour.hxx
@@ -28,9 +28,7 @@ namespace mgis {
 
   namespace behaviour {
 
-    /*!
-     * \brief structure describing a behaviour
-     */
+    //! \brief structure describing a behaviour
     struct MGIS_EXPORT Behaviour {
       /*!
        * \brief maximum number of behaviour options, whatever the kind of
diff --git a/include/MGIS/Behaviour/BehaviourData.hxx b/include/MGIS/Behaviour/BehaviourData.hxx
index db2196c54..7b514eff2 100644
--- a/include/MGIS/Behaviour/BehaviourData.hxx
+++ b/include/MGIS/Behaviour/BehaviourData.hxx
@@ -21,111 +21,120 @@
 #include "MGIS/Behaviour/State.hxx"
 #include "MGIS/Behaviour/BehaviourDataView.hxx"
 
-namespace mgis {
+namespace mgis::behaviour {
 
-  namespace behaviour {
-
-    // forward declaration
-    struct Behaviour;
-
-    /*!
-     * \brief structure in charge of containing the data required for a
-     * behaviour integration.
-     */
-    struct MGIS_EXPORT BehaviourData {
-      /*!
-       * \brief constructor from a behaviour
-       * \param[in] b: behaviour
-       */
-      BehaviourData(const Behaviour&);
-      //! \brief move constructor
-      BehaviourData(BehaviourData&&);
-      //! \brief copy constructor
-      BehaviourData(const BehaviourData&);
-      //! \brief move assignement
-      BehaviourData& operator=(BehaviourData&&);
-      //! \brief copy assignement
-      BehaviourData& operator=(const BehaviourData&);
-      //! \brief time increment
-      mgis::real dt;
-      /*!
-       * \brief the stiffness matrix.
-       *
-       * On input, the first element of K (K[0]) must contain the type of type
-       * of computation to be performed.
-       *
-       * Let Ke be equal to:
-       *
-       * - K[0] - 100 if K[0] is greater than 50
-       * - K[0] otherwise.
-       *
-       * If Ke is negative, only the prediction operator is computed and
-       * no behaviour integration is performed.
-       *
-       * Ke has the following meaning:
-       *
-       * - if Ke is lower than -2.5, the tangent operator must be
-       *   computed.
-       * - if Ke is in [-2.5:-1.5]: the secant operator must be
-       *   computed.
-       * - if Ke is in [-1.5:-0.5]: the elastic operator must be
-       *   computed.
-       * - if Ke is in [-0.5:0.5]: the behaviour integration is
-       *   performed, but no stiffness matrix.
-       * - if Ke is in [0.5:1.5]: the elastic operator must be
-       *   computed.
-       * - if Ke is in [1.5:2.5]: the secant operator must be
-       *   computed.
-       * - if Ke is in [2.5:3.5]: the secant operator must be
-       *   computed.
-       * - if Ke is greater than 3.5, the consistent tangent operator
-       *   must be computed.
-       */
-      std::vector<real> K;
-      //! \brief proposed time step increment increase factor
-      real rdt;
-      //! \brief speed of sound (only computed if requested)
-      real speed_of_sound = 0;
-      //! \brief state at the beginning of the time step
-      State s0;
-      //! \brief state at the end of the time step
-      State s1;
-    };  // end of struct BehaviourData
+  // forward declaration
+  struct Behaviour;
 
+  /*!
+   * \brief structure in charge of containing the data required for a
+   * behaviour integration.
+   */
+  struct MGIS_EXPORT BehaviourData {
     /*!
-     * \brief update the behaviour data by:
-     * - setting s1 equal to s0
-     * - filling the stiffness matrix with 0
-     */
-    MGIS_EXPORT void update(BehaviourData&);
-    /*!
-     * \brief revert the behaviour data by:
-     * - setting s1 equal to s0
-     * - filling the stiffness matrix with 0
+     * \brief constructor from a behaviour
+     * \param[in] b: behaviour
      */
-    MGIS_EXPORT void revert(BehaviourData&);
+    BehaviourData(const Behaviour&);
+    //! \brief move constructor
+    BehaviourData(BehaviourData&&);
+    //! \brief copy constructor
+    BehaviourData(const BehaviourData&);
+    //! \brief move assignement
+    BehaviourData& operator=(BehaviourData&&);
+    //! \brief copy assignement
+    BehaviourData& operator=(const BehaviourData&);
     /*!
-     * \brief make a view from a behaviour data
-     * \param[in] d: data
-     * \return the view
-     * \note the view has at most the same life time as the data.
+     * \brief a pointer to a buffer used to store error message
+     *
+     * By default, this is initialised to an internal buffer. This is
+     * **not** thread-safe.
+     *
+     * The user may use its own buffer by setting this pointer appropriately
+     * or set it to the null pointer. If not null, the pointer must
+     * point to a buffer which is at least 512 characters wide (longer
+     * error message are truncated).  The user must ensure
+     * thread-safety (i.e. each thread shall have its own buffer).
      */
-    MGIS_EXPORT BehaviourDataView make_view(BehaviourData&);
+    char* error_message;
+    //! \brief time increment
+    mgis::real dt;
     /*!
-     * \brief print a detailled (verbose) description of the data associated
-     * with an integration point using a markdown format
-     * \param[in] os: ouptut stream
-     * \param[in] b: behaviour
-     * \param[in] d: behaviour data
-     * \param[in] l: title level
+     * \brief the stiffness matrix.
+     *
+     * On input, the first element of K (K[0]) must contain the type of type
+     * of computation to be performed.
+     *
+     * Let Ke be equal to:
+     *
+     * - K[0] - 100 if K[0] is greater than 50
+     * - K[0] otherwise.
+     *
+     * If Ke is negative, only the prediction operator is computed and
+     * no behaviour integration is performed.
+     *
+     * Ke has the following meaning:
+     *
+     * - if Ke is lower than -2.5, the tangent operator must be
+     *   computed.
+     * - if Ke is in [-2.5:-1.5]: the secant operator must be
+     *   computed.
+     * - if Ke is in [-1.5:-0.5]: the elastic operator must be
+     *   computed.
+     * - if Ke is in [-0.5:0.5]: the behaviour integration is
+     *   performed, but no stiffness matrix.
+     * - if Ke is in [0.5:1.5]: the elastic operator must be
+     *   computed.
+     * - if Ke is in [1.5:2.5]: the secant operator must be
+     *   computed.
+     * - if Ke is in [2.5:3.5]: the secant operator must be
+     *   computed.
+     * - if Ke is greater than 3.5, the consistent tangent operator
+     *   must be computed.
      */
-    MGIS_EXPORT void print_markdown(std::ostream&,
-                                    const Behaviour&,
-                                    const BehaviourData&,
-                                    const mgis::size_type);
+    std::vector<real> K;
+    //! \brief proposed time step increment increase factor
+    real rdt;
+    //! \brief speed of sound (only computed if requested)
+    real speed_of_sound = 0;
+    //! \brief state at the beginning of the time step
+    State s0;
+    //! \brief state at the end of the time step
+    State s1;
+  };  // end of struct BehaviourData
 
-  }  // end of namespace behaviour
+  /*!
+   * \brief update the behaviour data by:
+   * - setting s1 equal to s0
+   * - filling the stiffness matrix with 0
+   */
+  MGIS_EXPORT void update(BehaviourData&);
+  /*!
+   * \brief revert the behaviour data by:
+   * - setting s1 equal to s0
+   * - filling the stiffness matrix with 0
+   */
+  MGIS_EXPORT void revert(BehaviourData&);
+  /*!
+   * \brief make a view from a behaviour data
+   * \param[in] d: data
+   * \return the view
+   * \note the view has at most the same life time as the data.
+   */
+  MGIS_EXPORT BehaviourDataView make_view(BehaviourData&);
+  /*!
+   * \brief print a detailled (verbose) description of the data associated
+   * with an integration point using a markdown format
+   * \param[in] os: ouptut stream
+   * \param[in] b: behaviour
+   * \param[in] d: behaviour data
+   * \param[in] l: title level
+   */
+  MGIS_EXPORT void print_markdown(std::ostream&,
+                                  const Behaviour&,
+                                  const BehaviourData&,
+                                  const mgis::size_type);
 
-}  // end of namespace mgis
+}  // end of namespace mgis::behaviour
 
 #endif /* LIB_MGIS_BEHAVIOUR_BEHAVIOURDATA_HXX */
diff --git a/include/MGIS/Behaviour/BehaviourDataView.hxx b/include/MGIS/Behaviour/BehaviourDataView.hxx
index 4d65c1a4b..1fa78a36b 100644
--- a/include/MGIS/Behaviour/BehaviourDataView.hxx
+++ b/include/MGIS/Behaviour/BehaviourDataView.hxx
@@ -32,6 +32,15 @@ typedef struct mgis_bv_BehaviourDataView mgis_bv_BehaviourDataView;
 
 //! \brief view on the data associated with an integration point
 struct mgis_bv_BehaviourDataView {
+  /*!
+   * \brief pointer to a buffer used to store error message
+   *
+   * \note This pointer can be nullptr.  If not null, the pointer must
+   * point to a buffer which is at least 512 characters wide (longer
+   * error message are truncated).  This caller must ensure
+   * thread-safety (i.e. each thread shall have its own buffer).
+   */
+  char* error_message;
   //! \brief time increment
   mgis_real dt;
   /*!
@@ -40,6 +49,8 @@ struct mgis_bv_BehaviourDataView {
    * On input, the first element of K (K[0]) must contain the type of type
    * of computation to be performed.
    *
+   * If K[0] is greater than 50, the speed of sound is computed.
+   *
    * Let Ke be equal to:
    *
    * - K[0] - 100 if K[0] is greater than 50
@@ -88,9 +99,9 @@ struct mgis_bv_BehaviourDataView {
    */
   mgis_real* K;
   //! \brief proposed time step increment increase factor
-  mgis_real rdt;
+  mgis_real* rdt;
   //! \brief speed of sound (only computed if requested)
-  mgis_real speed_of_sound;
+  mgis_real* speed_of_sound;
   //! \brief state at the beginning of the time step
   mgis_bv_InitialStateView s0;
   //! \brief state at the end of the time step
diff --git a/include/MGIS/Behaviour/Integrate.hxx b/include/MGIS/Behaviour/Integrate.hxx
index 192f2e2fd..fdfff2caa 100644
--- a/include/MGIS/Behaviour/Integrate.hxx
+++ b/include/MGIS/Behaviour/Integrate.hxx
@@ -14,6 +14,8 @@
 #ifndef LIB_MGIS_BEHAVIOUR_INTEGRATE_HXX
 #define LIB_MGIS_BEHAVIOUR_INTEGRATE_HXX
 
+#include <limits>
+#include <thread>
 #include <vector>
 #include "MGIS/Config.hxx"
 #include "MGIS/Behaviour/BehaviourDataView.hxx"
@@ -45,69 +47,170 @@ namespace mgis {
     };  // end of enum IntegrationType
 
     /*!
-     * \brief structure in charge of handling temporary memory access.
-     */
-    struct MGIS_EXPORT IntegrateWorkSpace {
-      /*!
-       * \brief constructor
-       * \param[in] b: behaviour
-       */
-      IntegrateWorkSpace(const Behaviour&);
-      IntegrateWorkSpace(IntegrateWorkSpace&&);
-      IntegrateWorkSpace(const IntegrateWorkSpace&);
-      IntegrateWorkSpace& operator=(IntegrateWorkSpace&&);
-      IntegrateWorkSpace& operator=(const IntegrateWorkSpace&);
-      //! material properties at the beginning of the time step
-      std::vector<real> mps0;
-      //! material properties at the end of the time step
-      std::vector<real> mps1;
-      //! external state variables at the beginning of the time step
-      std::vector<real> esvs0;
-      //! external state variables at the end of the time step
-      std::vector<real> esvs1;
-    };  // end of struct IntegrateWorkSpace
+   * \brief structure defining various option
+   */
+  struct BehaviourIntegrationOptions {
+    //! \brief type of integration to be performed
+    IntegrationType integration_type =
+        IntegrationType::INTEGRATION_CONSISTENT_TANGENT_OPERATOR;
+    //! \brief if true, the speed of sound shall be computed
+    bool compute_speed_of_sound = false;
+  };  // end of BehaviourIntegrationOptions
 
+  /*!
+   * \brief structure in charge of reporting the result of a behaviour
+   * integration.
+   */
+  struct MGIS_EXPORT BehaviourIntegrationResult {
+    //! \brief default constructor
+    BehaviourIntegrationResult();
+    //! \brief move constructor
+    BehaviourIntegrationResult(BehaviourIntegrationResult&&);
+    //! \brief copye constructor
+    BehaviourIntegrationResult(const BehaviourIntegrationResult&);
+    //! \brief move assignement
+    BehaviourIntegrationResult& operator=(BehaviourIntegrationResult&&);
+    //! \brief copy assignement
+    BehaviourIntegrationResult& operator=(const BehaviourIntegrationResult&);
+    //! \brief destructor
+    ~BehaviourIntegrationResult();
+    /*! \brief exit status
+     *
+     * The returned value has the following meaning:
+     * - -1: integration failed for at least one Gauss point
+     * -  0: all integrations succeeded but results are unreliable for at least
+     *       one Gauss point
+     * -  1: integration succeeded and results are reliable.
+     */
+    int exit_status = 1;
+    //! \brief proposed time step increase factor
+    mgis::real time_step_increase_factor = 10;
     /*!
-     * \brief return a thread-specific workspace associated with the given
-     * behaviour.
-     * \param[in] b: behaviour
+     * \brief number of the integration point that failed or number of
+     * integration that reported unreliable results.
      */
-    MGIS_EXPORT IntegrateWorkSpace& getIntegrateWorkSpace(const Behaviour&);
+    mgis::size_type n = std::numeric_limits<mgis::size_type>::max();
+    //! \brief error message, if any
+    std::string error_message;
+  };  // end of struct BehaviourIntegrationResult
 
-    /*!
-     * \brief integrate the behaviour. The returned value has the following
-     * meaning:
-     * - -1: integration failed
-     * -  0: integration succeeded but results are unreliable
-     * -  1: integration succeeded and results are reliable
-     *
-     * \param[in,out] d: behaviour data
-     * \param[in,out] b: behaviour
+  /*!
+   * \brief structure in charge of reporting the result of a behaviour
+   * integration.
+   */
+  struct MGIS_EXPORT MultiThreadedBehaviourIntegrationResult {
+    //! \brief default constructor
+    MultiThreadedBehaviourIntegrationResult();
+    //! \brief move constructor
+    MultiThreadedBehaviourIntegrationResult(
+        MultiThreadedBehaviourIntegrationResult&&);
+    //! \brief copye constructor
+    MultiThreadedBehaviourIntegrationResult(
+        const MultiThreadedBehaviourIntegrationResult&);
+    //! \brief move assignement
+    MultiThreadedBehaviourIntegrationResult& operator=(
+        MultiThreadedBehaviourIntegrationResult&&);
+    //! \brief copy assignement
+    MultiThreadedBehaviourIntegrationResult& operator=(
+        const MultiThreadedBehaviourIntegrationResult&);
+    //! \brief destructor
+    ~MultiThreadedBehaviourIntegrationResult();
+    /*! \brief exit status
      *
-
-     * \note: the type of integration to be performed, must be
-     * explicitely set in d.K[0], as follows (see the IntegrationType finite ):
-     * - d.K[0]<-2.5, one computes a prediction and request the tangent operator
-     * - -2.5<d.K[0]<-1.5: one computes a prediction and request the secant
-     operator
-     * - -1.5<d.K[0]<-0.5: one computes a prediction and request the elastic
-     operator
-     * - -0.5<d.K[0]< 0.5: one integrates the behaviour over the time step
-     *                     but does not compute an stiffness tensor
-     * -  0.5<d.K[0]< 1.5: one integrates the behaviour over the time step
-     *                     and computes an elastic stiffness
-     * -  1.5<d.K[0]< 2.5: one integrates the behaviour over the time step
-     *                     and computes a secant stiffness
-     * -  2.5<d.K[0]< 3.5: one integrates the behaviour over the time step
-     *                     and computes a tangent stiffness
-     * -  2.5<d.K[0]< 3.5: one integrates the behaviour over the time step
-     *                     and computes a consistent tangent stiffness
+     * The returned value has the following meaning:
+     * - -1: integration failed for at least one Gauss point
+     * -  0: all integrations succeeded but results are unreliable for at least
+     *       one Gauss point
+     * -  1: integration succeeded and results are reliable.
      */
+    int exit_status = 1;
+    //! \brief integration results per threads
+    std::vector<BehaviourIntegrationResult> results;
+  };  // end of struct MultiThreadedBehaviourIntegrationResult
+
+  /*!
+   * \brief integrate the behaviour. The returned value has the following
+   * meaning:
+   * - -1: integration failed
+   * -  0: integration succeeded but results are unreliable
+   * -  1: integration succeeded and results are reliable
+   *
+   * \param[in,out] d: behaviour data
+   * \param[in,out] b: behaviour
+   *
+
+   * \note: the type of integration to be performed, must be
+   * explicitely set in d.K[0], as follows (see the `IntegrationType` enum).
+   *
+   * Let Ke be equal to:
+   *
+   * - d.K[0] - 100 if d.K[0] is greater than 50
+   * - d.K[0] otherwise.
+   *
+   * If Ke is negative, only the prediction operator is computed and
+   * no behaviour integration is performed.
+   *
+   * Ke has the following meaning:
+   *
+   * - if Ke is lower than -2.5, the tangent operator must be
+   *   computed.
+   * - if Ke is in [-2.5:-1.5]: the secant operator must be
+   *   computed.
+   * - if Ke is in [-1.5:-0.5]: the elastic operator must be
+   *   computed.
+   * - if Ke is in [-0.5:0.5]: the behaviour integration is
+   *   performed, but no stiffness matrix.
+   * - if Ke is in [0.5:1.5]: the elastic operator must be
+   *   computed.
+   * - if Ke is in [1.5:2.5]: the secant operator must be
+   *   computed.
+   * - if Ke is in [2.5:3.5]: the secant operator must be
+   *   computed.
+   * - if Ke is greater than 3.5, the consistent tangent operator
+   *   must be computed.
+   */
     int integrate(BehaviourDataView&, const Behaviour&);
 
     /*!
-     * \brief integrate the behaviour for a range of integration points. The
-     * returned value has the following meaning:
+     * \brief integrate the behaviour for a range of integration points.
+   * \return the result of the behaviour integration.
+   * \param[in,out] m: material data manager
+   * \param[in] c: description of the operation to be performed
+   * \param[in] dt: time step
+   * \param[in] b: first index of the range
+   * \param[in] e: last index of the range
+   *
+   * \note if required, the memory associated with the tangent operator blocks
+   * is automatically allocated.
+   */
+  MGIS_EXPORT BehaviourIntegrationResult
+  integrate(MaterialDataManager&,
+            const BehaviourIntegrationOptions&,
+            const real,
+            const size_type,
+            const size_type);
+
+  /*!
+   * \brief integrate the behaviour over all integration points using a thread
+   * pool to parallelize the integration.
+   * \return the result of the behaviour integration.
+   * \param[in,out] p: thread pool
+   * \param[in,out] m: material data manager
+   * \param[in] c: description of the operation to be performed
+   * \param[in] dt: time step
+   *
+   * \note if required, the memory associated with the tangent operator blocks
+   * is automatically allocated.
+   */
+  MGIS_EXPORT MultiThreadedBehaviourIntegrationResult
+  integrate(mgis::ThreadPool&,
+            MaterialDataManager&,
+            const BehaviourIntegrationOptions&,
+            const real);
+
+  /*!
+   * \brief integrate the behaviour for a range of integration points.
+     * \return an exit status. The returned value has the following meaning:
      * - -1: integration failed for at least one Gauss point
      * -  0: all integrations succeeded but results are unreliable for at least
      *       one Gauss point
@@ -117,6 +220,9 @@ namespace mgis {
      * \param[in] dt: time step
      * \param[in] b: first index of the range
      * \param[in] e: last index of the range
+     *
+     * \note if required, the memory associated with the tangent operator blocks
+     * is automatically allocated.
      */
     MGIS_EXPORT int integrate(MaterialDataManager&,
                               const IntegrationType,
@@ -125,9 +231,9 @@ namespace mgis {
                               const size_type);
 
     /*!
-     * \brief integrate the behaviour for a range of integration points. The
-     * returned value has the following meaning:
-     * - -1: integration failed for at least one Gauss point
+     * \brief integrate the behaviour for a range of integration points.
+     * \return an exit status. The returned value has the following meaning:
+     * - -1: integration failed for at least one integration point
      * -  0: all integrations succeeded but results are unreliable for at least
      *       one Gauss point
      * -  1: integration succeeded and results are reliable.
@@ -135,6 +241,9 @@ namespace mgis {
      * \param[in,out] p: thread pool
      * \param[in,out] m: material data manager
      * \param[in] dt: time step
+     *
+     * \note if required, the memory associated with the tangent operator blocks
+     * is automatically allocated.
      */
     MGIS_EXPORT int integrate(mgis::ThreadPool&,
                               MaterialDataManager&,
diff --git a/include/MGIS/Behaviour/MaterialDataManager.hxx b/include/MGIS/Behaviour/MaterialDataManager.hxx
index 8dfa8574f..b6cd74595 100644
--- a/include/MGIS/Behaviour/MaterialDataManager.hxx
+++ b/include/MGIS/Behaviour/MaterialDataManager.hxx
@@ -15,6 +15,9 @@
 #ifndef LIB_MGIS_BEHAVIOUR_MATERIALDATAMANAGER_HXX
 #define LIB_MGIS_BEHAVIOUR_MATERIALDATAMANAGER_HXX
 
+#include <map>
+#include <thread>
+#include <memory>
 #include <vector>
 #include "MGIS/Config.hxx"
 #include "MGIS/Variant.hxx"
@@ -27,6 +30,35 @@ namespace mgis {
     // forward declaration
     struct Behaviour;
 
+    //! \brief structure in charge of handling temporary memory access.
+    struct MGIS_EXPORT BehaviourIntegrationWorkSpace {
+      /*!
+       * \brief constructor
+       * \param[in] b: behaviour
+       */
+      BehaviourIntegrationWorkSpace(const Behaviour&);
+      //! \brief move constructor
+      BehaviourIntegrationWorkSpace(BehaviourIntegrationWorkSpace&&);
+      //! \brief copye constructor
+      BehaviourIntegrationWorkSpace(const BehaviourIntegrationWorkSpace&);
+      //! \brief move assignement
+      BehaviourIntegrationWorkSpace& operator=(BehaviourIntegrationWorkSpace&&);
+      //! \brief copy assignement
+      BehaviourIntegrationWorkSpace& operator=(
+          const BehaviourIntegrationWorkSpace&);
+    //! \brief destructor
+    ~BehaviourIntegrationWorkSpace();
+      //! \brief a buffer to hold error messages
+      std::vector<char> error_message;
+      //! material properties at the beginning of the time step
+      std::vector<mgis::real> mps0;
+      //! material properties at the end of the time step
+      std::vector<mgis::real> mps1;
+      //! external state variables at the beginning of the time step
+      std::vector<mgis::real> esvs0;
+      //! external state variables at the end of the time step
+      std::vector<mgis::real> esvs1;
+    };  // end of struct BehaviourIntegrationWorkSpace
     /*!
      * \brief a structure in charge of holding information on how a material
      * data manager shall be initialized.
@@ -39,9 +71,15 @@ namespace mgis {
       /*!
        * \brief view to an externally allocated memory used to store the
        * tangent operator. If empty, the material data manager will
-       * initialize the required memory internally.
+       * initialize the required memory internally if required.
        */
       mgis::span<mgis::real> K;
+      /*!
+       * \brief view to an externally allocated memory used to store the
+       * speed_of_sound. If empty, the material data manager will
+       * initialize the required memory internally if required.
+       */
+      mgis::span<mgis::real> speed_of_sound;
       /*!
        * \brief object used to initalize the state manager associated with the
        * beginning of the time step.
@@ -80,14 +118,102 @@ namespace mgis {
       MaterialDataManager(const Behaviour&,
                           const size_type,
                           const MaterialDataManagerInitializer&);
-      //! destructor
+      /*!
+       * \brief set if the `MaterialDataManager` must take care of
+       * thread-safety. This flag is mostly used in members functions allocating
+       * memory. \param[in] bv: boolean
+       */
+      void setThreadSafe(const bool);
+      /*!
+       * \brief allocate the memory associated with the tangent operator blocks
+       * if required.
+       *
+       * This method is useless if the memory associated with the tangent
+       * operator blocks had previously been allocated or assigned to external
+       * memory (see the `MaterialDataManagerInitializer` structure).
+       *
+       * \note This method is thread-safe if the `thread_safe` is `true`.
+       * In this case, the memory allocation is guarded by a mutex.
+       * See the `setThreadSafe` method for details
+       */
+      void allocateArrayOfTangentOperatorBlocks();
+      /*!
+       * \brief use an externally allocated memory to store the tangent operator
+       * blocks.
+       *
+       * \param[in] m: memory view
+       *
+       * \note this method calls `releaseArrayOfTangentOperatorBlocks` before
+       * allocating the memory.
+       */
+      void useExternalArrayOfTangentOperatorBlocks(mgis::span<real>);
+      /*!
+       * \brief release the memory associated with the tangent operator blocks.
+       *
+       * If the memory associated with the tangent operator blocks was handled
+       * internally, this memory is freed.
+       * If an external memory buffer was used, reference to this buffer is
+       * removed.
+       */
+      void releaseArrayOfTangentOperatorBlocks();
+      /*!
+       * \brief allocate the memory associated with the speed of sound if
+       * required.
+       *
+       * This method is useless if the memory associated with the speed of sound
+       * had previously been allocated or assigned to external memory (see the
+       * `MaterialDataManagerInitializer` structure).
+       *
+       * \note This method is thread-safe if the `thread_safe` is `true`.
+       * In this case, the memory allocation is guarded by a mutex.
+       * See the `setThreadSafe` method for details
+       */
+      void allocateArrayOfSpeedOfSounds();
+      /*!
+       * \brief use an externally allocated memory to store the tangent operator
+       * blocks.
+       *
+       * \param[in] m: memory view
+       *
+       * \note this method calls `releaseArrayOfSpeedOfSounds` before
+       * allocating the memory.
+       */
+      void useExternalArrayOfSpeedOfSounds(mgis::span<real>);
+      /*!
+       * \brief release the memory associated with the tangent operator blocks.
+       *
+       * If the memory associated with the tangent operator blocks was handled
+       * internally, this memory is freed.
+       * If an external memory buffer was used, reference to this buffer is
+       * removed.
+       */
+      void releaseArrayOfSpeedOfSounds();
+      /*!
+       * \brief return a workspace associated with the given behaviour.
+       *
+       * \note This method returns a object per thread if the `thread_safe` is
+       * `true`.
+       */
+      BehaviourIntegrationWorkSpace& getBehaviourIntegrationWorkSpace();
+      /*!
+       * \brief clear behaviour integration workspaces.
+       *
+       * \note This is only interesting if you keep allocating different thread
+       * pools.
+       */
+      void releaseBehaviourIntegrationWorkspaces();
+      //! \brief destructor
       ~MaterialDataManager();
       //! \brief state at the beginning of the time step
       MaterialStateManager s0;
       //! \brief state at the end of the time step
       MaterialStateManager s1;
-      //! \brief view of the stiffness matrices.
+      //! \brief view of the stiffness matrices, if required.
       mgis::span<real> K;
+      //! \brief proposed time step increment increase factor
+      real rdt;
+      //! \brief view on the speed of sound.
+      mgis::span<real> speed_of_sound;
       //! \brief number of integration points
       const size_type n;
       /*!
@@ -95,12 +221,10 @@ namespace mgis {
        * size of K is K_stride times the number of integration points)
        */
       const size_type K_stride;
-      //! underlying behaviour
+      //! \brief underlying behaviour
       const Behaviour& b;
 
      private:
-      //! \brief values of the stiffness matrices, if hold internally.
-      std::vector<real> K_values;
       //! move constructor
       MaterialDataManager(MaterialDataManager&&) = delete;
       //! copy constructor
@@ -109,6 +233,17 @@ namespace mgis {
       MaterialDataManager& operator=(MaterialDataManager&&) = delete;
       //! copy assignement
       MaterialDataManager& operator=(const MaterialDataManager&) = delete;
+      //! \brief values of the stiffness matrices, if hold internally.
+      std::vector<real> K_values;
+      //! \brief values of the speed of sound, if hold internally.
+      std::vector<real> speed_of_sound_values;
+      //! \brief integration workspace for individual threads.
+      std::map<std::thread::id, std::unique_ptr<BehaviourIntegrationWorkSpace>>
+          iwks;
+      //! \brief a pointer to an integration workspace
+      std::unique_ptr<BehaviourIntegrationWorkSpace> iwk;
+      //! \brief boolean stating if thread safety must be unsured
+      bool thread_safe = true;
     };  // end of struct MaterialDataManager
 
     /*!
diff --git a/include/MGIS/Config-c.h b/include/MGIS/Config-c.h
index 4a8df41c2..334a72a3a 100644
--- a/include/MGIS/Config-c.h
+++ b/include/MGIS/Config-c.h
@@ -21,6 +21,21 @@
 #include <stddef.h>
 #endif /*  __cplusplus */
 
+/*!
+ * \brief define the API used by the behaviour support.
+ *
+ * The API mostly refers to the various functions provided by `MGIS` but also
+ * the data strutures used to manage behaviours (See `BehaviourData`,
+ * `MaterialDataManger`, etc...).
+ *
+ * - 0 corresponds to the initial api provided by the TFEL project
+ *   (versions 3.2.x, versions 3.3.x and 3.4.0)`.
+ * - 1 corresponds the initial api provided by the TFEL project (versions 3.4.1
+ *   and after that and 4.x).
+ *
+ */
+#define MGIS_BEHAVIOUR_API_VERSION 1
+
 /*!
  * Macro extracted from :
  * "Why is the new C++ visibility support so useful?"
@@ -53,10 +68,10 @@
 #endif /* defined _WIN32 || defined _WIN64 ||defined __CYGWIN__ */
 
 #ifdef MGIS_REAL_TYPE
-/*! alias to the numeric type used in the library */
+/*! \brief alias to the numeric type used in the library */
 typedef MGIS_REAL_TYPE mgis_real;
-#else   /* MGIS_REAL_TYPE */
-/*! alias to the numeric type used in the library */
+#else  /* MGIS_REAL_TYPE */
+/*! \brief alias to the numeric type used in the library */
 typedef double mgis_real;
 #endif /* MGIS_REAL_TYPE */
 
diff --git a/include/MGIS/LibrariesManager.hxx b/include/MGIS/LibrariesManager.hxx
index f18679c33..d193f8755 100644
--- a/include/MGIS/LibrariesManager.hxx
+++ b/include/MGIS/LibrariesManager.hxx
@@ -59,6 +59,15 @@ namespace mgis {
      * \param[in] n: entry point name
      */
     std::string getTFELVersion(const std::string &, const std::string &);
+    /*!
+     * \return the `TFEL`'s API version.
+     *
+     * This was added in TFEL 3.4.1 which broke the generic behaviour' API.
+     *
+     * \note If the corresponding symbol is not present, the API is assumed to
+     * conform to the TFEL' API prior to Version 3.4.1 and 0 is returned.
+     */
+    unsigned short getAPIVersion(const std::string &, const std::string &);
     /*!
      * \return the material knowledge type of an entry point in a
      * library. The returned value has the following meaning:
diff --git a/include/MGIS/ThreadedTaskResult.hxx b/include/MGIS/ThreadedTaskResult.hxx
index 3448f705e..9eee5ea69 100644
--- a/include/MGIS/ThreadedTaskResult.hxx
+++ b/include/MGIS/ThreadedTaskResult.hxx
@@ -18,6 +18,7 @@
 #include <exception>
 #include <type_traits>
 #include "MGIS/Config.hxx"
+#include "MGIS/Variant.hxx"
 
 namespace mgis {
 
@@ -85,22 +86,8 @@ namespace mgis {
     MGIS_INLINE ~ThreadedTaskResult();
 
    private:
-    /*!
-     * build an object
-     * \param[in] args: arguments to the constructor
-     */
-    template <typename... Args>
-    MGIS_INLINE void build(Args&&...);
-    //! \brief clear the underlying object destructor
-    MGIS_INLINE void clear();
-    //! \brief pointer the underlying object storage
-    MGIS_INLINE T* get_pointer();
-    //! \brief pointer the underlying object storage
-    MGIS_INLINE const T* get_pointer() const;
-    //! \brief if true, the ThreadedTaskResult contains a value
-    bool initialized = false;
-    //! \brief storage of the underlying value
-    typename std::aligned_storage<sizeof(T), alignof(T)>::type storage;
+    //! \brief result of the task
+    mgis::variant<mgis::monostate, T> result;
     //! exception ptr thrown during the task
     std::exception_ptr eptr;
   };
diff --git a/include/MGIS/ThreadedTaskResult.ixx b/include/MGIS/ThreadedTaskResult.ixx
index 1114ba087..5abdbd902 100644
--- a/include/MGIS/ThreadedTaskResult.ixx
+++ b/include/MGIS/ThreadedTaskResult.ixx
@@ -26,86 +26,41 @@ namespace mgis {
   template <typename T>
   template <typename... Args>
   ThreadedTaskResult<T>::ThreadedTaskResult(Args&&... args) {
-    this->build(args...);
+    this->result.template emplace<T>(std::forward<Args...>(args...));
   }  // end of ThreadedTaskResult<T>::ThreadedTaskResult
 
   template <typename T>
-  ThreadedTaskResult<T>::ThreadedTaskResult(ThreadedTaskResult&& src) {
-    this->eptr = std::move(src.eptr);
-    if (src.initialized) {
-      this->initialized = true;
-      this->build(std::move(*(src.get_pointer())));
-    }
-  }
+  ThreadedTaskResult<T>::ThreadedTaskResult(ThreadedTaskResult&&) = default;
 
   template <typename T>
-  ThreadedTaskResult<T>::ThreadedTaskResult(const ThreadedTaskResult& src) {
-    this->eptr = src.eptr;
-    if (src.initialized) {
-      this->initialized = true;
-      this->build(src.get_pointer());
-    }
-  }
+  ThreadedTaskResult<T>::ThreadedTaskResult(const ThreadedTaskResult&) =
+      default;
 
   template <typename T>
   ThreadedTaskResult<T>& ThreadedTaskResult<T>::operator=(
-      ThreadedTaskResult&& src) {
-    if (this != &src) {
-      this->eptr = std::move(src.eptr);
-      if (src.initialized) {
-        if (this->initialized) {
-          *(this->get_pointer()) = std::move(*(src.get_pointer()));
-        } else {
-          this->build(std::move(*(src.get_pointer())));
-        }
-      } else {
-        this->clear();
-      }
-    }
-    return *this;
-  }
+      ThreadedTaskResult&&) = default;
 
   template <typename T>
   ThreadedTaskResult<T>& ThreadedTaskResult<T>::operator=(
-      const ThreadedTaskResult& src) {
-    if (this != &src) {
-      this->eptr = src.eptr;
-      if (src.initialized) {
-        if (this->initialized) {
-          *(this->get_pointer()) = *(src.get_pointer());
-        } else {
-          this->build(*(src.get_pointer()));
-        }
-      } else {
-        this->clear();
-      }
-    }
-    return *this;
-  }
+      const ThreadedTaskResult&) = default;
 
   template <typename T>
   ThreadedTaskResult<T>& ThreadedTaskResult<T>::operator=(const T& src) {
-    if (this->initialized) {
-      *(this->get_pointer()) = src;
-    } else {
-      this->build(src);
-    }
+    this->result = src;
     return *this;
   }
 
   template <typename T>
   ThreadedTaskResult<T>& ThreadedTaskResult<T>::operator=(T&& src) {
-    if (this->initialized) {
-      *(this->get_pointer()) = std::forward<T>(src);
-    } else {
-      this->build(std::forward<T>(src));
-    }
+    this->result = std::forward<T>(src);
     return *this;
   }
 
   template <typename T>
   ThreadedTaskResult<T>::operator bool() const {
-    return this->initialized && (this->eptr == nullptr);
+    const auto initialized =
+        mgis::holds_alternative<mgis::monostate>(this->result);
+    return initialized && (this->eptr == nullptr);
   }
 
   template <typename T>
@@ -113,10 +68,12 @@ namespace mgis {
     if (this->eptr != nullptr) {
       this->rethrow();
     }
-    if (!this->initialized) {
+    const auto initialized =
+        !mgis::holds_alternative<mgis::monostate>(this->result);
+    if (!initialized) {
       ThreadedTaskResultBase::throwBadCastException();
     }
-    return *(this->get_pointer());
+    return mgis::get<T>(this->result);
   }
 
   template <typename T>
@@ -124,37 +81,27 @@ namespace mgis {
     if (this->eptr != nullptr) {
       this->rethrow();
     }
-    if (!this->initialized) {
+    const auto initialized =
+        !mgis::holds_alternative<mgis::monostate>(this->result);
+    if (initialized) {
       ThreadedTaskResultBase::throwBadCastException();
     }
-    return *(this->get_pointer());
+    return mgis::get<T>(this->result);
   }
 
   template <typename T>
   T* ThreadedTaskResult<T>::operator->() {
-    if (this->eptr != nullptr) {
-      this->rethrow();
-    }
-    if (!this->initialized) {
-      ThreadedTaskResultBase::throwBadCastException();
-    }
-    return this->get_pointer();
-  }
+    return &(*(*this));
+  } // end of operator->()
 
   template <typename T>
   const T* ThreadedTaskResult<T>::operator->() const {
-    if (this->eptr != nullptr) {
-      this->rethrow();
-    }
-    if (!this->initialized) {
-      ThreadedTaskResultBase::throwBadCastException();
-    }
-    return this->get_pointer();
-  }
+    return &(*(*this));
+  } // end of operator->()
 
   template <typename T>
   void ThreadedTaskResult<T>::setException(const std::exception_ptr& e) {
-    this->clear();
+    this->result = mgis::monostate{};
     this->eptr = e;
   }  // end of ThreadedTaskResult<T>::setException
 
@@ -167,37 +114,7 @@ namespace mgis {
   }  // end of ThreadedTaskResult<T>::setException
 
   template <typename T>
-  template <typename... Args>
-  void ThreadedTaskResult<T>::build(Args&&... args) {
-    this->clear();
-    new (this->get_pointer()) T(args...);
-    this->initialized = true;
-  }
-
-  template <typename T>
-  void ThreadedTaskResult<T>::clear() {
-    this->eptr = nullptr;
-    if (this->initialized) {
-      auto ptr = this->get_pointer();
-      ptr->~T();
-      this->initialized = false;
-    }
-  }  // end of clear
-
-  template <typename T>
-  T* ThreadedTaskResult<T>::get_pointer() {
-    return reinterpret_cast<T*>(::std::addressof(this->storage));
-  }  // end of get_pointer
-
-  template <typename T>
-  const T* ThreadedTaskResult<T>::get_pointer() const {
-    return reinterpret_cast<const T*>(::std::addressof(this->storage));
-  }  // end of get_pointer
-
-  template <typename T>
-  ThreadedTaskResult<T>::~ThreadedTaskResult() {
-    this->clear();
-  }
+  ThreadedTaskResult<T>::~ThreadedTaskResult() = default;
 
 }  // end of namespace mgis
 
diff --git a/src/Behaviour.cxx b/src/Behaviour.cxx
index f4ad94742..22d49310e 100644
--- a/src/Behaviour.cxx
+++ b/src/Behaviour.cxx
@@ -169,6 +169,15 @@ namespace mgis {
               " is not a behaviour");
       }
 
+      if (lm.getAPIVersion(l, b) != MGIS_BEHAVIOUR_API_VERSION) {
+        std::string msg("unmatched API version\n");
+        msg += "- the behaviour uses API version ";
+        msg += std::to_string(lm.getAPIVersion(l, b)) + "\n";
+        msg += "- mgis uses API version ";
+        msg += std::to_string(MGIS_BEHAVIOUR_API_VERSION);
+        raise(msg);
+      }
+
       d.tfel_version = lm.getTFELVersion(l, b);
       d.source = lm.getSource(l, b);
       d.btype = [&l, &b, &lm, &raise] {
diff --git a/src/BehaviourData.cxx b/src/BehaviourData.cxx
index 9eba53d63..0690e2930 100644
--- a/src/BehaviourData.cxx
+++ b/src/BehaviourData.cxx
@@ -37,7 +37,9 @@ namespace mgis {
     // declared in the class definition (again regardless of the order of the
     // mem-initializers).
     BehaviourData::BehaviourData(const Behaviour& b)
-        : dt(0), rdt(1), s0(b), s1(s0) {
+        : dt(0), rdt(1), speed_of_sound(0), s0(b), s1(s0) {
+      static char error_message_buffer[512];
+      this->error_message = error_message_buffer;
       this->K.resize(getTangentOperatorArraySize(b));
     }  // end of Behaviour::Behaviour
 
@@ -61,8 +63,10 @@ namespace mgis {
         return &v[0];
       };  // end of get_ptr
       BehaviourDataView v;
+      v.error_message = d.error_message;
       v.dt = d.dt;
-      v.rdt = d.rdt;
+      v.rdt = &(d.rdt);
+      v.speed_of_sound = &(d.speed_of_sound);
       v.K = get_ptr(d.K);
       v.s0 = make_view(static_cast<const State&>(d.s0));
       v.s1 = make_view(d.s1);
diff --git a/src/Integrate.cxx b/src/Integrate.cxx
index cd1fe05f9..4312b9f39 100644
--- a/src/Integrate.cxx
+++ b/src/Integrate.cxx
@@ -26,41 +26,72 @@ namespace mgis {
 
   namespace behaviour {
 
-    IntegrateWorkSpace::IntegrateWorkSpace(const Behaviour& b)
-        : mps0(getArraySize(b.mps, b.hypothesis)),
-          mps1(getArraySize(b.mps, b.hypothesis)),
-          esvs0(getArraySize(b.esvs, b.hypothesis)),
-          esvs1(getArraySize(b.esvs, b.hypothesis)) {
-    }  // end of IntegrateWorkSpace
-
-    IntegrateWorkSpace::IntegrateWorkSpace(IntegrateWorkSpace&&) = default;
-    IntegrateWorkSpace::IntegrateWorkSpace(const IntegrateWorkSpace&) = default;
-    IntegrateWorkSpace& IntegrateWorkSpace::operator=(IntegrateWorkSpace&&) =
-        default;
-    IntegrateWorkSpace& IntegrateWorkSpace::operator=(
-        const IntegrateWorkSpace&) = default;
-
-    IntegrateWorkSpace& getIntegrateWorkSpace(const Behaviour& b) {
-      static std::map<
-          const Behaviour*,
-          std::map<std::thread::id, std::shared_ptr<IntegrateWorkSpace>>>
-          m;
-      static std::mutex mt;
-      const auto id = std::this_thread::get_id();
-      std::lock_guard<std::mutex> lock(mt);
-      auto& mwks = m[&b];
-      auto p = mwks.find(id);
-      if (p == mwks.end()) {
-        p = mwks.insert({id, std::make_shared<IntegrateWorkSpace>(b)}).first;
+    static void allocate(MaterialDataManager& m,
+                         const BehaviourIntegrationOptions& opts) {
+      if (opts.integration_type !=
+          IntegrationType::INTEGRATION_NO_TANGENT_OPERATOR) {
+        m.allocateArrayOfTangentOperatorBlocks();
       }
-      return *(p->second);
-    }  // end of getIntegrateWorkSpace
+      if (opts.compute_speed_of_sound) {
+        m.allocateArrayOfSpeedOfSounds();
+      }
+    }  // end of allocate
 
-    int integrate(MaterialDataManager& m,
-                  const IntegrationType it,
-                  const real dt,
-                  const size_type b,
-                  const size_type e) {
+    static mgis::real encodeBehaviourIntegrationOptions(
+        const BehaviourIntegrationOptions& opts) {
+      if (opts.compute_speed_of_sound) {
+        return 100 + static_cast<int>(opts.integration_type);
+      }
+      return static_cast<int>(opts.integration_type);
+    }  // end of encodeBehaviourIntegrationOptions
+
+    BehaviourIntegrationResult::BehaviourIntegrationResult() = default;
+
+    BehaviourIntegrationResult::BehaviourIntegrationResult(
+        BehaviourIntegrationResult&&) = default;
+
+    BehaviourIntegrationResult::BehaviourIntegrationResult(
+        const BehaviourIntegrationResult&) = default;
+
+    BehaviourIntegrationResult& BehaviourIntegrationResult::operator=(
+        BehaviourIntegrationResult&&) = default;
+
+    BehaviourIntegrationResult& BehaviourIntegrationResult::operator=(
+        const BehaviourIntegrationResult&) = default;
+
+    BehaviourIntegrationResult::~BehaviourIntegrationResult() = default;
+
+    MultiThreadedBehaviourIntegrationResult::
+        MultiThreadedBehaviourIntegrationResult() = default;
+
+    MultiThreadedBehaviourIntegrationResult::
+        MultiThreadedBehaviourIntegrationResult(
+            MultiThreadedBehaviourIntegrationResult&&) = default;
+
+    MultiThreadedBehaviourIntegrationResult::
+        MultiThreadedBehaviourIntegrationResult(
+            const MultiThreadedBehaviourIntegrationResult&) = default;
+
+    MultiThreadedBehaviourIntegrationResult&
+    MultiThreadedBehaviourIntegrationResult::operator=(
+        MultiThreadedBehaviourIntegrationResult&&) = default;
+
+    MultiThreadedBehaviourIntegrationResult&
+    MultiThreadedBehaviourIntegrationResult::operator=(
+        const MultiThreadedBehaviourIntegrationResult&) = default;
+
+    MultiThreadedBehaviourIntegrationResult::
+        ~MultiThreadedBehaviourIntegrationResult() = default;
+
+    /*!
+     * \brief integration the behaviour over a range of integration points.
+     */
+    static BehaviourIntegrationResult integrate2(
+        MaterialDataManager& m,
+        const BehaviourIntegrationOptions& opts,
+        const real dt,
+        const size_type b,
+        const size_type e) {
       /*
        * \brief uniform values are treated immediatly. For spatially variable
        * fields, we return the information needed to evaluate them
@@ -115,12 +146,13 @@ namespace mgis {
           v[std::get<0>(ev)] = std::get<1>(ev)[i];
         }
       };  // end of eval
+
       // strides
       const auto g_stride = m.s0.gradients_stride;
       const auto t_stride = m.s0.thermodynamic_forces_stride;
       const auto isvs_stride = m.s0.internal_state_variables_stride;
       // workspace
-      auto& ws = getIntegrateWorkSpace(m.b);
+      auto& ws = m.getBehaviourIntegrationWorkSpace();
       // treating uniform values
       const auto vmps0 = dispatch(ws.mps0, m.s0.material_properties, m.b.mps);
       const auto vmps1 = dispatch(ws.mps1, m.s1.material_properties, m.b.mps);
@@ -131,17 +163,24 @@ namespace mgis {
       const auto computes_stored_energy = m.b.computesStoredEnergy;
       const auto computes_dissipated_energy = m.b.computesDissipatedEnergy;
       // loop over integration points
-      auto r = int{1};
-      real opts[Behaviour::nopts + 1];  // option passed to the behaviour
+      auto r = BehaviourIntegrationResult{};
+      auto rdt0 = r.time_step_increase_factor;
+      const real Ke = encodeBehaviourIntegrationOptions(opts);
+      real bopts[Behaviour::nopts + 1];  // option passed to the behaviour
       for (auto i = b; i != e; ++i) {
+        auto rdt = rdt0;
         BehaviourDataView v;
-        v.rdt = real(1);
+        v.error_message = ws.error_message.data();
+        v.error_message[0] = '\0';
+        v.rdt = &rdt;
         v.dt = dt;
-        if (it != IntegrationType::INTEGRATION_NO_TANGENT_OPERATOR) {
+        if (opts.integration_type !=
+            IntegrationType::INTEGRATION_NO_TANGENT_OPERATOR) {
           v.K = m.K.data() + m.K_stride * i;
         } else {
-          v.K = &opts[0];
+          v.K = &bopts[0];
         }
+        v.speed_of_sound = m.speed_of_sound.data() + i;
         eval(ws.mps0, vmps0, i);
         eval(ws.mps1, vmps1, i);
         eval(ws.esvs0, vesvs0, i);
@@ -174,55 +213,84 @@ namespace mgis {
         }
         v.s0.external_state_variables = ws.esvs0.data();
         v.s1.external_state_variables = ws.esvs1.data();
-        v.K[0] = static_cast<int>(it);
-        switch (integrate(v, m.b)) {
-          case 1:
-            r = std::min(r, 1);
-            break;
-          case 0:
-            r = 0;
-            break;
-          default:
-            return -1;
+        v.K[0] = Ke;
+        const auto ri = integrate(v, m.b);
+        r.exit_status = std::min(ri, r.exit_status);
+        r.time_step_increase_factor =
+            std::min(rdt, r.time_step_increase_factor);
+        if (r.exit_status == -1) {
+          r.n = i;
+          v.error_message[511] = '\0';
+          r.error_message = std::string(v.error_message);
+          return r;
         }
       }
       return r;
+    }  // end of integrate2
+
+    int integrate(MaterialDataManager& m,
+                  const IntegrationType it,
+                  const real dt,
+                  const size_type b,
+                  const size_type e) {
+      BehaviourIntegrationOptions opts;
+      opts.integration_type = it;
+      const auto r = integrate(m, opts, dt, b, e);
+      return r.exit_status;
+    }  // end of integrate
+
+    BehaviourIntegrationResult integrate(
+        MaterialDataManager& m,
+        const BehaviourIntegrationOptions& opts,
+        const real dt,
+        const size_type b,
+        const size_type e) {
+      allocate(m, opts);
+      return integrate2(m, opts, dt, b, e);
     }  // end of integrate
 
     int integrate(ThreadPool& p,
                   MaterialDataManager& m,
                   const IntegrationType it,
                   const real dt) {
+      BehaviourIntegrationOptions opts;
+      opts.integration_type = it;
+      const auto r = integrate(p, m, opts, dt);
+      return r.exit_status;
+    }  // end of integrate
+
+    MultiThreadedBehaviourIntegrationResult integrate(
+        ThreadPool& p,
+        MaterialDataManager& m,
+        const BehaviourIntegrationOptions& opts,
+        const real dt) {
+      m.setThreadSafe(true);
+      allocate(m, opts);
       // get number of threads
       const auto nth = p.getNumberOfThreads();
       const auto d = m.n / nth;
       const auto r = m.n % nth;
       size_type b = 0;
-      std::vector<std::future<ThreadedTaskResult<int>>> tasks;
+      std::vector<std::future<ThreadedTaskResult<BehaviourIntegrationResult>>>
+          tasks;
       tasks.reserve(nth);
       for (size_type i = 0; i != r; ++i) {
-        tasks.push_back(p.addTask(
-            [&m, it, dt, b, d] { return integrate(m, it, dt, b, b + d + 1); }));
+        tasks.push_back(p.addTask([&m, &opts, dt, b, d] {
+          return integrate2(m, opts, dt, b, b + d + 1);
+        }));
         b += d + 1;
       }
       for (size_type i = r; i != nth; ++i) {
-        tasks.push_back(p.addTask(
-            [&m, it, dt, b, d] { return integrate(m, it, dt, b, b + d); }));
+        tasks.push_back(p.addTask([&m, &opts, dt, b, d] {
+          return integrate2(m, opts, dt, b, b + d);
+        }));
         b += d;
       }
-      //      p.wait();
-      auto res = int{1};
+      auto res = MultiThreadedBehaviourIntegrationResult{};
       for (auto& t : tasks) {
-        switch (*(t.get())) {
-          case 1:
-            res = std::min(res, 1);
-            break;
-          case 0:
-            res = 0;
-            break;
-          default:
-            return -1;
-        }
+        const auto& ri = *(t.get());
+        res.exit_status = std::min(res.exit_status, ri.exit_status);
+        res.results.push_back(ri);
       }
       return res;
     }  // end of integrate
diff --git a/src/LibrariesManager.cxx b/src/LibrariesManager.cxx
index 0be82faa7..2df44df06 100644
--- a/src/LibrariesManager.cxx
+++ b/src/LibrariesManager.cxx
@@ -58,7 +58,7 @@ namespace mgis {
     }
     return "";
 #endif /* (defined _WIN32 || defined _WIN64) && (!defined __CYGWIN__) */
-  }    // end of  getErrorMessage
+  }    // end of getErrorMessage
 
   static LibrariesManager::libhandler load_library(const std::string &l) {
 #if (defined _WIN32 || defined _WIN64) && (!defined __CYGWIN__)
@@ -159,7 +159,7 @@ namespace mgis {
   LibrariesManager &LibrariesManager::get() {
     static LibrariesManager lm;
     return lm;
-  }  // end of LibrariesManager::get
+  }  // end of get
 
   LibrariesManager::LibrariesManager() = default;
 
@@ -174,7 +174,7 @@ namespace mgis {
           b + "' in library '" + l + "' for hypothesis '" + hn + "'");
     }
     return reinterpret_cast<mgis::behaviour::BehaviourFctPtr>(p);
-  }  // end of LibrariesManager::getBehaviour
+  }  // end of getBehaviour
 
   mgis::behaviour::RotateBehaviourGradientsFctPtr
   LibrariesManager::getRotateBehaviourGradientsFunction(const std::string &l,
@@ -193,7 +193,7 @@ namespace mgis {
           b + "' in library '" + l + "' for hypothesis '" + hn + "'");
     }
     return reinterpret_cast<mgis::behaviour::RotateBehaviourGradientsFctPtr>(p);
-  }  // end of LibrariesManager::getRotateBehaviourGradientsFunction
+  }  // end of getRotateBehaviourGradientsFunction
 
   mgis::behaviour::RotateArrayOfBehaviourGradientsFctPtr
   LibrariesManager::getRotateArrayOfBehaviourGradientsFunction(
@@ -213,7 +213,7 @@ namespace mgis {
     }
     return reinterpret_cast<
         mgis::behaviour::RotateArrayOfBehaviourGradientsFctPtr>(p);
-  }  // end of LibrariesManager::getRotateArrayOfBehaviourGradientsFunction
+  }  // end of getRotateArrayOfBehaviourGradientsFunction
 
   mgis::behaviour::RotateBehaviourThermodynamicForcesFctPtr
   LibrariesManager::getRotateBehaviourThermodynamicForcesFunction(
@@ -230,7 +230,7 @@ namespace mgis {
     }
     return reinterpret_cast<
         mgis::behaviour::RotateBehaviourThermodynamicForcesFctPtr>(p);
-  }  // end of LibrariesManager::getRotateBehaviourThermodynamicForcesFunction
+  }  // end of getRotateBehaviourThermodynamicForcesFunction
 
   mgis::behaviour::RotateArrayOfBehaviourThermodynamicForcesFctPtr
   LibrariesManager::getRotateArrayOfBehaviourThermodynamicForcesFunction(
@@ -248,8 +248,7 @@ namespace mgis {
     }
     return reinterpret_cast<
         mgis::behaviour::RotateArrayOfBehaviourThermodynamicForcesFctPtr>(p);
-  }  // end of
-     // LibrariesManager::getRotateArrayOfBehaviourThermodynamicForcesFunction
+  }  // end of getRotateArrayOfBehaviourThermodynamicForcesFunction
 
   mgis::behaviour::RotateBehaviourThermodynamicForcesFctPtr
   LibrariesManager::getRotateBehaviourThermodynamicForcesFunction(
@@ -281,7 +280,7 @@ namespace mgis {
     }
     return reinterpret_cast<
         mgis::behaviour::RotateBehaviourThermodynamicForcesFctPtr>(p);
-  }  // end of LibrariesManager::getRotateBehaviourThermodynamicForcesFunction
+  }  // end of getRotateBehaviourThermodynamicForcesFunction
 
   mgis::behaviour::RotateArrayOfBehaviourThermodynamicForcesFctPtr
   LibrariesManager::getRotateArrayOfBehaviourThermodynamicForcesFunction(
@@ -314,8 +313,7 @@ namespace mgis {
     }
     return reinterpret_cast<
         mgis::behaviour::RotateArrayOfBehaviourThermodynamicForcesFctPtr>(p);
-  }  // end of
-     // LibrariesManager::getRotateArrayOfBehaviourThermodynamicForcesFunction
+  }  // end of getRotateArrayOfBehaviourThermodynamicForcesFunction
 
   mgis::behaviour::RotateBehaviourTangentOperatorBlocksFctPtr
   LibrariesManager::getRotateBehaviourTangentOperatorBlocksFunction(
@@ -332,7 +330,7 @@ namespace mgis {
     }
     return reinterpret_cast<
         mgis::behaviour::RotateBehaviourTangentOperatorBlocksFctPtr>(p);
-  }  // end of LibrariesManager::getRotateBehaviourTangentOperatorBlocksFunction
+  }  // end of getRotateBehaviourTangentOperatorBlocksFunction
 
   mgis::behaviour::RotateArrayOfBehaviourTangentOperatorBlocksFctPtr
   LibrariesManager::getRotateArrayOfBehaviourTangentOperatorBlocksFunction(
@@ -350,8 +348,7 @@ namespace mgis {
     }
     return reinterpret_cast<
         mgis::behaviour::RotateArrayOfBehaviourTangentOperatorBlocksFctPtr>(p);
-  }  // end of
-     // LibrariesManager::getRotateArrayOfBehaviourTangentOperatorBlocksFunction
+  }  // end of getRotateArrayOfBehaviourTangentOperatorBlocksFunction
 
   mgis::behaviour::RotateBehaviourTangentOperatorBlocksFctPtr
   LibrariesManager::getRotateBehaviourTangentOperatorBlocksFunction(
@@ -383,7 +380,7 @@ namespace mgis {
     }
     return reinterpret_cast<
         mgis::behaviour::RotateBehaviourTangentOperatorBlocksFctPtr>(p);
-  }  // end of LibrariesManager::getRotateBehaviourTangentOperatorBlocksFunction
+  }  // end of getRotateBehaviourTangentOperatorBlocksFunction
 
   mgis::behaviour::RotateArrayOfBehaviourTangentOperatorBlocksFctPtr
   LibrariesManager::getRotateArrayOfBehaviourTangentOperatorBlocksFunction(
@@ -418,8 +415,7 @@ namespace mgis {
     }
     return reinterpret_cast<
         mgis::behaviour::RotateArrayOfBehaviourTangentOperatorBlocksFctPtr>(p);
-  }  // end of
-     // LibrariesManager::getRotateArrayOfBehaviourTangentOperatorBlocksFunction
+  }  // end of getRotateArrayOfBehaviourTangentOperatorBlocksFunction
 
   std::string LibrariesManager::getTFELVersion(const std::string &l,
                                                const std::string &n) {
@@ -428,12 +424,22 @@ namespace mgis {
       return "";
     }
     return *(static_cast<const char *const *>(p));
-  }  // end of LibrariesManager::getTFELVersion
+  }  // end of getTFELVersion
+
+  unsigned short LibrariesManager::getAPIVersion(const std::string &l,
+                                                 const std::string &n) {
+    const auto s = n + "_api_version";
+    const auto p = this->getSymbolAddress(l, s);
+    if (p == nullptr) {
+      return 0;
+    }
+    return *(static_cast<unsigned short *>(p));
+  }  // end of getAPIVersion
 
   unsigned short LibrariesManager::getMaterialKnowledgeType(
       const std::string &l, const std::string &b) {
     return *(this->extract<unsigned short>(l, b + "_mfront_mkt"));
-  }  // end of LibrariesManager::getMaterialKnowledgeType
+  }  // end of getMaterialKnowledgeType
 
   std::string LibrariesManager::getSource(const std::string &l,
                                           const std::string &n) {
@@ -442,7 +448,7 @@ namespace mgis {
       return "";
     }
     return *(static_cast<const char *const *>(p));
-  }  // end of LibrariesManager::getSource
+  }  // end of getSource
 
   std::string LibrariesManager::getInterface(const std::string &l,
                                              const std::string &n) {
@@ -451,27 +457,27 @@ namespace mgis {
       return "";
     }
     return *(static_cast<const char *const *>(p));
-  }  // end of LibrariesManager::getInterface
+  }  // end of getInterface
 
   unsigned short LibrariesManager::getBehaviourType(const std::string &l,
                                                     const std::string &b) {
     return *(this->extract<unsigned short>(l, b + "_BehaviourType"));
-  }  // end of LibrariesManager::getBehaviourType
+  }  // end of getBehaviourType
 
   unsigned short LibrariesManager::getBehaviourKinematic(const std::string &l,
                                                          const std::string &b) {
     return *(this->extract<unsigned short>(l, b + "_BehaviourKinematic"));
-  }  // end of LibrariesManager::getBehaviourKinematic
+  }  // end of getBehaviourKinematic
 
   unsigned short LibrariesManager::getBehaviourSymmetry(const std::string &l,
                                                         const std::string &b) {
     return *(this->extract<unsigned short>(l, b + "_SymmetryType"));
-  }  // end of LibrariesManager::getBehaviourSymmetry
+  }  // end of getBehaviourSymmetry
 
   unsigned short LibrariesManager::getElasticStiffnessSymmetry(
       const std::string &l, const std::string &b) {
     return *(this->extract<unsigned short>(l, b + "_ElasticSymmetryType"));
-  }  // end of LibrariesManager::getElasticStiffnessSymmetry
+  }  // end of getElasticStiffnessSymmetry
 
   bool LibrariesManager::requiresStiffnessTensor(const std::string &l,
                                                  const std::string &b,
@@ -480,7 +486,7 @@ namespace mgis {
     const auto bv =
         *(this->extract<unsigned short>(l, b + "_" + toString(h) + sn, b + sn));
     return bv == 1 ? true : false;
-  }  // end of LibrariesManager::requiresStiffnessTensor
+  }  // end of requiresStiffnessTensor
 
   bool LibrariesManager::computesStoredEnergy(const std::string &l,
                                               const std::string &b,
@@ -489,7 +495,7 @@ namespace mgis {
     const auto bv =
         *(this->extract<unsigned short>(l, b + "_" + toString(h) + sn, b + sn));
     return bv == 1 ? true : false;
-  }  // end of LibrariesManager::computesStoredEnergy
+  }  // end of computesStoredEnergy
 
   bool LibrariesManager::computesDissipatedEnergy(const std::string &l,
                                                   const std::string &b,
@@ -498,7 +504,7 @@ namespace mgis {
     const auto bv =
         *(this->extract<unsigned short>(l, b + "_" + toString(h) + sn, b + sn));
     return bv == 1 ? true : false;
-  }  // end of LibrariesManager::computesDissipatedEnergy
+  }  // end of computesDissipatedEnergy
 
   bool LibrariesManager::requiresThermalExpansionCoefficientTensor(
       const std::string &l, const std::string &b, const Hypothesis h) {
@@ -506,7 +512,7 @@ namespace mgis {
     const auto bv =
         *(this->extract<unsigned short>(l, b + "_" + toString(h) + sn, b + sn));
     return bv == 1 ? true : false;
-  }  // end of LibrariesManager::requiresThermalExpansionCoefficientTensor
+  }  // end of requiresThermalExpansionCoefficientTensor
 
   template <typename T>
   const T *LibrariesManager::extract(const std::string &l,
@@ -516,7 +522,7 @@ namespace mgis {
       raise("LibrariesManager::extract: could not load symbol '" + n + "'");
     }
     return static_cast<const T *const>(p);
-  }  // end of LibrariesManager::extract
+  }  // end of extract
 
   template <typename T>
   const T *LibrariesManager::extract(const std::string &l,
@@ -528,7 +534,7 @@ namespace mgis {
             "' nor '" + n2 + "'");
     }
     return static_cast<const T *const>(p);
-  }  // end of LibrariesManager::extract
+  }  // end of extract
 
   std::vector<std::string> LibrariesManager::getNames(const std::string &l,
                                                       const std::string &f,
@@ -542,12 +548,12 @@ namespace mgis {
                                                       f + '_' + n);
     std::copy(res, res + nb, std::back_inserter(vars));
     return vars;
-  }  // end of LibrariesManager::getNames
+  }  // end of getNames
 
   std::vector<std::string> LibrariesManager::getGradientsNames(
       const std::string &l, const std::string &b, const Hypothesis h) {
     return this->getNames(l, b, h, "Gradients");
-  }  // end of LibrariesManager::getGradientsNames
+  }  // end of getGradientsNames
 
   std::vector<int> LibrariesManager::getGradientsTypes(const std::string &l,
                                                        const std::string &b,
@@ -560,12 +566,12 @@ namespace mgis {
         l, b + "_" + hn + "_GradientsTypes", b + "_GradientsTypes");
     std::copy(res, res + nb, std::back_inserter(types));
     return types;
-  }  // end of LibrariesManager::getGradientsTypes
+  }  // end of getGradientsTypes
 
   std::vector<std::string> LibrariesManager::getThermodynamicForcesNames(
       const std::string &l, const std::string &b, const Hypothesis h) {
     return this->getNames(l, b, h, "ThermodynamicForces");
-  }  // end of LibrariesManager::getThermodynamicForcesNames
+  }  // end of getThermodynamicForcesNames
 
   std::vector<int> LibrariesManager::getThermodynamicForcesTypes(
       const std::string &l, const std::string &b, const Hypothesis h) {
@@ -579,7 +585,7 @@ namespace mgis {
                                  b + "_ThermodynamicForcesTypes");
     std::copy(res, res + nb, std::back_inserter(types));
     return types;
-  }  // end of LibrariesManager::getThermodynamicForcesTypes
+  }  // end of getThermodynamicForcesTypes
 
   std::vector<std::pair<std::string, std::string>>
   LibrariesManager::getTangentOperatorBlocksNames(const std::string &l,
@@ -596,17 +602,17 @@ namespace mgis {
       blocks.push_back({names[2 * i], names[2 * i + 1]});
     }
     return blocks;
-  }  // end of LibrariesManager::getTangentOperatorBlocksNames
+  }  // end of getTangentOperatorBlocksNames
 
   std::vector<std::string> LibrariesManager::getMaterialPropertiesNames(
       const std::string &l, const std::string &b, const Hypothesis h) {
     return this->getNames(l, b, h, "MaterialProperties");
-  }  // end of LibrariesManager::getMaterialPropertiesNames
+  }  // end of getMaterialPropertiesNames
 
   std::vector<std::string> LibrariesManager::getInternalStateVariablesNames(
       const std::string &l, const std::string &b, const Hypothesis h) {
     return this->getNames(l, b, h, "InternalStateVariables");
-  }  // end of LibrariesManager::getInternalStateVariablesNames
+  }  // end of getInternalStateVariablesNames
 
   std::vector<int> LibrariesManager::getInternalStateVariablesTypes(
       const std::string &l, const std::string &b, const Hypothesis h) {
@@ -620,16 +626,16 @@ namespace mgis {
         b + "_InternalStateVariablesTypes");
     std::copy(res, res + nb, std::back_inserter(types));
     return types;
-  }  // end of LibrariesManager::getInternalStateVariablesTypes
+  }  // end of getInternalStateVariablesTypes
 
   std::vector<std::string> LibrariesManager::getExternalStateVariablesNames(
       const std::string &l, const std::string &b, const Hypothesis h) {
     return this->getNames(l, b, h, "ExternalStateVariables");
-  }  // end of LibrariesManager::getMaterialPropertiesNames
+  }  // end of getMaterialPropertiesNames
 
   bool LibrariesManager::contains(const std::string &l, const std::string &n) {
     return this->getSymbolAddress(l, n) != nullptr;
-  }  // end of LibrariesManager::contains
+  }  // end of contains
 
   void *LibrariesManager::getSymbolAddress(const std::string &l,
                                            const std::string &n) {
@@ -671,7 +677,7 @@ namespace mgis {
       return lib;
     }
     return p->second;
-  }  // end of LibrariesManager::loadLibrary
+  }  // end of loadLibrary
 
   void LibrariesManager::setParameter(const std::string &l,
                                       const std::string &b,
@@ -695,7 +701,7 @@ namespace mgis {
           "call to the '" +
           b + "_setParameter' function failed");
     }
-  }  // end of LibrariesManager::setParameter
+  }  // end of setParameter
 
   void LibrariesManager::setParameter(const std::string &l,
                                       const std::string &b,
@@ -719,7 +725,7 @@ namespace mgis {
           "call to the '" +
           b + "_setParameter' function failed");
     }
-  }  // end of LibrariesManager::setParameter
+  }  // end of setParameter
 
   void LibrariesManager::setParameter(const std::string &l,
                                       const std::string &b,
@@ -744,12 +750,12 @@ namespace mgis {
           "call to the '" +
           b + "_setParameter' function failed");
     }
-  }  // end of LibrariesManager::setParameter
+  }  // end of setParameter
 
   std::vector<std::string> LibrariesManager::getParametersNames(
       const std::string &l, const std::string &b, const Hypothesis h) {
     return this->getNames(l, b, h, "Parameters");
-  }  // end of LibrariesManager::getMaterialPropertiesNames
+  }  // end of getMaterialPropertiesNames
 
   std::vector<int> LibrariesManager::getParametersTypes(const std::string &l,
                                                         const std::string &b,
@@ -762,7 +768,7 @@ namespace mgis {
         l, b + "_" + hn + "_ParametersTypes", b + "_ParametersTypes");
     std::copy(res, res + nb, std::back_inserter(types));
     return types;
-  }  // end of LibrariesManager::getInternalVariablesTypes
+  }  // end of getInternalVariablesTypes
 
   double LibrariesManager::getParameterDefaultValue(const std::string &l,
                                                     const std::string &b,
@@ -773,7 +779,7 @@ namespace mgis {
     return *(this->extract<double>(
         l, b + "_" + hn + "_" + pn + "_ParameterDefaultValue",
         b + "_" + pn + "_ParameterDefaultValue"));
-  }  // end of LibrariesManager::getParameterDefaultValue
+  }  // end of getParameterDefaultValue
 
   int LibrariesManager::getIntegerParameterDefaultValue(const std::string &l,
                                                         const std::string &b,
@@ -784,7 +790,7 @@ namespace mgis {
     return *(this->extract<int>(
         l, b + "_" + hn + "_" + pn + "_ParameterDefaultValue",
         b + "_" + pn + "_ParameterDefaultValue"));
-  }  // end of LibrariesManager::getIntegerParameterDefaultValue
+  }  // end of getIntegerParameterDefaultValue
 
   unsigned short LibrariesManager::getUnsignedShortParameterDefaultValue(
       const std::string &l,
@@ -796,7 +802,7 @@ namespace mgis {
     return *(this->extract<unsigned short>(
         l, b + "_" + hn + "_" + pn + "_ParameterDefaultValue",
         b + "_" + pn + "_ParameterDefaultValue"));
-  }  // end of LibrariesManager::getUnsignedShortParameterDefaultValue
+  }  // end of getUnsignedShortParameterDefaultValue
 
   bool LibrariesManager::hasBounds(const std::string &l,
                                    const std::string &b,
@@ -810,7 +816,7 @@ namespace mgis {
     const auto n4 = b + "_" + vn + "_UpperBound";
     return ((this->contains(l, n1)) || (this->contains(l, n2)) ||
             (this->contains(l, n3)) || (this->contains(l, n4)));
-  }  // end of LibrariesManager::hasBounds
+  }  // end of hasBounds
 
   bool LibrariesManager::hasLowerBound(const std::string &l,
                                        const std::string &b,
@@ -821,7 +827,7 @@ namespace mgis {
     const auto n1 = b + "_" + hn + "_" + vn + "_LowerBound";
     const auto n2 = b + "_" + vn + "_LowerBound";
     return ((this->contains(l, n1)) || (this->contains(l, n2)));
-  }  // end of LibrariesManager::hasLowerBound
+  }  // end of hasLowerBound
 
   bool LibrariesManager::hasUpperBound(const std::string &l,
                                        const std::string &b,
@@ -832,7 +838,7 @@ namespace mgis {
     const auto n1 = b + "_" + hn + "_" + vn + "_UpperBound";
     const auto n2 = b + "_" + vn + "_UpperBound";
     return ((this->contains(l, n1)) || (this->contains(l, n2)));
-  }  // end of LibrariesManager::hasUpperBound
+  }  // end of hasUpperBound
 
   long double LibrariesManager::getLowerBound(const std::string &l,
                                               const std::string &b,
@@ -843,7 +849,7 @@ namespace mgis {
     const auto n1 = b + "_" + hn + "_" + vn + "_LowerBound";
     const auto n2 = b + "_" + vn + "_LowerBound";
     return *(this->extract<long double>(l, n1, n2));
-  }  // end of LibrariesManager::getLowerBound
+  }  // end of getLowerBound
 
   long double LibrariesManager::getUpperBound(const std::string &l,
                                               const std::string &b,
@@ -854,7 +860,7 @@ namespace mgis {
     const auto n1 = b + "_" + hn + "_" + vn + "_UpperBound";
     const auto n2 = b + "_" + vn + "_UpperBound";
     return *(this->extract<long double>(l, n1, n2));
-  }  // end of LibrariesManager::getUpperBound
+  }  // end of getUpperBound
 
   bool LibrariesManager::hasPhysicalBounds(const std::string &l,
                                            const std::string &b,
@@ -868,7 +874,7 @@ namespace mgis {
     const auto n4 = b + "_" + vn + "_UpperPhysicalBound";
     return ((this->contains(l, n1)) || (this->contains(l, n2)) ||
             (this->contains(l, n3)) || (this->contains(l, n4)));
-  }  // end of LibrariesManager::hasPhysicalBounds
+  }  // end of hasPhysicalBounds
 
   bool LibrariesManager::hasLowerPhysicalBound(const std::string &l,
                                                const std::string &b,
@@ -879,7 +885,7 @@ namespace mgis {
     const auto n1 = b + "_" + hn + "_" + vn + "_LowerPhysicalBound";
     const auto n2 = b + "_" + vn + "_LowerPhysicalBound";
     return ((this->contains(l, n1)) || (this->contains(l, n2)));
-  }  // end of LibrariesManager::hasLowerPhysicalBound
+  }  // end of hasLowerPhysicalBound
 
   bool LibrariesManager::hasUpperPhysicalBound(const std::string &l,
                                                const std::string &b,
@@ -890,7 +896,7 @@ namespace mgis {
     const auto n1 = b + "_" + hn + "_" + vn + "_UpperPhysicalBound";
     const auto n2 = b + "_" + vn + "_UpperPhysicalBound";
     return ((this->contains(l, n1)) || (this->contains(l, n2)));
-  }  // end of LibrariesManager::hasUpperPhysicalBound
+  }  // end of hasUpperPhysicalBound
 
   long double LibrariesManager::getLowerPhysicalBound(const std::string &l,
                                                       const std::string &b,
@@ -901,7 +907,7 @@ namespace mgis {
     const auto n1 = b + "_" + hn + "_" + vn + "_LowerPhysicalBound";
     const auto n2 = b + "_" + vn + "_LowerPhysicalBound";
     return *(this->extract<long double>(l, n1, n2));
-  }  // end of LibrariesManager::getLowerPhysicalBound
+  }  // end of getLowerPhysicalBound
 
   long double LibrariesManager::getUpperPhysicalBound(const std::string &l,
                                                       const std::string &b,
@@ -912,7 +918,7 @@ namespace mgis {
     const auto n1 = b + "_" + hn + "_" + vn + "_UpperPhysicalBound";
     const auto n2 = b + "_" + vn + "_UpperPhysicalBound";
     return *(this->extract<long double>(l, n1, n2));
-  }  // end of LibrariesManager::getUpperPhysicalBound
+  }  // end of getUpperPhysicalBound
 
   LibrariesManager::~LibrariesManager() {
     for (const auto &l : this->libraries) {
@@ -922,6 +928,6 @@ namespace mgis {
       ::dlclose(l.second);
 #endif /* (defined _WIN32 || defined _WIN64) && (!defined __CYGWIN__) */
     }
-  }  // end of LibrariesManager::~LibrariesManager
+  }  // end of ~LibrariesManager
 
 }  // namespace mgis
diff --git a/src/MaterialDataManager.cxx b/src/MaterialDataManager.cxx
index 5008487af..afa7e1dc8 100644
--- a/src/MaterialDataManager.cxx
+++ b/src/MaterialDataManager.cxx
@@ -1,5 +1,5 @@
 /*!
- * \file   MaterialDataManager.cxx
+ * \file   src/MaterialDataManager.cxx
  * \brief
  * \author Thomas Helfer
  * \date   21/08/2018
@@ -12,6 +12,8 @@
  *   CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt).
  */
 
+#include <mutex>
+#include <thread>
 #include "MGIS/Raise.hxx"
 #include "MGIS/Behaviour/Behaviour.hxx"
 #include "MGIS/Behaviour/MaterialDataManager.hxx"
@@ -20,17 +22,32 @@ namespace mgis {
 
   namespace behaviour {
 
+    BehaviourIntegrationWorkSpace::BehaviourIntegrationWorkSpace(
+        const Behaviour& b)
+        : error_message(512),
+          mps0(getArraySize(b.mps, b.hypothesis)),
+          mps1(getArraySize(b.mps, b.hypothesis)),
+          esvs0(getArraySize(b.esvs, b.hypothesis)),
+          esvs1(getArraySize(b.esvs, b.hypothesis)) {
+    }  // end of BehaviourIntegrationWorkSpace
+
+    BehaviourIntegrationWorkSpace::BehaviourIntegrationWorkSpace(
+        BehaviourIntegrationWorkSpace&&) = default;
+    BehaviourIntegrationWorkSpace::BehaviourIntegrationWorkSpace(
+        const BehaviourIntegrationWorkSpace&) = default;
+    BehaviourIntegrationWorkSpace& BehaviourIntegrationWorkSpace::operator=(
+        BehaviourIntegrationWorkSpace&&) = default;
+    BehaviourIntegrationWorkSpace& BehaviourIntegrationWorkSpace::operator=(
+        const BehaviourIntegrationWorkSpace&) = default;
+  BehaviourIntegrationWorkSpace::~BehaviourIntegrationWorkSpace() = default;
+
     MaterialDataManager::MaterialDataManager(const Behaviour& behaviour,
                                              const size_type s)
         : s0(behaviour, s),
           s1(behaviour, s),
           n(s),
           K_stride(getTangentOperatorArraySize(behaviour)),
-          b(behaviour) {
-      constexpr const auto zero = real{0};
-      this->K_values.resize(this->n * this->K_stride, zero);
-      this->K = mgis::span<real>(this->K_values);
-    }  // end of Behaviour::Behaviour
+          b(behaviour) {}  // end of MaterialDataManager
 
     MaterialDataManager::MaterialDataManager(
         const Behaviour& behaviour,
@@ -41,21 +58,113 @@ namespace mgis {
           n(s),
           K_stride(getTangentOperatorArraySize(behaviour)),
           b(behaviour) {
-      if (i.K.empty()) {
+      if (!i.K.empty()) {
+        this->useExternalArrayOfTangentOperatorBlocks(i.K);
+      }
+      if (!i.speed_of_sound.empty()) {
+        this->useExternalArrayOfSpeedOfSounds(i.speed_of_sound);
+      }
+    }  // end of MaterialDataManager
+
+    static void allocateArrayWithoutSynchronization(
+        mgis::span<real>& v,
+        std::vector<mgis::real>& values,
+        const mgis::size_type s) {
+      if (v.empty()) {
         constexpr const auto zero = real{0};
-        this->K_values.resize(this->n * this->K_stride, zero);
-        this->K = mgis::span<real>(this->K_values);
+        values.resize(s, zero);
+        v = mgis::span<real>(values);
+      }
+    }  // end of allocateArrayWithoutSynchronization
+
+    static void allocateArrayWithSynchronization(
+        mgis::span<real>& v,
+        std::vector<mgis::real>& values,
+        const mgis::size_type s) {
+      static std::mutex mt;
+      std::lock_guard<std::mutex> lock(mt);
+      allocateArrayWithoutSynchronization(v, values, s);
+    }  // end of allocateArrayWithSynchronization
+
+    void MaterialDataManager::setThreadSafe(const bool bv) {
+      this->thread_safe = bv;
+    }  // end of setThreadSafe
+    void MaterialDataManager::allocateArrayOfTangentOperatorBlocks() {
+      if (this->thread_safe) {
+        allocateArrayWithSynchronization(this->K, this->K_values,
+                                         this->n * this->K_stride);
+      } else {
+        allocateArrayWithoutSynchronization(this->K, this->K_values,
+                                            this->n * this->K_stride);
+      }
+    }  // end of allocateArrayOfTangentOperatorBlocks
+
+    void MaterialDataManager::releaseArrayOfTangentOperatorBlocks() {
+      this->K = mgis::span<real>();
+      this->K_values.clear();
+    }  // end of releaseArrayOfTangentOperatorBlocks
+
+    void MaterialDataManager::useExternalArrayOfTangentOperatorBlocks(
+        mgis::span<real> m) {
+      if (m.size() != this->n * this->K_stride) {
+        mgis::raise(
+            "MaterialDataManager::useExternalArrayOfTangentOperatorBlocks: "
+            "the external memory has not been allocated properly");
+      }
+      this->releaseArrayOfTangentOperatorBlocks();
+      this->K = m;
+    }  // end of useExternalArrayOfTangentOperatorBlocks
+    void MaterialDataManager::allocateArrayOfSpeedOfSounds() {
+      if (this->thread_safe) {
+        allocateArrayWithSynchronization(this->speed_of_sound,
+                                         this->speed_of_sound_values, this->n);
       } else {
-        if (i.K.size() != this->n * this->K_stride) {
-          mgis::raise(
-              "MaterialDataManager::MaterialDataManager: "
-              "the consistent tangent operator has not been allocated "
-              "properly");
+        allocateArrayWithoutSynchronization(
+            this->speed_of_sound, this->speed_of_sound_values, this->n);
+      }
+    }  // end of allocateArrayOfSpeedOfSounds
+
+    void MaterialDataManager::releaseArrayOfSpeedOfSounds() {
+      this->speed_of_sound = mgis::span<real>();
+      this->speed_of_sound_values.clear();
+    }  // end of releaseArrayOfSpeedOfSounds
+
+    void MaterialDataManager::useExternalArrayOfSpeedOfSounds(
+        mgis::span<real> m) {
+      if (m.size() != this->n) {
+        mgis::raise(
+            "MaterialDataManager::useExternalArrayOfSpeedOfSounds: "
+            "the external memory has not been allocated properly");
+      }
+      this->releaseArrayOfSpeedOfSounds();
+      this->speed_of_sound = m;
+    }  // end of useExternalArrayOfSpeedOfSounds
+
+    BehaviourIntegrationWorkSpace&
+    MaterialDataManager::getBehaviourIntegrationWorkSpace() {
+      if (this->thread_safe) {
+        static std::mutex mt;
+        std::lock_guard<std::mutex> lock(mt);
+        const auto id = std::this_thread::get_id();
+        auto p = this->iwks.find(id);
+        if (p == this->iwks.end()) {
+          std::unique_ptr<BehaviourIntegrationWorkSpace> wk(
+              new BehaviourIntegrationWorkSpace(b));
+          p = this->iwks.insert({id, std::move(wk)}).first;
         }
-        this->K = i.K;
+        return *(p->second);
+      }
+      if (this->iwk == nullptr) {
+        this->iwk = std::unique_ptr<BehaviourIntegrationWorkSpace>(
+            new BehaviourIntegrationWorkSpace(b));
       }
-    }  // end of Behaviour::Behaviour
+      return *(this->iwk);
+    }  // end of getBehaviourIntegrationWorkSpace
 
+    void MaterialDataManager::releaseBehaviourIntegrationWorkspaces() {
+      this->iwk.reset();
+      this->iwks.clear();
+    }  // end of releaseBehaviourIntegrationWorkspaces
     MaterialDataManager::~MaterialDataManager() = default;
 
     void update(MaterialDataManager& m) {
diff --git a/src/ThreadedTaskResult.cxx b/src/ThreadedTaskResult.cxx
index d8034daf8..ecf8d53a2 100644
--- a/src/ThreadedTaskResult.cxx
+++ b/src/ThreadedTaskResult.cxx
@@ -20,7 +20,7 @@
 namespace mgis {
 
   void ThreadedTaskResultBase::throwBadCastException() {
-    throw(std::bad_cast());
+    raise<std::bad_cast>();
   }  // end of ThreadPool::ResultBase::throwBadCastException
 
   void ThreadedTaskResultBase::throwNullException() {
diff --git a/tests/IntegrateTest.cxx b/tests/IntegrateTest.cxx
index f00bf797a..68171ef1b 100644
--- a/tests/IntegrateTest.cxx
+++ b/tests/IntegrateTest.cxx
@@ -44,8 +44,8 @@ int main(const int argc, const char* const* argv) {
     auto p = std::array<real, 21>{};
     p[0] = d.s0.internal_state_variables[o];
     for (size_type i = 0; i != 20; ++i) {
+      d.rdt = 1;
       auto v = make_view(d);
-      v.rdt = 1;
       integrate(v, b);
       update(d);
       d.s1.gradients[0] += de;