diff --git a/ProcessLib/SmallDeformation/ConstitutiveRelations/StressData.h b/ProcessLib/SmallDeformation/ConstitutiveRelations/StressData.h
index 751db9aca34c4541338510fcf8169deb2faecf33..90cf7d5df202042707bddf6439c368dda3045ea8 100644
--- a/ProcessLib/SmallDeformation/ConstitutiveRelations/StressData.h
+++ b/ProcessLib/SmallDeformation/ConstitutiveRelations/StressData.h
@@ -17,10 +17,7 @@ namespace ProcessLib::SmallDeformation
 template <int DisplacementDim>
 struct StressData
 {
-    // Stress is stateful for some constitutive settings and therefore must be
-    // initialized to something valid, e.g., zero.
-    // TODO find a better solution for that.
-    KelvinVector<DisplacementDim> sigma = KVzero<DisplacementDim>();
+    KelvinVector<DisplacementDim> sigma = KVnan<DisplacementDim>();
 
     static auto reflect()
     {
diff --git a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
index 47283639c14b35a85d3cd666f729f50070a4a315..a90d95c6c5dab7f93f0d37038134405883cf70d4 100644
--- a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
+++ b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
@@ -52,16 +52,22 @@ struct SmallDeformationLocalAssemblerInterface
         unsigned const n_integration_points =
             integration_method_.getNumberOfPoints();
 
+        current_states_.resize(n_integration_points);
+        prev_states_.resize(n_integration_points);
+        output_data_.resize(n_integration_points);
+
         material_states_.reserve(n_integration_points);
         for (unsigned ip = 0; ip < n_integration_points; ++ip)
         {
             material_states_.emplace_back(
                 solid_material_.createMaterialStateVariables());
-        }
 
-        current_states_.resize(n_integration_points);
-        prev_states_.resize(n_integration_points);
-        output_data_.resize(n_integration_points);
+            // Set initial stress field to zero. Might be overwritten by
+            // integration point data or initial stress.
+            this->current_states_[ip].stress_data.sigma.noalias() =
+                MathLib::KelvinVector::KelvinVectorType<
+                    DisplacementDim>::Zero();
+        }
     }
     /// Returns number of read integration points.
     std::size_t setIPDataInitialConditions(std::string_view name,