diff --git a/MaterialLib/FractureModels/MohrCoulomb.cpp b/MaterialLib/FractureModels/MohrCoulomb.cpp
index 8a403be9230eb131fd7d39da4a79bc521299c8f0..dd68dfc9d618997450fe4d49362b3f58801e8756 100644
--- a/MaterialLib/FractureModels/MohrCoulomb.cpp
+++ b/MaterialLib/FractureModels/MohrCoulomb.cpp
@@ -66,8 +66,7 @@ void MohrCoulomb<DisplacementDim>::computeConstitutiveRelation(
     MaterialPropertyValues const mat(_mp, t, x);
     Eigen::VectorXd const dw = w - w_prev;
 
-    Eigen::MatrixXd Ke(2,2);
-    Ke.setZero();
+    Eigen::MatrixXd Ke = Eigen::MatrixXd::Zero(2, 2);
     Ke(0,0) = mat.Ks;
     Ke(1,1) = mat.Kn;
 
diff --git a/MaterialLib/SolidModels/Lubby2-impl.h b/MaterialLib/SolidModels/Lubby2-impl.h
index 8ece8e1dd01c3a3d0138e0f1b16fc3d5fe0474d3..11bb4d6f9a747028a8db10ac12a37bafa019ca52 100644
--- a/MaterialLib/SolidModels/Lubby2-impl.h
+++ b/MaterialLib/SolidModels/Lubby2-impl.h
@@ -29,8 +29,9 @@ calculatedGdEBurgers()
 {
     Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
                   Lubby2<DisplacementDim>::KelvinVectorSize>
-        dGdE;
-    dGdE.setZero();
+        dGdE =
+            Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
+                          Lubby2<DisplacementDim>::KelvinVectorSize>::Zero();
     dGdE.template topLeftCorner<Lubby2<DisplacementDim>::KelvinVectorSize,
                                 Lubby2<DisplacementDim>::KelvinVectorSize>()
         .diagonal()
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index 31797fd333dc8989a0fb9b436ae8cbf5cb4eacb8..6c92ed322e6d93f07d039aff40f81536af0c34e2 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -291,16 +291,19 @@ public:
                 displacement_size + pressure_size>>(
             local_rhs_data, displacement_size + pressure_size);
 
-        typename ShapeMatricesTypePressure::NodalMatrixType laplace_p;
-        laplace_p.setZero(pressure_size, pressure_size);
+        typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
+            ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+                                                             pressure_size);
 
-        typename ShapeMatricesTypePressure::NodalMatrixType storage_p;
-        storage_p.setZero(pressure_size, pressure_size);
+        typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
+            ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+                                                             pressure_size);
 
         typename ShapeMatricesTypeDisplacement::template MatrixType<
             displacement_size, pressure_size>
-            Kup;
-        Kup.setZero(displacement_size, pressure_size);
+            Kup = ShapeMatricesTypeDisplacement::template MatrixType<
+                displacement_size, pressure_size>::Zero(displacement_size,
+                                                        pressure_size);
 
         double const& dt = _process_data.dt;
 
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
index b7a6c4c60ad02f774a957dd919ddcff425eea137..81cca9a982186f4151608cb18c619361813fe2a7 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
@@ -163,23 +163,24 @@ assembleBlockMatricesWithJacobian(
     // in a displacement vector
     auto const index_normal = GlobalDim - 1;
 
-    typename ShapeMatricesTypePressure::NodalMatrixType laplace_p;
-    laplace_p.setZero(pressure_size, pressure_size);
+    typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
+        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+                                                         pressure_size);
 
-    typename ShapeMatricesTypePressure::NodalMatrixType storage_p;
-    storage_p.setZero(pressure_size, pressure_size);
+    typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
+        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+                                                         pressure_size);
 
     typename ShapeMatricesTypeDisplacement::template MatrixType<
         displacement_size, pressure_size>
-        Kgp;
-    Kgp.setZero(displacement_size, pressure_size);
+        Kgp = ShapeMatricesTypeDisplacement::template MatrixType<
+            displacement_size, pressure_size>::Zero(displacement_size,
+                                                    pressure_size);
 
     using GlobalDimMatrix = Eigen::Matrix<double, GlobalDim, GlobalDim>;
     using GlobalDimVector = Eigen::Matrix<double, GlobalDim, 1>;
-    GlobalDimMatrix local_k_tensor;
-    local_k_tensor.setZero();
-    GlobalDimMatrix local_dk_db_tensor;
-    local_dk_db_tensor.setZero();
+    GlobalDimMatrix local_k_tensor = GlobalDimMatrix::Zero();
+    GlobalDimMatrix local_dk_db_tensor = GlobalDimMatrix::Zero();
 
     auto const& gravity_vec = _process_data.specific_body_force;
 
