diff --git a/.gitmodules b/.gitmodules
index 36e7c5e246cd3144c82581e6220736ae84698819..84dceb42bea9b8e8076166c3d2e10c0ddd0c119c 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -37,3 +37,6 @@
 [submodule "ThirdParty/iphreeqc/src"]
 	path = ThirdParty/iphreeqc/src
 	url = https://github.com/ufz/iphreeqc.git
+[submodule "ThirdParty/MGIS"]
+	path = ThirdParty/MGIS
+	url = https://github.com/ufz/MFrontGenericInterfaceSupport.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index d998bdb9fc5be797b888a11ad75462657815fa79..5c2b0d86235e823cf9d7f692c1c2906e754ea4c4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -46,6 +46,7 @@ set(OGS_LIBS
     petsc
     lis
     cvode
+    tfel
     CACHE INTERNAL "")
 foreach(LIB ${OGS_LIBS})
     set(OGS_LIB_${LIB}
@@ -255,7 +256,9 @@ endif()
 option(OGS_USE_MFRONT
        "Enable solid material models by MFront (https://tfel.sourceforge.net)"
        OFF)
-
+if(OGS_USE_MFRONT)
+    add_definitions(-DOGS_USE_MFRONT)
+endif()
 # ---- Definitions ----
 if(OGS_USE_LIS)
     add_definitions(-DUSE_LIS)
diff --git a/Documentation/ProjectFile/material/solid/constitutive_relation/MFront/c_MFront.md b/Documentation/ProjectFile/material/solid/constitutive_relation/MFront/c_MFront.md
index d840e4e555e06432eacfe0eeca91c80f51771f39..24c2a9d9f93011f860c19fc65b65a2fab175170d 100644
--- a/Documentation/ProjectFile/material/solid/constitutive_relation/MFront/c_MFront.md
+++ b/Documentation/ProjectFile/material/solid/constitutive_relation/MFront/c_MFront.md
@@ -3,7 +3,11 @@ with MFront's generic interface.
 
 See http://tfel.sourceforge.net/ https://github.com/thelfer/MFrontGenericInterfaceSupport
 
-This requires MFront 3.2.0 or newer (because of MFront's "generic" interface).
+This requires MFront version 3.2.1 (which must correspond to the MFront's
+"generic" interface), which must be available on your system: see
+http://tfel.sourceforge.net page for download and installation intstructions.
+
+\note The MFront library is distributed under GPL license.
 
 \attention
 OpenGeoSys initializes the internal state of the MFront material model to
diff --git a/Jenkinsfile b/Jenkinsfile
index 8e9a9ebe6dc26455f65084fab780e46a87d3536e..c3feecd8b11d5112ea288c799b7ed29488c48bd6 100644
--- a/Jenkinsfile
+++ b/Jenkinsfile
@@ -101,6 +101,7 @@ pipeline {
           }
           environment {
             OMP_NUM_THREADS = '1'
+            LD_LIBRARY_PATH = "$WORKSPACE/build/lib"
           }
           steps {
             script {
@@ -111,7 +112,9 @@ pipeline {
                   '-DOGS_CPU_ARCHITECTURE=generic ' +
                   '-DOGS_BUILD_UTILS=ON ' +
                   '-DOGS_CONAN_BUILD=outdated ' +
-                  '-DOGS_USE_CVODE=ON '
+                  '-DOGS_USE_CVODE=ON ' +
+                  '-DOGS_USE_MFRONT=ON ' +
+                  '-DOGS_USE_PYTHON=ON '
               }
               build {
                 target="package"
@@ -139,6 +142,10 @@ pipeline {
               ])
               recordIssues enabledForFailure: true, filters: [
                 excludeFile('.*qrc_icons\\.cpp.*'), excludeFile('.*QVTKWidget.*'),
+                excludeFile('.*\\.conan.*/TFEL/.*'),
+                excludeFile('.*ThirdParty/MGIS/src.*'),
+                excludeFile('.*MaterialLib/SolidModels/MFront/.*\\.mfront'),
+                excludeFile('.*MaterialLib/SolidModels/MFront/.*\\.hxx'),
                 excludeMessage('.*tmpnam.*')],
                 tools: [gcc4(name: 'GCC', pattern: 'build/build*.log')],
                 unstableTotalAll: 1
diff --git a/MaterialLib/CMakeLists.txt b/MaterialLib/CMakeLists.txt
index 6eb3884711e3f306ee990f9a14ad5092a49c9acc..548d6d176f0737e7a5dcde115963843ba3447ec2 100644
--- a/MaterialLib/CMakeLists.txt
+++ b/MaterialLib/CMakeLists.txt
@@ -29,7 +29,9 @@ append_source_files(SOURCES
                     PorousMedium/UnsaturatedProperty/RelativePermeability)
 append_source_files(SOURCES TwoPhaseModels)
 
-add_subdirectory(SolidModels/MFront)
+if(OGS_USE_MFRONT)
+    add_subdirectory(SolidModels/MFront)
+endif()
 
 add_library(MaterialLib ${SOURCES})
 if(BUILD_SHARED_LIBS)
@@ -40,9 +42,11 @@ include(GenerateExportHeader)
 generate_export_header(MaterialLib)
 target_include_directories(MaterialLib PUBLIC ${CMAKE_CURRENT_BINARY_DIR})
 
-target_link_libraries(MaterialLib
-                      PUBLIC MaterialLib_SolidModels_MFront
-                      PRIVATE MathLib MeshLib ParameterLib)
+target_link_libraries(MaterialLib PRIVATE MathLib MeshLib ParameterLib)
+
+if(OGS_USE_MFRONT)
+    target_link_libraries(MaterialLib PUBLIC MaterialLib_SolidModels_MFront)
+endif()
 
 if(OGS_USE_PCH)
     cotire(MaterialLib)
diff --git a/MaterialLib/SolidModels/CreateConstitutiveRelation.cpp b/MaterialLib/SolidModels/CreateConstitutiveRelation.cpp
index 17a786b940597932b1600ceaaf614c959a50fe85..be94fe171f77e266a8be004252f68242409ffc9b 100644
--- a/MaterialLib/SolidModels/CreateConstitutiveRelation.cpp
+++ b/MaterialLib/SolidModels/CreateConstitutiveRelation.cpp
@@ -16,7 +16,9 @@
 #include "CreateLinearElasticIsotropic.h"
 #include "CreateLinearElasticOrthotropic.h"
 #include "CreateLubby2.h"
+#ifdef OGS_USE_MFRONT
 #include "MFront/CreateMFront.h"
+#endif  // OGS_USE_MFRONT
 
 #include "MechanicsBase.h"
 
@@ -71,8 +73,14 @@ createConstitutiveRelation(
     }
     if (type == "MFront")
     {
+#ifdef OGS_USE_MFRONT
         return MaterialLib::Solids::MFront::createMFront<DisplacementDim>(
             parameters, config);
+#else   // OGS_USE_MFRONT
+        OGS_FATAL(
+            "OGS is compiled without MFront support. See OGS_USE_MFRONT CMake "
+            "option.");
+#endif  // OGS_USE_MFRONT
     }
     OGS_FATAL("Cannot construct constitutive relation of given type '%s'.",
               type.c_str());
diff --git a/MaterialLib/SolidModels/MFront/BDT.mfront b/MaterialLib/SolidModels/MFront/BDT.mfront
new file mode 100644
index 0000000000000000000000000000000000000000..88bd6d008346ff9920eae568e2987929752e0861
--- /dev/null
+++ b/MaterialLib/SolidModels/MFront/BDT.mfront
@@ -0,0 +1,161 @@
+@Parser Implicit;
+@Behaviour BDT;
+@Algorithm NewtonRaphson_NumericalJacobian;
+@Theta 1.;
+
+@MaterialProperty real young;
+young.setGlossaryName("YoungModulus");
+@MaterialProperty real nu;
+nu.setGlossaryName("PoissonRatio");
+@MaterialProperty thermalexpansion alpha;
+alpha.setGlossaryName("ThermalExpansion");
+@MaterialProperty real f_c;
+@MaterialProperty real m_0;
+@MaterialProperty real alpha_p;
+@MaterialProperty real n_exp_T;
+@MaterialProperty real q_h0;
+@MaterialProperty real chi_h;
+@MaterialProperty real alpha_d;
+@MaterialProperty real h_d;
+@MaterialProperty real Qact;
+@MaterialProperty real A_creep;
+@MaterialProperty real n_creep;
+@MaterialProperty real El_1;
+@MaterialProperty real El_2;
+@MaterialProperty real El_3;
+@MaterialProperty real at_1;
+@MaterialProperty real at_2;
+@MaterialProperty real at_3;
+
+@StateVariable real pla_mult;
+@StateVariable real a;
+@StateVariable real dam;
+@StateVariable real Rv;
+@StateVariable real triax_p;
+@StateVariable real vp;
+
+@LocalVariable stress lambda;
+@LocalVariable stress mu;
+@LocalVariable stress Fel;
+@LocalVariable StressStensor s0;
+
+// @Parameter Tref;
+// Tref.setDefaultValue(293.15);
+
+@InitLocalVariables
+{
+    auto const T_Celsius = T - 273.15;
+
+    const auto y2 = El_1 * power<2>(T_Celsius) + El_2 * T_Celsius + El_3;
+    lambda = computeLambda(young * y2, nu);
+    mu = computeMu(young * y2, nu);
+    StressStensor sigel(lambda * trace(eel + deto) * Stensor::Id() +
+                        2 * mu * (eel + deto));  // StresssStensor in Tutorial
+    const auto s_dev = deviator(sigel);
+    // J2 = (s_dev|s_dev)/2.; //defines double contraction
+    const stress seq = sigmaeq(s_dev);
+    const stress I1 = trace(sigel);
+    const auto aux = (seq + I1) / (3. * f_c);
+
+    const auto aux_ex = (1. - 1. / n_exp_T);
+    const auto aux_arg = 1 + pow(alpha_p * (T_Celsius - 10.), n_exp_T);
+
+    Rv = max(at_1 * T_Celsius + at_2 * I1 / 3. + at_3, 1.0e-4);
+    auto Rpel = -(1. - q_h0) * pow((a + vp) / Rv, 2) +
+                2. * (1. - q_h0) * (a + vp) / Rv + q_h0;
+    if ((a + vp) > Rv)
+        Rpel = 1.0;
+
+    const auto q_h = Rpel / (pow(aux_arg, aux_ex));
+
+    Fel = power<2>((1. - q_h) * power<2>(aux) + seq / f_c) +
+          (m_0 * aux - 1) * power<2>(q_h);
+}
+
+@ComputeStress
+{
+    sig = (1 - dam) * (lambda * trace(eel) * Stensor::Id() + 2 * mu * eel);
+}
+
+@Integrator
+{
+    auto const T_Celsius = T - 273.15;
+
+    Stensor nvp = Stensor(0.);
+    Stensor s_dev = deviator(sig / (1 - dam));
+    const stress seq = sigmaeq(s_dev);
+
+    if (seq > 1.e-15)
+    {
+        nvp = 1.5 * s_dev / seq;
+    }
+    fvp -=
+        dt * A_creep * exp(-Qact / 8.3144598 / T_Celsius) * pow(seq, n_creep);
+
+    if (Fel > 0)
+    {
+        Stensor nq = s_dev;
+        Stensor np = 1.0 * Stensor::Id();
+        const auto aux = (seq + trace(sig / (1 - dam))) / (3. * f_c);
+
+        const auto aux_ex = (1. - 1. / n_exp_T);
+        const auto aux_arg = 1 + pow(alpha_p * (T_Celsius - 10), n_exp_T);
+        real a_;
+        real vp_;
+
+        a_ = a + theta * da;
+        vp_ = vp + theta * dvp;
+        Rv = max(at_1 * T_Celsius + at_2 * trace(sig / (1 - dam)) / 3. + at_3,
+                 1.0e-4);
+        auto Rp_ = -(1. - q_h0) * pow((a_ + vp_) / Rv, 2) +
+                   2. * (1. - q_h0) * (a_ + vp_) / Rv + q_h0;
+        if ((a_ + vp_) > Rv)
+            Rp_ = 1.0;
+        const auto q_h = Rp_ / (pow(aux_arg, aux_ex));
+
+        const auto yield = power<2>((1. - q_h) * power<2>(aux) + seq / f_c) +
+                           (m_0 * aux - 1) * power<2>(q_h);
+
+        const auto big_aux = power<2>(aux) * (1 - q_h) + seq / f_c;
+        const auto dev_flow =
+            power<2>(q_h) * m_0 / (2. * f_c * seq) +
+            2. * ((1 - q_h) * aux / (f_c * seq) + 3. / (2. * f_c * seq)) *
+                big_aux;
+        const auto iso_flow = power<2>(q_h) * m_0 / (3. * f_c) +
+                              4. * (1. - q_h) * aux / (3. * f_c) * big_aux;
+
+        feel = deel - deto + dpla_mult * (nq * dev_flow + np * iso_flow) +
+               dvp * nvp;
+        fpla_mult = yield / young;
+        fa = da - sqrt(2. / 3. * (dpla_mult * (nq * dev_flow + np * iso_flow)) |
+                       (dpla_mult * (nq * dev_flow + np * iso_flow)));
+        triax_p = iso_flow / dev_flow;
+
+        dam = max(min(1 - exp(-(((a + vp) - Rv) / alpha_d)), 1.0), 0.0);
+        sig = (1 - dam) * sig;
+    }
+    else
+    {
+        feel = deel - deto + dvp * nvp;
+    }
+}
+
+@TangentOperator
+{
+    if ((smt == ELASTIC) || (smt == SECANTOPERATOR))
+    {
+        computeElasticStiffness<N, Type>::exe(Dt, lambda, mu);
+    }
+    else if (smt == CONSISTENTTANGENTOPERATOR)
+    {
+        StiffnessTensor De;
+        Stensor4 Je;
+        computeElasticStiffness<N, Type>::exe(De, lambda, mu);
+        getPartialJacobianInvert(Je);
+        Dt = De * Je;
+    }
+    else
+    {
+        return false;
+    }
+}
diff --git a/Tests/Data/ThermoMechanics/BDT/bdt.mtest b/MaterialLib/SolidModels/MFront/BDT.mtest
similarity index 96%
rename from Tests/Data/ThermoMechanics/BDT/bdt.mtest
rename to MaterialLib/SolidModels/MFront/BDT.mtest
index d2c5b2b1a296aaa40abcf1f041b2a27391ad0b2e..4b058e896512b0e6729dd67e0c945ea61c21eafe 100644
--- a/Tests/Data/ThermoMechanics/BDT/bdt.mtest
+++ b/MaterialLib/SolidModels/MFront/BDT.mtest
@@ -1,4 +1,4 @@
-@Behaviour<generic> 'src/libBehaviour.so' 'BDT';
+@Behaviour<generic> 'libOgsMFrontBehaviour.so' 'BDT';
 @MaterialProperty<constant> 'YoungModulus' 4.4619E+04;
 @MaterialProperty<constant> 'PoissonRatio' 3.0000E-01;
 @MaterialProperty<constant> 'ThermalExpansion' 1.0000E-05;
diff --git a/MaterialLib/SolidModels/MFront/CMakeLists.txt b/MaterialLib/SolidModels/MFront/CMakeLists.txt
index 05b2a276d8e8d58fe98966a0899431742a2d5554..54df69a8d44b3b6ac41bde131ebbb00dc5e7c38a 100644
--- a/MaterialLib/SolidModels/MFront/CMakeLists.txt
+++ b/MaterialLib/SolidModels/MFront/CMakeLists.txt
@@ -1,8 +1,6 @@
 set(SOURCES CreateMFront.cpp CreateMFront.h)
 
-if(OGS_USE_MFRONT)
-    list(APPEND SOURCES MFront.cpp MFront.h)
-endif()
+list(APPEND SOURCES MFront.cpp MFront.h)
 
 add_library(MaterialLib_SolidModels_MFront ${SOURCES})
 
@@ -12,13 +10,21 @@ if(BUILD_SHARED_LIBS)
 endif()
 
 target_link_libraries(MaterialLib_SolidModels_MFront
-                      PUBLIC BaseLib NumLib logog
+                      PUBLIC BaseLib NumLib logog OgsMFrontBehaviour
                       PRIVATE MathLib MeshLib)
 
-if(OGS_USE_MFRONT)
-    target_include_directories(MaterialLib_SolidModels_MFront
-                               PUBLIC ${MGIS_INCLUDE_DIR})
-    target_link_libraries(MaterialLib_SolidModels_MFront PUBLIC ${MGIS_LIBRARY})
-    target_compile_definitions(MaterialLib_SolidModels_MFront PRIVATE
-                               OGS_USE_MFRONT)
-endif()
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH}
+                      "${PROJECT_SOURCE_DIR}/ThirdParty/MGIS")
+include(cmake/modules/tfel)
+mfront_behaviours_check_library(OgsMFrontBehaviour
+                                BDT
+                                DruckerPrager
+                                Elasticity
+                                StandardElasticityBrick)
+
+target_include_directories(MaterialLib_SolidModels_MFront
+                           PUBLIC ThirdParty/MGIS/include)
+target_link_libraries(MaterialLib_SolidModels_MFront
+                      PUBLIC MFrontGenericInterface)
+target_compile_definitions(MaterialLib_SolidModels_MFront PRIVATE
+                           OGS_USE_MFRONT)
diff --git a/MaterialLib/SolidModels/MFront/CreateMFront.cpp b/MaterialLib/SolidModels/MFront/CreateMFront.cpp
index 97446167c46e28f6e5408233cdc9f3806f12bb17..8bfc0e88870ceee0cc780279ff333e7238b7325a 100644
--- a/MaterialLib/SolidModels/MFront/CreateMFront.cpp
+++ b/MaterialLib/SolidModels/MFront/CreateMFront.cpp
@@ -9,8 +9,6 @@
 
 #include "CreateMFront.h"
 
-#ifdef OGS_USE_MFRONT
-
 #include "BaseLib/FileTools.h"
 #include "ParameterLib/Utils.h"
 
@@ -35,6 +33,17 @@ void varInfo(std::string const& msg,
              mgis::behaviour::getVariableOffset(vars, var.name, hypothesis));
     }
 }
