diff --git a/opm/models/io/vtktpsamodule.hpp b/opm/models/io/vtktpsamodule.hpp
index 6602c3bb021..10689f0b720 100644
--- a/opm/models/io/vtktpsamodule.hpp
+++ b/opm/models/io/vtktpsamodule.hpp
@@ -28,6 +28,8 @@
 #include <dune/common/dynmatrix.hh>
 #include <dune/common/fvector.hh>
 
+#include <opm/common/utility/SymmTensor.hpp>
+#include <opm/common/utility/VoigtArray.hpp>
 #include <opm/material/densead/Math.hpp>
 
 #include <opm/models/io/baseoutputmodule.hh>
@@ -64,7 +66,7 @@ class VtkTpsaModule : public BaseOutputModule<TypeTag>
     using VectorBuffer = typename ParentType::VectorBuffer;
     using TensorBuffer = typename ParentType::TensorBuffer;
 
-    using SymTensor = Dune::FieldVector<Scalar, 6>;
+    using SymTensor = SymmTensor<Scalar>;
     using Tensor = Dune::DynamicMatrix<Scalar>;
 
 public:
@@ -279,14 +281,15 @@ class VtkTpsaModule : public BaseOutputModule<TypeTag>
     static void setTensorFromVoigt_(Tensor& tensor, const SymTensor& symTensor)
     {
         // Diagonal terms
+        constexpr auto& ind = SymTensor::diag_indices;
         for (std::size_t i = 0; i < 3; ++i) {
-            tensor[i][i] = symTensor[i];
+            tensor[i][i] = symTensor[ind[i]];
         }
 
         // Off-diagonal terms
-        tensor[0][1] = symTensor[5];
-        tensor[0][2] = symTensor[4];
-        tensor[1][2] = symTensor[3];
+        tensor[0][1] = symTensor[VoigtIndex::XY];
+        tensor[0][2] = symTensor[VoigtIndex::XZ];
+        tensor[1][2] = symTensor[VoigtIndex::YZ];
         for (std::size_t i = 0; i < 3; ++i) {
             for (std::size_t j = 0; j < 3; ++j) {
                 if (i > j) {
diff --git a/opm/models/tpsa/tpsamodel.hpp b/opm/models/tpsa/tpsamodel.hpp
index 9cb89562019..14b7ce1f4ed 100644
--- a/opm/models/tpsa/tpsamodel.hpp
+++ b/opm/models/tpsa/tpsamodel.hpp
@@ -32,6 +32,8 @@
 
 #include <opm/grid/utility/ElementChunks.hpp>
 
+#include <opm/common/utility/SymmTensor.hpp>
+#include <opm/common/utility/VoigtArray.hpp>
 #include <opm/material/common/MathToolbox.hpp>
 
 #include <opm/models/parallel/threadmanager.hpp>
@@ -77,7 +79,7 @@ class TpsaModel
     using MaterialState = MaterialStateTPSA<Evaluation>;
 
     using DimVector = Dune::FieldVector<Scalar, dimWorld>;
-    using SymTensor = Dune::FieldVector<Scalar, 6>;
+    using SymTensor = SymmTensor<Scalar>;
     using PotForceVector = Dune::BlockVector<Scalar>;
 
 public:
@@ -491,7 +493,7 @@ class TpsaModel
         // the total stress
         SymTensor linStressTensor = stress(globalIdx, false);
         const auto potForce = mechPotentialForce(globalIdx);
-        for (unsigned dirIdx = 0; dirIdx < 3; ++dirIdx) {
+        for (const auto& dirIdx : SymTensor::diag_indices) {
             linStressTensor[dirIdx] -= potForce;
         }
         return -1.0 * linStressTensor;
@@ -591,12 +593,12 @@ class TpsaModel
 
             // Reconstructed stress tensor at cell center
             const auto& stress = lsq.x();
-            stressOutput[0] = stress[0]; // XX
-            stressOutput[1] = stress[1]; // YY
-            stressOutput[2] = stress[2]; // ZZ
-            stressOutput[3] = stress[5]; // YZ
-            stressOutput[4] = stress[4]; // XZ
-            stressOutput[5] = stress[3]; // XY
+            stressOutput[VoigtIndex::XX] = stress[0]; // XX
+            stressOutput[VoigtIndex::YY] = stress[1]; // YY
+            stressOutput[VoigtIndex::ZZ] = stress[2]; // ZZ
+            stressOutput[VoigtIndex::YZ] = stress[5]; // YZ
+            stressOutput[VoigtIndex::XZ] = stress[4]; // XZ
+            stressOutput[VoigtIndex::XY] = stress[3]; // XY
         }
         return stressOutput;
     }
@@ -612,10 +614,10 @@ class TpsaModel
     {
         // Deviatoric stress
         auto stressDev = this->linstress(globalIdx);
-        Scalar traceStress = 1.0 / 3.0 * (stressDev[0] + stressDev[1] + stressDev[2]);
-        stressDev[0] -= traceStress;
-        stressDev[1] -= traceStress;
-        stressDev[2] -= traceStress;
+        Scalar traceStress = 1.0 / 3.0 * stressDev.trace();
+        stressDev[VoigtIndex::XX] -= traceStress;
+        stressDev[VoigtIndex::YY] -= traceStress;
+        stressDev[VoigtIndex::ZZ] -= traceStress;
 
         // Deviatoric strain
         auto& problem = simulator_.problem();
@@ -627,9 +629,9 @@ class TpsaModel
         Scalar strainVolTerm = traceStress / (3.0 * lameParam + 2 * sMod);
 
         SymTensor strainVol;
-        strainVol[0] = strainVolTerm;
-        strainVol[1] = strainVolTerm;
-        strainVol[2] = strainVolTerm;
+        strainVol[VoigtIndex::XX] = strainVolTerm;
+        strainVol[VoigtIndex::YY] = strainVolTerm;
+        strainVol[VoigtIndex::ZZ] = strainVolTerm;
 
         // Total
         SymTensor strainOutput = strainDev + strainVol;
diff --git a/opm/simulators/flow/MechContainer.cpp b/opm/simulators/flow/MechContainer.cpp
index ea5a3e21b9f..7ed9de4a629 100644
--- a/opm/simulators/flow/MechContainer.cpp
+++ b/opm/simulators/flow/MechContainer.cpp
@@ -78,9 +78,9 @@ allocate(const std::size_t bufferSize,
 template<class Scalar>
 void MechContainer<Scalar>::
 assignDelStress(const unsigned globalDofIdx,
-                const Dune::FieldVector<Scalar,6>& delStress)
+                const SymmTensor<Scalar>& delStress)
 {
-    this->delstress_.assign(globalDofIdx, VoigtContainer<Scalar>(delStress));
+    this->delstress_.assign(globalDofIdx, delStress);
 }
 
 template<class Scalar>
@@ -96,17 +96,17 @@ assignDisplacement(const unsigned globalDofIdx,
 template<class Scalar>
 void MechContainer<Scalar>::
 assignFracStress(const unsigned globalDofIdx,
-                 const Dune::FieldVector<Scalar,6>& fracStress)
+                 const SymmTensor<Scalar>& fracStress)
 {
-    this->fracstress_.assign(globalDofIdx, VoigtContainer<Scalar>(fracStress));
+    this->fracstress_.assign(globalDofIdx, fracStress);
 }
 
 template<class Scalar>
 void MechContainer<Scalar>::
 assignLinStress(const unsigned globalDofIdx,
-                const Dune::FieldVector<Scalar,6>& linStress)
+                const SymmTensor<Scalar>& linStress)
 {
-    this->linstress_.assign(globalDofIdx, VoigtContainer<Scalar>(linStress));
+    this->linstress_.assign(globalDofIdx, linStress);
 }
 
 template<class Scalar>
@@ -124,17 +124,17 @@ assignPotentialForces(const unsigned globalDofIdx,
 template<class Scalar>
 void MechContainer<Scalar>::
 assignStrain(const unsigned globalDofIdx,
-             const Dune::FieldVector<Scalar,6>& strain)
+             const SymmTensor<Scalar>& strain)
 {
-    this->strain_.assign(globalDofIdx, VoigtContainer<Scalar>(strain));
+    this->strain_.assign(globalDofIdx, strain);
 }
 
 template<class Scalar>
 void MechContainer<Scalar>::
 assignStress(const unsigned globalDofIdx,
-             const Dune::FieldVector<Scalar,6>& stress)
+             const SymmTensor<Scalar>& stress)
 {
-    this->stress_.assign(globalDofIdx, VoigtContainer<Scalar>(stress));
+    this->stress_.assign(globalDofIdx, stress);
 }
 
 template<class Scalar>
diff --git a/opm/simulators/flow/MechContainer.hpp b/opm/simulators/flow/MechContainer.hpp
index d03e0b7e3a7..e1fcc9ff86f 100644
--- a/opm/simulators/flow/MechContainer.hpp
+++ b/opm/simulators/flow/MechContainer.hpp
@@ -28,6 +28,7 @@
 
 #include <dune/common/fvector.hh>
 
+#include <opm/common/utility/SymmTensor.hpp>
 #include <opm/common/utility/VoigtArray.hpp>
 
 #include <array>
@@ -53,7 +54,7 @@ class MechContainer
                             const Dune::FieldVector<Scalar,3>& disp);
 
     void assignDelStress(const unsigned globalDofIdx,
-                         const Dune::FieldVector<Scalar,6>& delStress);
+                         const SymmTensor<Scalar>& delStress);
 
     void assignPotentialForces(const unsigned globalDofIdx,
                                const Scalar force,
@@ -61,16 +62,16 @@ class MechContainer
                                const Scalar tempForce);
 
     void assignFracStress(const unsigned globalDofIdx,
-                          const Dune::FieldVector<Scalar,6>& fracStress);
+                          const SymmTensor<Scalar>& fracStress);
 
     void assignLinStress(const unsigned globalDofIdx,
-                         const Dune::FieldVector<Scalar,6>& linStress);
+                         const SymmTensor<Scalar>& linStress);
 
     void assignStrain(const unsigned globalDofIdx,
-                      const Dune::FieldVector<Scalar,6>& strain);
+                      const SymmTensor<Scalar>& strain);
 
     void assignStress(const unsigned globalDofIdx,
-                      const Dune::FieldVector<Scalar,6>& stress);
+                      const SymmTensor<Scalar>& stress);
 
     void outputRestart(data::Solution& sol);