@@ -309,8 +310,7 @@ computeSecondaryVariableConcreteWithVector(
 
     using GlobalDimMatrix = Eigen::Matrix<double, GlobalDim, GlobalDim>;
     using GlobalDimVector = Eigen::Matrix<double, GlobalDim, 1>;
-    GlobalDimMatrix local_k_tensor;
-    local_k_tensor.setZero();
+    GlobalDimMatrix local_k_tensor = GlobalDimMatrix::Zero();
 
     auto const& gravity_vec = _process_data.specific_body_force;
 
@@ -364,10 +364,8 @@ computeSecondaryVariableConcreteWithVector(
 
     double ele_b = 0;
     double ele_k = 0;
-    Eigen::Vector2d ele_sigma_eff;
-    ele_sigma_eff.setZero();
-    Eigen::Vector2d ele_w;
-    ele_w.setZero();
+    Eigen::Vector2d ele_sigma_eff = Eigen::Vector2d::Zero();
+    Eigen::Vector2d ele_w = Eigen::Vector2d::Zero();
     double ele_Fs = - std::numeric_limits<double>::max();
     for (auto const& ip : _ip_data)
     {
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
index e2cdc3eca41e577fa220cafc993eea637ccd55b7..196c2403f3e5c513000774df3af4236eb1efdfb1 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
@@ -163,16 +163,19 @@ assembleBlockMatricesWithJacobian(
 {
     assert (this->_element.getDimension() == GlobalDim);
 
-    typename ShapeMatricesTypePressure::NodalMatrixType laplace_p;
-    laplace_p.setZero(pressure_size, pressure_size);
+    typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
+        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+                                                         pressure_size);
 
-    typename ShapeMatricesTypePressure::NodalMatrixType storage_p;
-    storage_p.setZero(pressure_size, pressure_size);
+    typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
+        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+                                                         pressure_size);
 
     typename ShapeMatricesTypeDisplacement::template MatrixType<
         displacement_size, pressure_size>
-        Kup;
-    Kup.setZero(displacement_size, pressure_size);
+        Kup = ShapeMatricesTypeDisplacement::template MatrixType<
+            displacement_size, pressure_size>::Zero(displacement_size,
+                                                    pressure_size);
 
     double const& dt = _process_data.dt;
     auto const& gravity_vec = _process_data.specific_body_force;
@@ -363,12 +366,9 @@ computeSecondaryVariableConcreteWithBlockVectors(
         q.noalias() = - k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
     }
 
-    Eigen::Vector3d ele_stress;
-    ele_stress.setZero();
-    Eigen::Vector3d ele_strain;
-    ele_strain.setZero();
-    Eigen::Vector3d ele_velocity;
-    ele_velocity.setZero();
+    Eigen::Vector3d ele_stress = Eigen::Vector3d::Zero();
+    Eigen::Vector3d ele_strain = Eigen::Vector3d::Zero();
+    Eigen::Vector3d ele_velocity = Eigen::Vector3d::Zero();
 
     for (auto const& ip_data : _ip_data)
     {
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h
index 00549b3d140bf6f38c2adf9b64571182c57ed227..8fbf7a92b69e2f9c6df5b6351120b16db8076184 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h
@@ -100,8 +100,8 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
     auto Met = local_M.template block<temperature_size, temperature_size>(
         temperature_matrix_index, temperature_matrix_index);
 
-    NodalMatrixType laplace_operator;
-    laplace_operator.setZero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
+    NodalMatrixType laplace_operator =
+        NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
 
     auto Kgp =
         local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
index ee734e2f5e2f0a41acb0643c0e61f95d6c6b90c1..1ec88af26ab1f69a30efdd554754847b4ea41616 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
@@ -69,8 +69,8 @@ void TwoPhaseFlowWithPPLocalAssembler<
     auto Mlpc = local_M.template block<cap_pressure_size, cap_pressure_size>(
         cap_pressure_matrix_index, cap_pressure_matrix_index);
 
-    NodalMatrixType laplace_operator;
-    laplace_operator.setZero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
+    NodalMatrixType laplace_operator =
+        NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
 
     auto Kgp =
         local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler-impl.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler-impl.h
index 7211da2309c26958eb80beae0690542040541750..14d0447ae3d31441891c4179af1af04afedaad44 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler-impl.h
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler-impl.h
@@ -51,8 +51,8 @@ void TwoPhaseFlowWithPrhoLocalAssembler<
     auto Mlx = local_M.template block<cap_pressure_size, cap_pressure_size>(
         cap_pressure_matrix_index, cap_pressure_matrix_index);
 
-    NodalMatrixType laplace_operator;
-    laplace_operator.setZero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
+    NodalMatrixType laplace_operator =
+        NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
 
     auto Kgp =
         local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
diff --git a/Tests/NumLib/TestFe.cpp b/Tests/NumLib/TestFe.cpp
index 13a71fefa27771f9c5b45044b2afdf73a8d0def2..e2b8cb390a38234bf32285da57d0ff1b2ec9eb47 100644
--- a/Tests/NumLib/TestFe.cpp
+++ b/Tests/NumLib/TestFe.cpp
@@ -191,8 +191,7 @@ TYPED_TEST(NumLibFemIsoTest, CheckMassMatrix)
     FeType fe(*this->mesh_element);
 
     // evaluate a mass matrix M = int{ N^T D N }dA_e
-    NodalMatrix M(this->e_nnodes, this->e_nnodes);
-    M.setZero();
+    NodalMatrix M = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
     ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
     for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
         shape.setZero();
@@ -216,8 +215,7 @@ TYPED_TEST(NumLibFemIsoTest, CheckLaplaceMatrix)
     FeType fe(*this->mesh_element);
 
     // evaluate a Laplace matrix K = int{ dNdx^T D dNdx }dA_e
-    NodalMatrix K(this->e_nnodes, this->e_nnodes);
-    K.setZero();
+    NodalMatrix K = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
     ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
     for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
         shape.setZero();
@@ -241,10 +239,8 @@ TYPED_TEST(NumLibFemIsoTest, CheckMassLaplaceMatrices)
     FeType fe(*this->mesh_element);
 
     // evaluate both mass and laplace matrices at once
-    NodalMatrix M(this->e_nnodes, this->e_nnodes);
-    M.setZero();
-    NodalMatrix K(this->e_nnodes, this->e_nnodes);
-    K.setZero();
+    NodalMatrix M = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
+    NodalMatrix K = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
     ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
     for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
         shape.setZero();