+
+/// Prints info about MFront parameters.
+void varInfo(std::string const& msg, std::vector<std::string> const& parameters)
+{
+    INFO("#%s: %lu (array size %lu).", msg.c_str(), parameters.size());
+    // mgis::behaviour::getArraySize(vars, hypothesis));
+    for (auto const& parameter : parameters)
+    {
+        INFO("  --> with name `%s'.", parameter.c_str());
+    }
+}
 }  // anonymous namespace
 
 namespace MaterialLib
@@ -53,10 +62,14 @@ std::unique_ptr<MechanicsBase<DisplacementDim>> createMFront(
     //! \ogs_file_param{material__solid__constitutive_relation__type}
     config.checkConfigParameter("type", "MFront");
 
-    auto const lib_path = BaseLib::joinPaths(
-        BaseLib::getProjectDirectory(),
-        //! \ogs_file_param{material__solid__constitutive_relation__MFront__library}
-        config.getConfigParameter<std::string>("library"));
+    //! \ogs_file_param{material__solid__constitutive_relation__MFront__library}
+    auto const library_name =
+        config.getConfigParameterOptional<std::string>("library");
+    auto const lib_path =
+        library_name
+            ? BaseLib::joinPaths(BaseLib::getProjectDirectory(), *library_name)
+            : "libOgsMFrontBehaviour.so";
+
     auto const behaviour_name =
         //! \ogs_file_param{material__solid__constitutive_relation__MFront__behaviour}
         config.getConfigParameter<std::string>("behaviour");
@@ -99,7 +112,9 @@ std::unique_ptr<MechanicsBase<DisplacementDim>> createMFront(
 
     // TODO read parameters from prj file, not yet (2018-11-05) supported by
     // MGIS library.
-    varInfo("Parameters", behaviour.parameters, hypothesis);
+    varInfo("Real-valued parameters", behaviour.params);
+    varInfo("Integer parameters", behaviour.iparams);
+    varInfo("Unsigned parameters", behaviour.usparams);
 
     std::vector<ParameterLib::Parameter<double> const*> material_properties;
 
@@ -182,28 +197,6 @@ std::unique_ptr<MechanicsBase<DisplacementDim>> createMFront(
 }  // namespace Solids
 }  // namespace MaterialLib
 
-#else  // OGS_USE_MFRONT
-
-namespace MaterialLib
-{
-namespace Solids
-{
-namespace MFront
-{
-template <int DisplacementDim>
-std::unique_ptr<MechanicsBase<DisplacementDim>> createMFront(
-    std::vector<
-        std::unique_ptr<ParameterLib::ParameterBase>> const& /*parameters*/,
-    BaseLib::ConfigTree const& /*config*/)
-{
-    OGS_FATAL("OpenGeoSys has not been build with MFront support.");
-}
-}  // namespace MFront
-}  // namespace Solids
-}  // namespace MaterialLib
-
-#endif  // OGS_USE_MFRONT
-
 namespace MaterialLib
 {
 namespace Solids
diff --git a/MaterialLib/SolidModels/MFront/DruckerPrager.mfront b/MaterialLib/SolidModels/MFront/DruckerPrager.mfront
new file mode 100644
index 0000000000000000000000000000000000000000..5f138f4ea68b1ed99fe66288d2b182b4c05c630b
--- /dev/null
+++ b/MaterialLib/SolidModels/MFront/DruckerPrager.mfront
@@ -0,0 +1,83 @@
+@Parser Implicit;
+@Behaviour DruckerPrager;
+@Algorithm NewtonRaphson_NumericalJacobian;
+@Theta 1.;
+
+@MaterialProperty real young;
+young.setGlossaryName("YoungModulus");
+@MaterialProperty real nu;
+nu.setGlossaryName("PoissonRatio");
+@MaterialProperty real kappa;
+@MaterialProperty real beta;
+
+@StateVariable real pla_mult;
+@StateVariable real a;
+
+@LocalVariable stress lambda;
+@LocalVariable stress mu;
+@LocalVariable stress Fel;
+@LocalVariable StressStensor s0;
+
+@InitLocalVariables
+{
+    lambda = computeLambda(young, nu);
+    mu = computeMu(young, nu);
+    StressStensor sigel(lambda * trace(eel + deto) * Stensor::Id() +
+                        2 * mu * (eel + deto));
+    const auto s_dev = deviator(sigel);
+    const stress J2 = (s_dev | s_dev) / 2.;
+    const stress I1 = trace(sigel);
+
+    Fel = sqrt(J2) + beta * I1 - kappa;
+}
+
+@ComputeStress
+{
+    sig = (lambda * trace(eel) * Stensor::Id() + 2 * mu * eel);
+}
+
+@Integrator
+{
+    Stensor s_dev = deviator(sig);
+
+    if (Fel > 0)
+    {
+        Stensor nq = sqrt(3.) / 2. / sigmaeq(sig) * s_dev;
+        Stensor np = 1.0 * Stensor::Id();
+
+        const auto yield =
+            sqrt((s_dev | s_dev) / 2.) + beta * trace(sig) - kappa;
+
+        const auto dev_flow = 1.0;
+        const auto iso_flow = beta;
+
+        feel = deel - deto + dpla_mult * (nq * dev_flow + np * iso_flow);
+        fpla_mult = yield / young;
+        fa = da - sqrt(2. / 3. * (dpla_mult * (nq * dev_flow + np * iso_flow)) |
+                       (dpla_mult * (nq * dev_flow + np * iso_flow)));
+    }
+    else
+    {
+        feel = deel - deto;
+    }
+}
+
+@TangentOperator
+{
+    if ((smt == ELASTIC) || (smt == SECANTOPERATOR))
+    {
+        computeElasticStiffness<N, Type>::exe(Dt, lambda, mu);
+    }
+    else if (smt == CONSISTENTTANGENTOPERATOR)
+    {
+        StiffnessTensor De;
+        Stensor4 Je;
+        computeElasticStiffness<N, Type>::exe(De, lambda, mu);
+        getPartialJacobianInvert(Je);
+        Dt = De * Je;
+    }
+    else
+    {
+        return false;
+    }
+}
diff --git a/MaterialLib/SolidModels/MFront/Elasticity.mfront b/MaterialLib/SolidModels/MFront/Elasticity.mfront
new file mode 100644
index 0000000000000000000000000000000000000000..b257b785756efae7393e42d147fdd75e0e419ca2
--- /dev/null
+++ b/MaterialLib/SolidModels/MFront/Elasticity.mfront
@@ -0,0 +1,36 @@
+@Behaviour Elasticity;
+@Author Helfer Thomas;
+@Date 23/11/06;
+@Description{
+    A very first test
+    (the most simple one in fact).
+}
+
+@UseQt true;
+@ProvidesSymmetricTangentOperator;
+
+@MaterialProperty stress young;
+young.setGlossaryName("YoungModulus");
+@MaterialProperty real nu;
+nu.setGlossaryName("PoissonRatio");
+
+@LocalVariable stress lambda;
+@LocalVariable stress mu;
+
+@Includes
+{
+#include "TFEL/Material/Lame.hxx"
+}
+
+@Integrator
+{
+    using namespace tfel::material::lame;
+    lambda = computeLambda(young, nu);
+    mu = computeMu(young, nu);
+    sig = lambda * trace(eto + deto) * StrainStensor::Id() +
+          2 * mu * (eto + deto);
+    if (computeTangentOperator_)
+    {
+        Dt = lambda * Stensor4::IxI() + 2 * mu * Stensor4::Id();
+    }
+}
diff --git a/MaterialLib/SolidModels/MFront/MFront.cpp b/MaterialLib/SolidModels/MFront/MFront.cpp
index 95e076e543a7471881fbb3b5925015cfd1f5f52f..3051d48f936b7a937be052821e9d6b4b90acc85d 100644
--- a/MaterialLib/SolidModels/MFront/MFront.cpp
+++ b/MaterialLib/SolidModels/MFront/MFront.cpp
@@ -222,6 +222,18 @@ MFront<DisplacementDim>::MFront(
                 mgis::behaviour::getVariableSize(
                     _behaviour.thermodynamic_forces[0], hypothesis));
     }
+
+    if (_behaviour.mps.size() != _material_properties.size())
+    {
+        ERR("There are %d material properties in the loaded behaviour:",
+            _behaviour.mps.size());
+        for (auto const& mp : _behaviour.mps)
+        {
+            ERR("\t%s", mp.name.c_str());
+        }
+        OGS_FATAL("But the number of passed material properties is %d.",
+                  _material_properties.size());
+    }
 }
 
 template <int DisplacementDim>
