diff --git a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
index e29ecb9d75b2a5137fc7be8095d8cba0730fa5b6..9bb2c3edfc4a6d3b5c67dac8f773b19a19da7c8b 100644
--- a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
+++ b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
@@ -24,6 +24,13 @@ struct SmallDeformationLocalAssemblerInterface
     : public ProcessLib::LocalAssemblerInterface,
       public NumLib::ExtrapolatableElement
 {
+    virtual std::vector<double> const& getIntPtSigma(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    // TODO remove the component-wise methods.
     virtual std::vector<double> const& getIntPtSigmaXX(
         const double /*t*/,
         GlobalVector const& /*current_solution*/,
@@ -60,6 +67,13 @@ struct SmallDeformationLocalAssemblerInterface
         NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
         std::vector<double>& cache) const = 0;
 
+    virtual std::vector<double> const& getIntPtEpsilon(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    // TODO remove the component-wise methods
     virtual std::vector<double> const& getIntPtEpsilonXX(
         const double /*t*/,
         GlobalVector const& /*current_solution*/,
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index 201b6516c68b6249df90c28d879147fd6da2bced..f1bf21ba09763ec29b96cb3214fe2cbdbab5378c 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
@@ -242,6 +242,46 @@ public:
         return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
     }
 
+    std::vector<double> const& getIntPtSigma(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const override
+    {
+        using KelvinVectorType = typename BMatricesType::KelvinVectorType;
+        auto const kelvin_vector_size =
+            KelvinVectorDimensions<DisplacementDim>::value;
+        auto const num_intpts = _ip_data.size();
+
+        cache.clear();
+        auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+            double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
+            cache, kelvin_vector_size, num_intpts);
+
+        // TODO make a general implementation for converting KelvinVectors
+        // back to symmetric rank-2 tensors.
+        for (unsigned ip = 0; ip < num_intpts; ++ip)
+        {
+            auto const& sigma = _ip_data[ip].sigma;
+
+            for (typename KelvinVectorType::Index component = 0;
+                 component < kelvin_vector_size && component < 3;
+                 ++component)
+            {  // xx, yy, zz components
+                cache_mat(component, ip) = sigma[component];
+            }
+            for (typename KelvinVectorType::Index component = 3;
+                 component < kelvin_vector_size;
+                 ++component)
+            {  // mixed xy, yz, xz components
+                cache_mat(component, ip) = sigma[component] / std::sqrt(2);
+            }
+        }
+
+        return cache;
+    }
+
+    // TODO remove the component-wise methods.
     std::vector<double> const& getIntPtSigmaXX(
         const double /*t*/,
         GlobalVector const& /*current_solution*/,
@@ -298,6 +338,46 @@ public:
         return getIntPtSigma(cache, 5);
     }
 
+    virtual std::vector<double> const& getIntPtEpsilon(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const override
+    {
+        using KelvinVectorType = typename BMatricesType::KelvinVectorType;
+        auto const kelvin_vector_size =
+            KelvinVectorDimensions<DisplacementDim>::value;
+        auto const num_intpts = _ip_data.size();
+
+        cache.clear();
+        auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+            double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
+            cache, kelvin_vector_size, num_intpts);
+
+        // TODO make a general implementation for converting KelvinVectors
+        // back to symmetric rank-2 tensors.
+        for (unsigned ip = 0; ip < num_intpts; ++ip)
+        {
+            auto const& eps = _ip_data[ip].eps;
+
+            for (typename KelvinVectorType::Index component = 0;
+                 component < kelvin_vector_size && component < 3;
+                 ++component)
+            {  // xx, yy, zz components
+                cache_mat(component, ip) = eps[component];
+            }
+            for (typename KelvinVectorType::Index component = 3;
+                 component < kelvin_vector_size;
+                 ++component)
+            {  // mixed xy, yz, xz components
+                cache_mat(component, ip) = eps[component] / std::sqrt(2);
+            }
+        }
+
+        return cache;
+    }
+
+    // TODO remove the component-wise methods
     std::vector<double> const& getIntPtEpsilonXX(
         const double /*t*/,
         GlobalVector const& /*current_solution*/,
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h b/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h
index 633f1e00c49c8c4c0c49f9cf36c84526a35f5d10..dd9363a9f20e5ed6d2de6b53c7d9d6e34501286e 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h
@@ -69,6 +69,14 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
             NumLib::ComponentOrder::BY_LOCATION);
     _nodal_forces->resize(DisplacementDim * mesh.getNumberOfNodes());
 
+    Base::_secondary_variables.addSecondaryVariable(
+        "sigma",
+        makeExtrapolator(
+            ProcessLib::KelvinVectorType<DisplacementDim>::RowsAtCompileTime,
+            getExtrapolator(), _local_assemblers,
+            &LocalAssemblerInterface::getIntPtSigma));
+
+    // TODO remove the component-wise methods
     Base::_secondary_variables.addSecondaryVariable(
         "sigma_xx",
         makeExtrapolator(1, getExtrapolator(), _local_assemblers,
@@ -102,6 +110,14 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
                              &LocalAssemblerInterface::getIntPtSigmaYZ));
     }
 
+    Base::_secondary_variables.addSecondaryVariable(
+        "epsilon",
+        makeExtrapolator(
+            ProcessLib::KelvinVectorType<DisplacementDim>::RowsAtCompileTime,
+            getExtrapolator(), _local_assemblers,
+            &LocalAssemblerInterface::getIntPtEpsilon));
+
+    // TODO remove the component-wise methods
     Base::_secondary_variables.addSecondaryVariable(
         "epsilon_xx",
         makeExtrapolator(1, getExtrapolator(), _local_assemblers,
@@ -134,6 +150,7 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
                              &LocalAssemblerInterface::getIntPtEpsilonXZ));
     }
 
+    // enable output of internal variables defined by material models
     auto const internal_variables =
         _process_data.material->getInternalVariables();
     for (auto const& internal_variable : internal_variables)