diff --git a/ProcessLib/PhaseField/PhaseFieldProcess-impl.h b/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
index a6d92bfab2c1c1796c93bf34603f3e2e6369e94d..4a8e70af6eb2f3ad68149c41518a809c70358258 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
@@ -40,6 +40,8 @@ PhaseFieldProcess<DisplacementDim>::PhaseFieldProcess(
               use_monolithic_scheme),
       _process_data(std::move(process_data))
 {
+    _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
+        mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
 }
 
 template <int DisplacementDim>
@@ -158,6 +160,8 @@ void PhaseFieldProcess<DisplacementDim>::initializeConcreteProcess(
         mesh.getElements(), dof_table, _local_assemblers,
         mesh.isAxiallySymmetric(), integration_order, _process_data);
 
+    _nodal_forces->resize(DisplacementDim * mesh.getNumberOfNodes());
+
     Base::_secondary_variables.addSecondaryVariable(
         "sigma_xx",
         makeExtrapolator(1, getExtrapolator(), _local_assemblers,
@@ -327,6 +331,13 @@ void PhaseFieldProcess<DisplacementDim>::assembleWithJacobianConcreteProcess(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, dof_tables, t, x, xdot, dxdot_dx, dx_dx, M, K, b,
         Jac, _coupled_solutions);
+
+    if (_use_monolithic_scheme || (_coupled_solutions->process_id == 0))
+    {
+        b.copyValues(*_nodal_forces);
+        std::transform(_nodal_forces->begin(), _nodal_forces->end(),
+                       _nodal_forces->begin(), [](double val) { return -val; });
+    }
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess.h b/ProcessLib/PhaseField/PhaseFieldProcess.h
index 009cb18506cbb8b51d94625b6e1aed3aa00c25d3..9860b946153879cc31fd2a4a7b147b5505ee1c8e 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess.h
+++ b/ProcessLib/PhaseField/PhaseFieldProcess.h
@@ -114,6 +114,8 @@ private:
     std::unique_ptr<NumLib::LocalToGlobalIndexMap>
         _local_to_global_index_map_single_component;
 
+    MeshLib::PropertyVector<double>* _nodal_forces = nullptr;
+
     /// Sparsity pattern for the phase field equation, and it is initialized
     ///  only if the staggered scheme is used.
     GlobalSparsityPattern _sparsity_pattern_with_single_component;