@@ -249,33 +261,35 @@ MFront<DisplacementDim>::integrateStress(
 {
     assert(
         dynamic_cast<MaterialStateVariables const*>(&material_state_variables));
-    auto& d =
-        static_cast<MaterialStateVariables const&>(material_state_variables)
-            ._data;
+    // New state, copy of current one, packed in unique_ptr for return.
+    auto state = std::make_unique<MaterialStateVariables>(
+        static_cast<MaterialStateVariables const&>(material_state_variables));
+    auto& behaviour_data = state->_behaviour_data;
 
     // TODO add a test of material behaviour where the value of dt matters.
-    d.dt = dt;
-    d.rdt = 1.0;
-    d.K[0] = 4.0;  // if K[0] is greater than 3.5, the consistent tangent
-                   // operator must be computed.
+    behaviour_data.dt = dt;
+    behaviour_data.rdt = 1.0;
+    behaviour_data.K[0] = 4.0;  // if K[0] is greater than 3.5, the consistent
+                                // tangent operator must be computed.
 
     // evaluate parameters at (t, x)
     {
-        auto out = d.s1.material_properties.begin();
+        auto out = behaviour_data.s1.material_properties.begin();
         for (auto* param : _material_properties)
         {
             auto const& vals = (*param)(t, x);
             out = std::copy(vals.begin(), vals.end(), out);
         }
+        assert(out == behaviour_data.s1.material_properties.end());
     }
 
-    if (!d.s1.external_state_variables.empty())
+    if (!behaviour_data.s1.external_state_variables.empty())
     {
         // assuming that there is only temperature
-        d.s1.external_state_variables[0] = T;
+        behaviour_data.s1.external_state_variables[0] = T;
     }
 
-    auto v = mgis::behaviour::make_view(d);
+    auto v = mgis::behaviour::make_view(behaviour_data);
 
     auto const eps_MFront = OGSToMFront(eps);
     for (auto i = 0; i < KelvinVector::SizeAtCompileTime; ++i)
@@ -292,26 +306,66 @@ MFront<DisplacementDim>::integrateStress(
     KelvinVector sigma;
     for (auto i = 0; i < KelvinVector::SizeAtCompileTime; ++i)
     {
-        sigma[i] = d.s1.thermodynamic_forces[i];
+        sigma[i] = behaviour_data.s1.thermodynamic_forces[i];
     }
     sigma = MFrontToOGS(sigma);
 
     // TODO row- vs. column-major storage order. This should only matter for
     // anisotropic materials.
-    if (d.K.size() !=
+    if (behaviour_data.K.size() !=
         KelvinMatrix::RowsAtCompileTime * KelvinMatrix::ColsAtCompileTime)
         OGS_FATAL("Stiffness matrix has wrong size.");
 
-    KelvinMatrix C = MFrontToOGS(Eigen::Map<KelvinMatrix>(d.K.data()));
+    KelvinMatrix C =
+        MFrontToOGS(Eigen::Map<KelvinMatrix>(behaviour_data.K.data()));
+
+    return boost::make_optional(
+        std::make_tuple<typename MFront<DisplacementDim>::KelvinVector,
+                        std::unique_ptr<typename MechanicsBase<
+                            DisplacementDim>::MaterialStateVariables>,
+                        typename MFront<DisplacementDim>::KelvinMatrix>(
+            std::move(sigma),
+            std::move(state),
+            std::move(C)));
+}
 
-    // TODO avoid copying the state
-    auto state_copy = std::make_unique<MaterialStateVariables>(
-        static_cast<MaterialStateVariables const&>(material_state_variables));
-    std::unique_ptr<
-        typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
-        state_upcast(state_copy.release());
+template <int DisplacementDim>
+std::vector<typename MechanicsBase<DisplacementDim>::InternalVariable>
+MFront<DisplacementDim>::getInternalVariables() const
+{
+    std::vector<typename MechanicsBase<DisplacementDim>::InternalVariable>
+        internal_variables;
+
+    for (auto const& iv : _behaviour.isvs)
+    {
+        auto const name = iv.name;
+        auto const offset = mgis::behaviour::getVariableOffset(
+            _behaviour.isvs, name, _behaviour.hypothesis);
+        auto const size =
+            mgis::behaviour::getVariableSize(iv, _behaviour.hypothesis);
+
+        typename MechanicsBase<DisplacementDim>::InternalVariable new_variable{
+            name, static_cast<unsigned>(size),
+            [offset, size](
+                typename MechanicsBase<
+                    DisplacementDim>::MaterialStateVariables const& state,
+                std::vector<double>& cache) -> std::vector<double> const& {
+                assert(dynamic_cast<MaterialStateVariables const*>(&state) !=
+                       nullptr);
+                auto const& internal_state_variables =
+                    static_cast<MaterialStateVariables const&>(state)
+                        ._behaviour_data.s1.internal_state_variables;
+
+                cache.resize(size);
+                std::copy_n(internal_state_variables.data() + offset,
+                            size,
+                            begin(cache));
+                return cache;
+            }};
+        internal_variables.push_back(new_variable);
+    }
 
-    return {std::make_tuple(std::move(sigma), std::move(state_upcast), std::move(C))};
+    return internal_variables;
 }
 
 template <int DisplacementDim>
diff --git a/MaterialLib/SolidModels/MFront/MFront.h b/MaterialLib/SolidModels/MFront/MFront.h
index 5f886b22d690d0c33a80f17ce300c8f89e4a6afc..2afe85184c91cd211cf43fae0a106c142d7894a6 100644
--- a/MaterialLib/SolidModels/MFront/MFront.h
+++ b/MaterialLib/SolidModels/MFront/MFront.h
@@ -44,26 +44,19 @@ public:
         : public MechanicsBase<DisplacementDim>::MaterialStateVariables
     {
         explicit MaterialStateVariables(mgis::behaviour::Behaviour const& b)
-            : _data{b}
+            : _behaviour_data{b}
         {
         }
 
         MaterialStateVariables(MaterialStateVariables const&) = default;
+        MaterialStateVariables(MaterialStateVariables&&) = delete;
 
-        void pushBackState() override { mgis::behaviour::update(_data); }
-
-        MaterialStateVariables& operator=(MaterialStateVariables const&) =
-            default;
-
-        typename MechanicsBase<DisplacementDim>::MaterialStateVariables&
-        operator=(typename MechanicsBase<DisplacementDim>::
-                      MaterialStateVariables const& state) noexcept override
+        void pushBackState() override
         {
-            return operator=(static_cast<MaterialStateVariables const&>(state));
+            mgis::behaviour::update(_behaviour_data);
         }
 
-        // TODO fix: this has to be mutable.
-        mutable mgis::behaviour::BehaviourData _data;
+        mgis::behaviour::BehaviourData _behaviour_data;
     };
 
     using KelvinVector =
@@ -94,6 +87,9 @@ public:
             material_state_variables,
         double const T) const override;
 
+    std::vector<typename MechanicsBase<DisplacementDim>::InternalVariable>
+    getInternalVariables() const override;
+
     double computeFreeEnergyDensity(
         double const t,
         ParameterLib::SpatialPosition const& x,
diff --git a/MaterialLib/SolidModels/MFront/StandardElasticityBrick.mfront b/MaterialLib/SolidModels/MFront/StandardElasticityBrick.mfront
new file mode 100644
index 0000000000000000000000000000000000000000..c8005f444e28b882da1d55937e52d50d4b45645d
--- /dev/null
+++ b/MaterialLib/SolidModels/MFront/StandardElasticityBrick.mfront
@@ -0,0 +1,24 @@
+@DSL Implicit;
+@Behaviour StandardElasticityBrick;
+@Author Thomas Nagel;
+@Date 05/02/2019;
+
+@Algorithm NewtonRaphson;
+//@CompareToNumericalJacobian true;
+//@NumericallyComputedJacobianBlocks {dfeel_ddeel};
+// remove the above blocks once an analytical one is provided.
+//@PerturbationValueForNumericalJacobianComputation 1.e-7;
+//@JacobianComparisonCriterion 1.e-6; // adjust to your needs
+
+
+@Brick StandardElasticity;
+
+@Theta 1.;       // time integration scheme
+@Epsilon 1e-14;  // tolerance of local stress integration algorithm
+@ModellingHypotheses{".+"};
+
+@RequireStiffnessTensor<UnAltered>;
+
+@Integrator
+{
+}
diff --git a/ProcessLib/SmallDeformation/Tests.cmake b/ProcessLib/SmallDeformation/Tests.cmake
index fdb2896cf040e780e20170ac323ea90b1618440a..8f7010f8aa1bdfd0fac7fa1e0ad95ae665b9cab4 100644
--- a/ProcessLib/SmallDeformation/Tests.cmake
+++ b/ProcessLib/SmallDeformation/Tests.cmake
@@ -73,7 +73,7 @@ AddTest(
     # See also the prj file.
     DIFF_DATA
     cube_1e0_dp_ref_created_with_OGS_Ehlers.vtu cube_1e0_dp_pcs_0_ts_203_t_5.100000.vtu displacement displacement 1e-14 0
-    cube_1e0_dp_ref_created_with_OGS_Ehlers.vtu cube_1e0_dp_pcs_0_ts_203_t_5.100000.vtu sigma sigma 1e-13 0
+    cube_1e0_dp_ref_created_with_OGS_Ehlers.vtu cube_1e0_dp_pcs_0_ts_203_t_5.100000.vtu sigma sigma 2e-13 0
     cube_1e0_dp_ref_created_with_OGS_Ehlers.vtu cube_1e0_dp_pcs_0_ts_203_t_5.100000.vtu epsilon epsilon 1e-14 0
 )
 
@@ -87,8 +87,8 @@ AddTest(
     TESTER vtkdiff
     REQUIREMENTS OGS_USE_MFRONT AND NOT OGS_USE_MPI
     DIFF_DATA
-    ../../ring_plane_strain_1e4_solution.vtu ring_plane_strain_pcs_0_ts_1_t_1.000000.vtu displacement displacement 1e-16 0
-    ../../ring_plane_strain_1e4_solution.vtu ring_plane_strain_pcs_0_ts_1_t_1.000000.vtu sigma sigma 1e-15 0
+    ../../ring_plane_strain_pcs_0_ts_1_t_1.000000.vtu ring_plane_strain_pcs_0_ts_1_t_1.000000.vtu displacement displacement 1e-16 0
+    ../../ring_plane_strain_pcs_0_ts_1_t_1.000000.vtu ring_plane_strain_pcs_0_ts_1_t_1.000000.vtu sigma sigma 1e-15 0
 )
 
 AddTest(
@@ -177,4 +177,4 @@ AddTest(
     REQUIREMENTS NOT OGS_USE_MPI
     DIFF_DATA
     m1_3Dtopload_pcs_0_ts_1_t_1.000000.vtu m1_3Dtopload_pcs_0_ts_1_t_1.000000.vtu displacement displacement 10e-12 0.0
-)
\ No newline at end of file
+)
diff --git a/ProcessLib/ThermoMechanics/Tests.cmake b/ProcessLib/ThermoMechanics/Tests.cmake
index d91d9d238a630c8638c1e5e6fdccfb0e699b45e2..af6efa001f8e5b7bc1a4174befb4fbd3aa2599de 100644
--- a/ProcessLib/ThermoMechanics/Tests.cmake
+++ b/ProcessLib/ThermoMechanics/Tests.cmake
@@ -331,14 +331,14 @@ AddTest(
     TESTER vtkdiff
     REQUIREMENTS OGS_USE_MFRONT AND NOT (OGS_USE_LIS OR OGS_USE_MPI)
     DIFF_DATA
-    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_51_t_300.000000.vtu     epsilon_300   epsilon  1e-8   0
-    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_51_t_300.000000.vtu     sigma_300     sigma    1e-3   0
-    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_151_t_900.000000.vtu    epsilon_900   epsilon  1e-9   0
-    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_151_t_900.000000.vtu    sigma_900     sigma    1e-3   0
-    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_251_t_1500.000000.vtu   epsilon_1500  epsilon  1e-8   0
-    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_251_t_1500.000000.vtu   sigma_1500    sigma    1e-3   0
-    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_501_t_3000.000000.vtu   epsilon_3000  epsilon  1e-10  0
-    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_501_t_3000.000000.vtu   sigma_3000    sigma    1e-3   0
-    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_1001_t_6000.000000.vtu  epsilon_6000  epsilon  1e-7   0
-    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_1001_t_6000.000000.vtu  sigma_6000    sigma    1e-3   0
+    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_51_t_300.000000.vtu     epsilon_300   epsilon  1e-4   0
+    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_51_t_300.000000.vtu     sigma_300     sigma    2e+1   0
+    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_151_t_900.000000.vtu    epsilon_900   epsilon  1e-4   0
+    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_151_t_900.000000.vtu    sigma_900     sigma    1e+0   0
+    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_251_t_1500.000000.vtu   epsilon_1500  epsilon  1e-4   0
+    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_251_t_1500.000000.vtu   sigma_1500    sigma    1e+0   0
+    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_501_t_3000.000000.vtu   epsilon_3000  epsilon  1e-4  0
+    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_501_t_3000.000000.vtu   sigma_3000    sigma    1e+0   0
+    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_1001_t_6000.000000.vtu  epsilon_6000  epsilon  1e-4   0
+    bdt_ref.vtu  cube_1e0_bdt_pcs_0_ts_1001_t_6000.000000.vtu  sigma_6000    sigma    1e+0   0
 )
diff --git a/Tests/Data/Mechanics/Ehlers/MFront/DruckerPrager.mfront b/Tests/Data/Mechanics/Ehlers/MFront/DruckerPrager.mfront
deleted file mode 100644
index b6a5571748a8b3ed43b1f3b402da321050af8aab..0000000000000000000000000000000000000000
--- a/Tests/Data/Mechanics/Ehlers/MFront/DruckerPrager.mfront
+++ /dev/null
@@ -1,72 +0,0 @@
-@Parser Implicit;
-@Behaviour DruckerPrager;
-@Algorithm NewtonRaphson_NumericalJacobian;
-@Theta 1.;
-
-@MaterialProperty real young;
-young.setGlossaryName("YoungModulus");
-@MaterialProperty real nu;
-nu.setGlossaryName("PoissonRatio");
-@MaterialProperty real kappa;
-@MaterialProperty real beta;
-
-@StateVariable real pla_mult;
-@StateVariable real a;
-
-@LocalVariable stress lambda;
-@LocalVariable stress mu;
-@LocalVariable stress Fel;
-@LocalVariable StressStensor s0;
-
-@InitLocalVariables {
-  lambda = computeLambda(young, nu);
-  mu = computeMu(young, nu);
-  StressStensor sigel(lambda * trace(eel + deto) * Stensor::Id() +
-                      2 * mu * (eel + deto));
-  const auto s_dev = deviator(sigel);
-  const stress J2 = (s_dev|s_dev)/2.;
-  const stress I1 = trace(sigel);
-
-  Fel = sqrt(J2) + beta * I1 - kappa;
-}
-
-@ComputeStress {
-  sig = (lambda * trace(eel) * Stensor::Id() + 2 * mu * eel);
-}
-
-@Integrator {
-  Stensor s_dev = deviator(sig);
-
-  if (Fel > 0) {
-    Stensor nq = sqrt(3.)/2./sigmaeq(sig)*s_dev;
-    Stensor np = 1.0 * Stensor::Id();
-
-    const auto yield = sqrt((s_dev|s_dev)/2.) + beta * trace(sig) - kappa;
-
-    const auto dev_flow = 1.0;
-    const auto iso_flow = beta;
-
-    feel =
-        deel - deto + dpla_mult * (nq * dev_flow + np * iso_flow);
-    fpla_mult = yield / young;
-    fa = da - sqrt(2. / 3. * (dpla_mult * (nq * dev_flow + np * iso_flow)) |
-                   (dpla_mult * (nq * dev_flow + np * iso_flow)));
-
-  } else {
-    feel = deel - deto ;
-  }
-}
-
-@TangentOperator {
-  if ((smt == ELASTIC) || (smt == SECANTOPERATOR)) {
-    computeElasticStiffness<N, Type>::exe(Dt, lambda, mu);
-  } else if (smt == CONSISTENTTANGENTOPERATOR) {
-    StiffnessTensor De;
-    Stensor4 Je;
-    computeElasticStiffness<N, Type>::exe(De, lambda, mu);
-    getPartialJacobianInvert(Je);
-    Dt = De * Je;
-  } else {
-    return false;
-  }
-}
diff --git a/Tests/Data/Mechanics/Ehlers/MFront/DruckerPrager.so b/Tests/Data/Mechanics/Ehlers/MFront/DruckerPrager.so
deleted file mode 100755
index a311106f008e899cd4c8fe9eb3146de01d08ee50..0000000000000000000000000000000000000000
Binary files a/Tests/Data/Mechanics/Ehlers/MFront/DruckerPrager.so and /dev/null differ
diff --git a/Tests/Data/Mechanics/Ehlers/MFront/cube_1e0_dp.prj b/Tests/Data/Mechanics/Ehlers/MFront/cube_1e0_dp.prj
index 936c3b4fbad43f6ce2affa1bf0feccb39afe8c05..37a43241e6dd7abe8426e7e5b183e3fcfecbf0e7 100644
--- a/Tests/Data/Mechanics/Ehlers/MFront/cube_1e0_dp.prj
+++ b/Tests/Data/Mechanics/Ehlers/MFront/cube_1e0_dp.prj
@@ -9,7 +9,6 @@
             <integration_order>2</integration_order>
             <constitutive_relation>
                 <type>MFront</type>
-                <library>DruckerPrager.so</library>
                 <behaviour>DruckerPrager</behaviour>
                 <material_properties>
                     <material_property name="YoungModulus" parameter="E"/>
diff --git a/Tests/Data/Mechanics/Linear/MFront/axisymm_ring/ring_plane_strain.prj b/Tests/Data/Mechanics/Linear/MFront/axisymm_ring/ring_plane_strain.prj
index d7b240863fa18011b2f0d7614626ceb45c3dd42f..c0eaff4ea8c9968008798916de08ea76e94bba3d 100644
--- a/Tests/Data/Mechanics/Linear/MFront/axisymm_ring/ring_plane_strain.prj
+++ b/Tests/Data/Mechanics/Linear/MFront/axisymm_ring/ring_plane_strain.prj
@@ -9,7 +9,6 @@
             <integration_order>2</integration_order>
             <constitutive_relation>
                 <type>MFront</type>
-                <library>../disc_with_hole/libBehaviour.so</library>
                 <behaviour>Elasticity</behaviour>
                 <material_properties>
                     <material_property name="YoungModulus" parameter="E"/>
diff --git a/Tests/Data/Mechanics/Linear/MFront/cube_1x1x1.gml b/Tests/Data/Mechanics/Linear/MFront/cube_1x1x1.gml
new file mode 100644
index 0000000000000000000000000000000000000000..92bb259f46c3e9d4d808fd638556cfa1fd337e77
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/MFront/cube_1x1x1.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:eb8c02563eaa7fc99ab48b12feef92962f2265821911831f05ced5673847d7a4
+size 1602
diff --git a/Tests/Data/Mechanics/Linear/MFront/cube_1x1x1_hex_1e0.vtu b/Tests/Data/Mechanics/Linear/MFront/cube_1x1x1_hex_1e0.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..7e182287a44b78741a25cbcfabc2e0109d918d54
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/MFront/cube_1x1x1_hex_1e0.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:baa51eea253d1fe35591bc41b3723e19f1102bb6812ecef5d6aaa5766a53cdc7
+size 1437
diff --git a/Tests/Data/Mechanics/Linear/MFront/cube_standard_elasticity_brick.prj b/Tests/Data/Mechanics/Linear/MFront/cube_standard_elasticity_brick.prj
new file mode 100644
index 0000000000000000000000000000000000000000..fd8e8f9eb0dedcd3c970bb9b50d244b6ee250da7
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/MFront/cube_standard_elasticity_brick.prj
@@ -0,0 +1,193 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <mesh>cube_1x1x1_hex_1e0.vtu</mesh>
+    <geometry>cube_1x1x1.gml</geometry>
+    <processes>
+        <process>
+            <name>SD</name>
+            <type>SMALL_DEFORMATION</type>
+            <integration_order>2</integration_order>
+            <constitutive_relation>
+                <type>MFront</type>
+                <behaviour>StandardElasticityBrick</behaviour>
+                <material_properties>
+                    <material_property name="YoungModulus" parameter="E"/>
+                    <material_property name="PoissonRatio" parameter="nu"/>
+                </material_properties>
+            </constitutive_relation>
+            <solid_density>rho_sr</solid_density>
+            <specific_body_force>0 0 0</specific_body_force>
+            <reference_temperature>293.15</reference_temperature>
+            <process_variables>
+                <process_variable>displacement</process_variable>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+            </secondary_variables>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="SD">
+                <nonlinear_solver>basic_newton</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstol>1e-14</abstol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0</t_initial>
+                    <t_end>1</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>4</repeat>
+                            <delta_t>0.25</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>cube_1e0</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>1</repeat>
+                    <each_steps>10000000</each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable>displacement</variable>
+                <variable>sigma</variable>
+            </variables>
+        </output>
+    </time_loop>
+    <parameters>
+        <parameter>
+            <name>E</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>nu</name>
+            <type>Constant</type>
+            <value>.3</value>
+        </parameter>
+        <parameter>
+            <name>Cohesion</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>FrictionAngle</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>DilatancyAngle</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>TransitionAngle</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>TensionCutOffParameter</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>rho_sr</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>displacement0</name>
+            <type>Constant</type>
+            <values>0 0 0</values>
+        </parameter>
+        <parameter>
+            <name>dirichlet0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>dirichlet1</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>neumann_force</name>
+            <type>Constant</type>
+            <values>0.01</values>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>3</components>
+            <order>1</order>
+            <initial_condition>displacement0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>front</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>2</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>top</geometry>
+                    <type>Neumann</type>
+                    <component>2</component>
+                    <parameter>neumann_force</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_newton</name>
+            <type>Newton</type>
+            <max_iter>5</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+            <eigen>
+                <solver_type>SparseLU</solver_type>
+                <precon_type>DIAGONAL</precon_type>
+                <max_iteration_step>10000</max_iteration_step>
+                <error_tolerance>1e-16</error_tolerance>
+            </eigen>
+            <petsc>
+                <prefix>sd</prefix>
+                <parameters>-sd_ksp_type cg -sd_pc_type bjacobi -sd_ksp_rtol 1e-16 -sd_ksp_max_it 10000</parameters>
+            </petsc>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Mechanics/Linear/MFront/disc_with_hole/disc_with_hole.prj b/Tests/Data/Mechanics/Linear/MFront/disc_with_hole/disc_with_hole.prj
index 92c23e313d3882fb97de74f35b0a215c02be6203..b1d4633edbf7ae526ae3b18ce67d0916bf7f7ac6 100644
--- a/Tests/Data/Mechanics/Linear/MFront/disc_with_hole/disc_with_hole.prj
+++ b/Tests/Data/Mechanics/Linear/MFront/disc_with_hole/disc_with_hole.prj
@@ -9,7 +9,6 @@
             <integration_order>2</integration_order>
             <constitutive_relation>
                 <type>MFront</type>
-                <library>libBehaviour.so</library>
                 <behaviour>Elasticity</behaviour>
                 <material_properties>
                     <material_property name="YoungModulus" parameter="E"/>
diff --git a/Tests/Data/Mechanics/Linear/MFront/disc_with_hole/elasticity.mfront b/Tests/Data/Mechanics/Linear/MFront/disc_with_hole/elasticity.mfront
deleted file mode 100644
index a985c58315e3a4ffeed3048da9945e4b7f76892c..0000000000000000000000000000000000000000
--- a/Tests/Data/Mechanics/Linear/MFront/disc_with_hole/elasticity.mfront
+++ /dev/null
@@ -1,32 +0,0 @@
-@Behaviour Elasticity;
-@Author Helfer Thomas;
-@Date 23/11/06;
-@Description{
- A very first test
- (the most simple one in fact).
-}
-
-@UseQt true;
-@ProvidesSymmetricTangentOperator;
-
-@MaterialProperty stress young;
-young.setGlossaryName("YoungModulus");
-@MaterialProperty real   nu;
-nu.setGlossaryName("PoissonRatio");
-
-@LocalVariable stress lambda;
-@LocalVariable stress mu;
-
-@Includes{
-#include"TFEL/Material/Lame.hxx"
-}
-
-@Integrator{
-  using namespace tfel::material::lame;
-  lambda = computeLambda(young,nu);
-  mu = computeMu(young,nu);
-  sig = lambda*trace(eto+deto)*StrainStensor::Id()+2*mu*(eto+deto);
-  if(computeTangentOperator_){
-    Dt = lambda*Stensor4::IxI()+2*mu*Stensor4::Id();
-  }
-}
diff --git a/Tests/Data/Mechanics/Linear/MFront/disc_with_hole/libBehaviour.so b/Tests/Data/Mechanics/Linear/MFront/disc_with_hole/libBehaviour.so
deleted file mode 100755
index 93632cd232f2d286a7400d749ef622ae8b81f758..0000000000000000000000000000000000000000
Binary files a/Tests/Data/Mechanics/Linear/MFront/disc_with_hole/libBehaviour.so and /dev/null differ
diff --git a/Tests/Data/ThermoMechanics/BDT/bdt.mfront b/Tests/Data/ThermoMechanics/BDT/bdt.mfront
deleted file mode 100644
index 8ba480df97ec7010735d29d06f577e49dfd3b58a..0000000000000000000000000000000000000000
--- a/Tests/Data/ThermoMechanics/BDT/bdt.mfront
+++ /dev/null
@@ -1,148 +0,0 @@
-@Parser Implicit;
-@Behaviour BDT;
-@Algorithm NewtonRaphson_NumericalJacobian;
-@Theta 1.;
-
-@MaterialProperty real young;
-young.setGlossaryName("YoungModulus");
-@MaterialProperty real nu;
-nu.setGlossaryName("PoissonRatio");
-@MaterialProperty thermalexpansion alpha;
-alpha.setGlossaryName("ThermalExpansion");
-@MaterialProperty real f_c;
-@MaterialProperty real m_0;
-@MaterialProperty real alpha_p;
-@MaterialProperty real n_exp_T;
-@MaterialProperty real q_h0;
-@MaterialProperty real chi_h;
-@MaterialProperty real alpha_d;
-@MaterialProperty real h_d;
-@MaterialProperty real Qact;
-@MaterialProperty real A_creep;
-@MaterialProperty real n_creep;
-@MaterialProperty real El_1;
-@MaterialProperty real El_2;
-@MaterialProperty real El_3;
-@MaterialProperty real at_1;
-@MaterialProperty real at_2;
-@MaterialProperty real at_3;
-
-@StateVariable real pla_mult;
-@StateVariable real a;
-@StateVariable real dam;
-@StateVariable real Rv;
-@StateVariable real triax_p;
-@StateVariable real vp;
-
-@LocalVariable stress lambda;
-@LocalVariable stress mu;
-@LocalVariable stress Fel;
-@LocalVariable StressStensor s0;
-
-// @Parameter Tref;
-// Tref.setDefaultValue(293.15);
-
-@InitLocalVariables {
-  auto const T_Celsius = T - 273.15;
-
-  const auto y2 = El_1 * power<2>(T_Celsius) + El_2 * T_Celsius + El_3;
-  lambda = computeLambda(young * y2, nu);
-  mu = computeMu(young * y2, nu);
-  StressStensor sigel(lambda * trace(eel + deto) * Stensor::Id() +
-                      2 * mu * (eel + deto)); // StresssStensor in Tutorial
-  const auto s_dev = deviator(sigel);
-  // J2 = (s_dev|s_dev)/2.; //defines double contraction
-  const stress seq = sigmaeq(s_dev);
-  const stress I1 = trace(sigel);
-  const auto aux = (seq + I1) / (3. * f_c);
-
-  const auto aux_ex = (1. - 1. / n_exp_T);
-  const auto aux_arg = 1 + pow(alpha_p * (T_Celsius - 10.), n_exp_T);
-
-  Rv = max(at_1 * T_Celsius + at_2 * I1 / 3. + at_3, 1.0e-4);
-  auto Rpel = -(1. - q_h0) * pow((a + vp) / Rv, 2) +
-              2. * (1. - q_h0) * (a + vp) / Rv + q_h0;
-  if ((a + vp) > Rv)
-    Rpel = 1.0;
-
-  const auto q_h = Rpel / (pow(aux_arg, aux_ex));
-
-  Fel = power<2>((1. - q_h) * power<2>(aux) + seq / f_c) +
-        (m_0 * aux - 1) * power<2>(q_h);
-}
-
-@ComputeStress {
-  sig = (1 - dam) * (lambda * trace(eel) * Stensor::Id() + 2 * mu * eel);
-}
-
-@Integrator {
-  auto const T_Celsius = T - 273.15;
-
-  Stensor nvp = Stensor(0.);
-  Stensor s_dev = deviator(sig / (1 - dam));
-  const stress seq = sigmaeq(s_dev);
-
-  if (seq > 1.e-15) {
-    nvp = 1.5 * s_dev / seq;
-  }
-  fvp -= dt * A_creep * exp(-Qact / 8.3144598 / T_Celsius) * pow(seq, n_creep);
-
-  if (Fel > 0) {
-    const real pla_mult_ = pla_mult + theta * dpla_mult;
-    Stensor nq = s_dev;
-    Stensor np = 1.0 * Stensor::Id();
-    const auto aux = (seq + trace(sig / (1 - dam))) / (3. * f_c);
-
-    const auto aux_ex = (1. - 1. / n_exp_T);
-    const auto aux_arg = 1 + pow(alpha_p * (T_Celsius - 10), n_exp_T);
-    real a_;
-    real dam_;
-    real vp_;
-
-    a_ = a + theta * da;
-    dam_ = dam + theta * ddam;
-    vp_ = vp + theta * dvp;
-    Rv = max(at_1 * T_Celsius + at_2 * trace(sig / (1 - dam)) / 3. + at_3, 1.0e-4);
-    auto Rp_ = -(1. - q_h0) * pow((a_ + vp_) / Rv, 2) +
-               2. * (1. - q_h0) * (a_ + vp_) / Rv + q_h0;
-    if ((a_ + vp_) > Rv)
-      Rp_ = 1.0;
-    const auto q_h = Rp_ / (pow(aux_arg, aux_ex));
-
-    const auto yield = power<2>((1. - q_h) * power<2>(aux) + seq / f_c) +
-                       (m_0 * aux - 1) * power<2>(q_h);
-
-    const auto big_aux = power<2>(aux) * (1 - q_h) + seq / f_c;
-    const auto dev_flow =
-        power<2>(q_h) * m_0 / (2. * f_c * seq) +
-        2. * ((1 - q_h) * aux / (f_c * seq) + 3. / (2. * f_c * seq)) * big_aux;
-    const auto iso_flow = power<2>(q_h) * m_0 / (3. * f_c) +
-                          4. * (1. - q_h) * aux / (3. * f_c) * big_aux;
-
-    feel =
-        deel - deto + dpla_mult * (nq * dev_flow + np * iso_flow) + dvp * nvp;
-    fpla_mult = yield / young;
-    fa = da - sqrt(2. / 3. * (dpla_mult * (nq * dev_flow + np * iso_flow)) |
-                   (dpla_mult * (nq * dev_flow + np * iso_flow)));
-    triax_p = iso_flow / dev_flow;
-
-    dam = max(min(1 - exp(-(((a + vp) - Rv) / alpha_d)), 1.0), 0.0);
-    sig = (1 - dam) * sig;
-  } else {
-    feel = deel - deto + dvp * nvp;
-  }
-}
-
-@TangentOperator {
-  if ((smt == ELASTIC) || (smt == SECANTOPERATOR)) {
-    computeElasticStiffness<N, Type>::exe(Dt, lambda, mu);
-  } else if (smt == CONSISTENTTANGENTOPERATOR) {
-    StiffnessTensor De;
-    Stensor4 Je;
-    computeElasticStiffness<N, Type>::exe(De, lambda, mu);
-    getPartialJacobianInvert(Je);
-    Dt = De * Je;
-  } else {
-    return false;
-  }
-}
diff --git a/Tests/Data/ThermoMechanics/BDT/cube_1e0_bdt.prj b/Tests/Data/ThermoMechanics/BDT/cube_1e0_bdt.prj
index 0719b16c569438c1321f7112e1117f9a3560a042..9c4f219694cd179f2361e3365ec6c8a9cbb6065f 100644
--- a/Tests/Data/ThermoMechanics/BDT/cube_1e0_bdt.prj
+++ b/Tests/Data/ThermoMechanics/BDT/cube_1e0_bdt.prj
@@ -9,7 +9,6 @@
             <integration_order>2</integration_order>
             <constitutive_relation>
                 <type>MFront</type>
-                <library>libBehaviour.so</library>
                 <behaviour>BDT</behaviour>
                 <material_properties>
                     <material_property name="YoungModulus" parameter="E"/>
@@ -54,9 +53,9 @@
             <process ref="SD">
                 <nonlinear_solver>basic_newton</nonlinear_solver>
                 <convergence_criterion>
-                    <type>DeltaX</type>
+                    <type>PerComponentDeltaX</type>
                     <norm_type>NORM2</norm_type>
-                    <abstol>1e-14</abstol>
+                    <abstols>1e-16 1e-3 1e-16 1e-3</abstols>
                 </convergence_criterion>
                 <time_discretization>
                     <type>BackwardEuler</type>
diff --git a/Tests/Data/ThermoMechanics/CreepBGRa/Verification/m2_2Dload/m2_2Dload.prj b/Tests/Data/ThermoMechanics/CreepBGRa/Verification/m2_2Dload/m2_2Dload.prj
index 9a084853deff5868946cab4f32f360f1132cd71b..708265a68b92b95a01b32e8e15d964a0f184c53a 100644
--- a/Tests/Data/ThermoMechanics/CreepBGRa/Verification/m2_2Dload/m2_2Dload.prj
+++ b/Tests/Data/ThermoMechanics/CreepBGRa/Verification/m2_2Dload/m2_2Dload.prj
@@ -270,7 +270,7 @@
         <vtkdiff>
             <file>m2_2Dload_pcs_0_ts_6_t_0.500000.vtu</file>
             <field>sigma</field>
-            <absolute_tolerance>1e-11</absolute_tolerance>
+            <absolute_tolerance>2e-11</absolute_tolerance>
             <relative_tolerance>0</relative_tolerance>
         </vtkdiff>
         <vtkdiff>
diff --git a/Tests/Data/ThermoMechanics/LinearMFront/cube_1e0_lin.prj b/Tests/Data/ThermoMechanics/LinearMFront/cube_1e0_lin.prj
index 1641b3dd5698f5f99177c49bf23d5da4962c4d15..cf3dcc9c588b8ce3fc497e154377a559e8d53513 100644
--- a/Tests/Data/ThermoMechanics/LinearMFront/cube_1e0_lin.prj
+++ b/Tests/Data/ThermoMechanics/LinearMFront/cube_1e0_lin.prj
@@ -9,7 +9,6 @@
             <integration_order>2</integration_order>
             <constitutive_relation>
                 <type>MFront</type>
-                <library>libBehaviour.so</library>
                 <behaviour>Elasticity</behaviour>
                 <material_properties>
                     <material_property name="YoungModulus" parameter="E"/>
diff --git a/Tests/Data/ThermoMechanics/LinearMFront/elasticity.mfront b/Tests/Data/ThermoMechanics/LinearMFront/elasticity.mfront
deleted file mode 120000
index 5b5479bafded1f2fccc30bd08d5805ef3dad2aaa..0000000000000000000000000000000000000000
--- a/Tests/Data/ThermoMechanics/LinearMFront/elasticity.mfront
+++ /dev/null
@@ -1 +0,0 @@
-../../Mechanics/Linear/MFront/disc_with_hole/elasticity.mfront
\ No newline at end of file
diff --git a/Tests/Data/ThermoMechanics/LinearMFront/libBehaviour.so b/Tests/Data/ThermoMechanics/LinearMFront/libBehaviour.so
deleted file mode 100755
index 93632cd232f2d286a7400d749ef622ae8b81f758..0000000000000000000000000000000000000000
Binary files a/Tests/Data/ThermoMechanics/LinearMFront/libBehaviour.so and /dev/null differ
diff --git a/Tests/MaterialLib/MFront.cpp b/Tests/MaterialLib/MFront.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5243dbd36d114ac8d720d3ff81043165c7f1bac4
--- /dev/null
+++ b/Tests/MaterialLib/MFront.cpp
@@ -0,0 +1,149 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/LICENSE.txt
+ */
+
+#ifdef OGS_USE_MFRONT
+
+#include <gtest/gtest.h>
+#include <MGIS/Behaviour/Integrate.hxx>
+
+#include "MaterialLib/SolidModels/MFront/MFront.h"
+#include "ParameterLib/ConstantParameter.h"
+
+using namespace MaterialLib::Solids;
+template <int Dim>
+using KelvinVector = MathLib::KelvinVector::KelvinVectorType<Dim>;
+
+constexpr mgis::behaviour::Hypothesis hypothesis(int dim)
+{
+    switch (dim)
+    {
+        case 2:
+            return mgis::behaviour::Hypothesis::PLANESTRAIN;
+        case 3:
+            return mgis::behaviour::Hypothesis::TRIDIMENSIONAL;
+    }
+    // workaround for simple throw, because of gcc bugs 67371, 66026, and 80061.
+    return (2 < dim && dim > 3) ? throw
+                                // unreachable
+                                : mgis::behaviour::Hypothesis::PLANESTRAIN;
+}
+
+template <int Dim>
+struct StandardElasticityBrickBehaviour
+{
+    static std::unique_ptr<MFront::MFront<Dim>> createConstitutiveRelation()
+    {
+        auto behaviour =
+            mgis::behaviour::load("libOgsMFrontBehaviour.so",
+                                  "StandardElasticityBrick", hypothesis(Dim));
+
+        using P = ParameterLib::ConstantParameter<double>;
+        // Parameters used by mfront model in the order of appearence in the
+        // .mfront file.
+        static P const young_modulus("", 1e11);
+        static P const poisson_ratio("", 1e9);
+        std::vector<ParameterLib::Parameter<double> const*> parameters{
+            &young_modulus, &poisson_ratio};
+
+        auto result = std::make_unique<MFront::MFront<Dim>>(
+            std::move(behaviour), std::move(parameters));
+        return result;
+    }
+};
+template <int Dim>
+struct ElasticBehaviour
+{
+    static std::unique_ptr<MFront::MFront<Dim>> createConstitutiveRelation()
+    {
+        auto behaviour = mgis::behaviour::load("libOgsMFrontBehaviour.so",
+                                               "Elasticity", hypothesis(Dim));
+
+        using P = ParameterLib::ConstantParameter<double>;
+        // Parameters used by mfront model in the order of appearence in the
+        // .mfront file.
+        static P const young_modulus("", 1e11);
+        static P const poisson_ratio("", 1e9);
+        std::vector<ParameterLib::Parameter<double> const*> parameters{
+            &young_modulus, &poisson_ratio};
+
+        auto result = std::make_unique<MFront::MFront<Dim>>(
+            std::move(behaviour), std::move(parameters));
+        return result;
+    }
+};
+
+template <int Dim, typename TestBehaviour>
+struct MaterialLib_SolidModelsMFront : public testing::Test
+{
+    MaterialLib_SolidModelsMFront()
+    {
+        constitutive_relation = TestBehaviour::createConstitutiveRelation();
+    }
+
+    KelvinVector<Dim> const eps_prev = KelvinVector<Dim>::Zero();
+    KelvinVector<Dim> const eps = KelvinVector<Dim>::Zero();
+    KelvinVector<Dim> const sigma_prev = KelvinVector<Dim>::Zero();
+
+    double t = 0;
+    ParameterLib::SpatialPosition x;
+    double dt = 0;
+    double T = 0;
+
+    std::unique_ptr<MechanicsBase<Dim>> constitutive_relation;
+};
+
+template <typename TestBehaviour>
+using MaterialLib_SolidModelsMFront2 =
+    MaterialLib_SolidModelsMFront<2, TestBehaviour>;
+
+template <typename TestBehaviour>
+using MaterialLib_SolidModelsMFront3 =
+    MaterialLib_SolidModelsMFront<3, TestBehaviour>;
+
+template <int Dim>
+using TestBehaviourTypes =
+    ::testing::Types<StandardElasticityBrickBehaviour<Dim>,
+                     ElasticBehaviour<Dim>>;
+
+TYPED_TEST_CASE(MaterialLib_SolidModelsMFront2, TestBehaviourTypes<2>);
+TYPED_TEST_CASE(MaterialLib_SolidModelsMFront3, TestBehaviourTypes<3>);
+
+TYPED_TEST(MaterialLib_SolidModelsMFront2, IntegrateZeroDisplacement)
+{
+    ASSERT_TRUE(this->constitutive_relation != nullptr);
+    auto state = this->constitutive_relation->createMaterialStateVariables();
+
+    auto solution = this->constitutive_relation->integrateStress(
+        this->t, this->x, this->dt, this->eps_prev, this->eps, this->sigma_prev,
+        *state, this->T);
+
+    ASSERT_TRUE(solution != boost::none);
+    state = std::move(std::get<1>(*solution));
+    ASSERT_TRUE(state != nullptr);
+    state.reset(nullptr);
+    ASSERT_TRUE(state == nullptr);
+}
+
+TYPED_TEST(MaterialLib_SolidModelsMFront3, IntegrateZeroDisplacement)
+{
+    ASSERT_TRUE(this->constitutive_relation != nullptr);
+    auto state = this->constitutive_relation->createMaterialStateVariables();
+
+    auto solution = this->constitutive_relation->integrateStress(
+        this->t, this->x, this->dt, this->eps_prev, this->eps, this->sigma_prev,
+        *state, this->T);
+
+    ASSERT_TRUE(solution != boost::none);
+    state = std::move(std::get<1>(*solution));
+    ASSERT_TRUE(state != nullptr);
+    state.reset(nullptr);
+    ASSERT_TRUE(state == nullptr);
+}
+#endif  // OGS_USE_MFRONT
diff --git a/ThirdParty/CMakeLists.txt b/ThirdParty/CMakeLists.txt
index 5226f73019986b5b47736ac9ceebe932928593ef..665147b6f28e4708138897b4c8ef97ff408eab98 100644
--- a/ThirdParty/CMakeLists.txt
+++ b/ThirdParty/CMakeLists.txt
@@ -27,6 +27,15 @@ option(OGS_USE_OPTIONAL_SUBMODULES "Option for enabling optional submodules" OFF
 #iphreeqc
 add_subdirectory(iphreeqc)
 
+# MFrontGenericInterfaceSupport
+if (OGS_USE_MFRONT)
+    set(enable-doxygen-doc OFF CACHE INTERNAL "")
+    set(enable-fortran-bindings OFF CACHE INTERNAL "")
+    set(CMAKE_CXX_STANDARD 11)
+    add_subdirectory(MGIS)
+    set(CMAKE_CXX_STANDARD 17)
+endif()
+
 # VtkFbxConverter
 if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/VtkFbxConverter/CMakeLists.txt
     AND OGS_USE_OPTIONAL_SUBMODULES)
diff --git a/ThirdParty/MGIS b/ThirdParty/MGIS
new file mode 160000
index 0000000000000000000000000000000000000000..bdfcaa13f8ad00a3ac27a264e8f6e1f7c535a5ee
--- /dev/null
+++ b/ThirdParty/MGIS
@@ -0,0 +1 @@
+Subproject commit bdfcaa13f8ad00a3ac27a264e8f6e1f7c535a5ee
diff --git a/scripts/cmake/ConanSetup.cmake b/scripts/cmake/ConanSetup.cmake
index 3235621be37cb108584af904b40c6d11436c51f2..c5cae781a10cd208c6ea13d318e059f9d67b1d9b 100644
--- a/scripts/cmake/ConanSetup.cmake
+++ b/scripts/cmake/ConanSetup.cmake
@@ -61,6 +61,10 @@ if(OGS_USE_CVODE)
     set(CONAN_REQUIRES ${CONAN_REQUIRES} cvode/2.8.2@bilke/stable)
 endif()
 
+if(OGS_USE_MFRONT)
+    set(CONAN_REQUIRES ${CONAN_REQUIRES} tfel/3.2.1@bilke/testing)
+endif()
+
 if(OGS_BUILD_GUI)
     set(CONAN_REQUIRES ${CONAN_REQUIRES}
         shapelib/1.3.0@bilke/stable
diff --git a/scripts/cmake/Find.cmake b/scripts/cmake/Find.cmake
index c24ab8dc8bd24fea60e7f775a285addbfa79b6cd..254955b29e1de49396a15a94976ca9faf7c572a8 100644
--- a/scripts/cmake/Find.cmake
+++ b/scripts/cmake/Find.cmake
@@ -165,7 +165,3 @@ if(OGS_USE_CVODE)
     find_package(CVODE REQUIRED)
     add_definitions(-DCVODE_FOUND)
 endif()
-
-if(OGS_USE_MFRONT)
-    find_package(MGIS REQUIRED)
-endif()
diff --git a/scripts/cmake/FindMGIS.cmake b/scripts/cmake/FindMGIS.cmake
deleted file mode 100644
index 05d0e7d48aa7967c006466a1960f97042f4abc7f..0000000000000000000000000000000000000000
--- a/scripts/cmake/FindMGIS.cmake
+++ /dev/null
@@ -1,16 +0,0 @@
-# Find the MGIS includes and libraries
-#
-#    MGIS_INCLUDE_DIR - Where to find MGIS headers
-#    MGIS_LIBRARY     - The MGIS library to link against.
-#    MGIS_FOUND       - Do not attempt to use if "no" or undefined.
-#
-# MGIS, the MFront Generic Interface Support library
-# See https://github.com/thelfer/MFrontGenericInterfaceSupport
-#     http://tfel.sourceforge.net/generic-behaviours-interface.html
-
-find_path(MGIS_INCLUDE_DIR MGIS)
-
-find_library(MGIS_LIBRARY MFrontGenericInterface)
-
-include(FindPackageHandleStandardArgs)
-FIND_PACKAGE_HANDLE_STANDARD_ARGS(MGIS DEFAULT_MSG MGIS_LIBRARY MGIS_INCLUDE_DIR)
diff --git a/scripts/cmake/SubmoduleSetup.cmake b/scripts/cmake/SubmoduleSetup.cmake
index 25877cf25b8907642c9b147fc4b195641d92b6c4..cb64ecd4dcd7d71280e71d2802f430090da9cf20 100644
--- a/scripts/cmake/SubmoduleSetup.cmake
+++ b/scripts/cmake/SubmoduleSetup.cmake
@@ -26,6 +26,9 @@ endif()
 if(OGS_USE_PYTHON)
     list(APPEND REQUIRED_SUBMODULES ThirdParty/pybind11)
 endif()
+if (OGS_USE_MFRONT)
+    list(APPEND REQUIRED_SUBMODULES ThirdParty/MGIS)
+endif()
 
 # Sync submodules, which is required when a submodule changed its URL
 if(OGS_SYNC_SUBMODULES)
diff --git a/scripts/cmake/test/Test.cmake b/scripts/cmake/test/Test.cmake
index 7d1a81663347f375c2ab503539531cca45d4dbb9..972a2d92be6007bb4c0a651d4a6d45dbce5ded76 100644
--- a/scripts/cmake/test/Test.cmake
+++ b/scripts/cmake/test/Test.cmake
@@ -55,6 +55,23 @@ if(CMAKE_CONFIGURATION_TYPES)
     set(CONFIG_PARAMETER --build-config "$<CONFIGURATION>")
 endif()
 add_custom_target(ctest-cleanup ${CMAKE_COMMAND} -E remove -f Tests/ctest.log)
+
+set(test_dependencies ogs vtkdiff)
+if(OGS_BUILD_UTILS)
+    list(APPEND test_dependencies partmesh MapGeometryToMeshSurface)
+endif()
+if(OGS_USE_MFRONT)
+    list(APPEND test_dependencies
+        MFrontGenericBehaviourInterfaceTest
+        MFrontGenericBehaviourInterfaceTest2
+        BoundsCheckTest
+        ParameterTest
+        IntegrateTest
+        IntegrateTest2
+        IntegrateTest3
+        BehaviourTest)
+endif()
+
 add_custom_target(
     ctest
     COMMAND ${CMAKE_CTEST_COMMAND} -T Test
@@ -63,12 +80,9 @@ add_custom_target(
     --exclude-regex LARGE
     ${CONFIG_PARAMETER} --parallel ${NUM_CTEST_PROCESSORS}
     --timeout 900 # 15 minutes
-    DEPENDS ogs vtkdiff ctest-cleanup
+    DEPENDS ${test_dependencies} ctest-cleanup
     USES_TERMINAL
 )
-if(OGS_BUILD_UTILS)
-    add_dependencies(ctest partmesh MapGeometryToMeshSurface)
-endif()
 
 add_custom_target(
     ctest-serial
@@ -78,14 +92,12 @@ add_custom_target(
     --exclude-regex LARGE
     ${CONFIG_PARAMETER}
     --timeout 900 # 15 minutes
-    DEPENDS ogs vtkdiff ctest-cleanup
+    DEPENDS ${test_dependencies} ctest-cleanup
     USES_TERMINAL
 )
-if(OGS_BUILD_UTILS)
-    add_dependencies(ctest-serial partmesh MapGeometryToMeshSurface)
-endif()
 
 add_custom_target(ctest-large-cleanup ${CMAKE_COMMAND} -E remove -f Tests/ctest-large.log)
+
 add_custom_target(
     ctest-large
     COMMAND ${CMAKE_CTEST_COMMAND} -T Test
@@ -94,12 +106,9 @@ add_custom_target(
     --tests-regex LARGE
     ${CONFIG_PARAMETER} --parallel ${NUM_CTEST_LARGE_PROCESSORS}
     --timeout 3600
-    DEPENDS ogs vtkdiff ctest-large-cleanup
+    DEPENDS ${test_dependencies} ctest-large-cleanup
     USES_TERMINAL
 )
-if(OGS_BUILD_UTILS)
-    add_dependencies(ctest-large partmesh MapGeometryToMeshSurface)
-endif()
 
 add_custom_target(
     ctest-large-serial
@@ -109,18 +118,13 @@ add_custom_target(
     --tests-regex LARGE
     ${CONFIG_PARAMETER}
     --timeout 3600
-    DEPENDS ogs vtkdiff ctest-large-cleanup
+    DEPENDS ${test_dependencies} ctest-large-cleanup
     USES_TERMINAL
 )
-if(OGS_BUILD_UTILS)
-    add_dependencies(ctest-large-serial partmesh MapGeometryToMeshSurface)
-endif()
+
 set_directory_properties(PROPERTIES
     ADDITIONAL_MAKE_CLEAN_FILES ${PROJECT_BINARY_DIR}/Tests/Data
 )
 
-set_target_properties(ctest PROPERTIES FOLDER Testing)
-set_target_properties(ctest-large PROPERTIES FOLDER Testing)
-set_target_properties(ctest-large-serial PROPERTIES FOLDER Testing)
-set_target_properties(ctest-cleanup PROPERTIES FOLDER Testing)
-set_target_properties(ctest-large-cleanup PROPERTIES FOLDER Testing)
+set_target_properties(ctest ctest-large ctest-large-serial ctest-cleanup ctest-large-cleanup
+    PROPERTIES FOLDER Testing)