diff --git a/src/scheme/ScalarDiamondScheme.cpp b/src/scheme/ScalarDiamondScheme.cpp
index 32eafcc7575d8c89b887e326fc2a57f4fff38574..07561e6f5abfee4461cd72d2e634bce94f7f12e8 100644
--- a/src/scheme/ScalarDiamondScheme.cpp
+++ b/src/scheme/ScalarDiamondScheme.cpp
@@ -281,13 +281,10 @@ class ScalarDiamondSchemeHandler::ScalarDiamondScheme : public ScalarDiamondSche
                       const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f,
                       const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
   {
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != dual_mu->mesh()) {
-      throw NormalError("diffusion coefficient is not defined on the dual mesh!");
-    }
-
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != dual_mub->mesh()) {
-      throw NormalError("boundary diffusion coefficient is not defined on the dual mesh!");
-    }
+    Assert(DualMeshManager::instance().getDiamondDualMesh(*mesh) == dual_mu->mesh(),
+           "diffusion coefficient is not defined on the dual mesh!");
+    Assert(DualMeshManager::instance().getDiamondDualMesh(*mesh) == dual_mub->mesh(),
+           "boundary diffusion coefficient is not defined on the dual mesh!");
 
     using BoundaryCondition = std::variant<DirichletBoundaryCondition, FourierBoundaryCondition,
                                            NeumannBoundaryCondition, SymmetryBoundaryCondition>;
@@ -781,6 +778,13 @@ ScalarDiamondSchemeHandler::ScalarDiamondSchemeHandler(
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
 {
   const std::shared_ptr i_mesh = getCommonMesh({alpha, f});
+  if (not i_mesh) {
+    throw NormalError("primal discrete functions are not defined on the same mesh");
+  }
+  const std::shared_ptr i_dual_mesh = getCommonMesh({dual_mub, dual_mu});
+  if (not i_dual_mesh) {
+    throw NormalError("dual discrete functions are not defined on the same mesh");
+  }
   checkDiscretizationType({alpha, dual_mub, dual_mu, f}, DiscreteFunctionType::P0);
 
   switch (i_mesh->dimension()) {
@@ -788,9 +792,14 @@ ScalarDiamondSchemeHandler::ScalarDiamondSchemeHandler(
     using MeshType                   = Mesh<Connectivity<1>>;
     using DiscreteScalarFunctionType = DiscreteFunctionP0<1, double>;
 
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
+      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
+    }
+
     m_scheme =
-      std::make_unique<ScalarDiamondScheme<1>>(std::dynamic_pointer_cast<const MeshType>(i_mesh),
-                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+      std::make_unique<ScalarDiamondScheme<1>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
                                                std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
                                                std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
                                                std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
@@ -801,9 +810,14 @@ ScalarDiamondSchemeHandler::ScalarDiamondSchemeHandler(
     using MeshType                   = Mesh<Connectivity<2>>;
     using DiscreteScalarFunctionType = DiscreteFunctionP0<2, double>;
 
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
+      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
+    }
+
     m_scheme =
-      std::make_unique<ScalarDiamondScheme<2>>(std::dynamic_pointer_cast<const MeshType>(i_mesh),
-                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+      std::make_unique<ScalarDiamondScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
                                                std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
                                                std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
                                                std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
@@ -814,9 +828,14 @@ ScalarDiamondSchemeHandler::ScalarDiamondSchemeHandler(
     using MeshType                   = Mesh<Connectivity<3>>;
     using DiscreteScalarFunctionType = DiscreteFunctionP0<3, double>;
 
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
+      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
+    }
+
     m_scheme =
-      std::make_unique<ScalarDiamondScheme<3>>(std::dynamic_pointer_cast<const MeshType>(i_mesh),
-                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+      std::make_unique<ScalarDiamondScheme<3>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
                                                std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
                                                std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
                                                std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
diff --git a/src/scheme/VectorDiamondScheme.cpp b/src/scheme/VectorDiamondScheme.cpp
index 32cfbdab5904f2f8ceaa025f2539490f6e52a03c..368713589c323baa59201ea312348837293eba5d 100644
--- a/src/scheme/VectorDiamondScheme.cpp
+++ b/src/scheme/VectorDiamondScheme.cpp
@@ -289,9 +289,8 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
     Assert(dual_lambda->mesh() == dual_mu->mesh());
     Assert(dual_lambdab->mesh() == dual_mu->mesh());
     Assert(dual_mub->mesh() == dual_mu->mesh());
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != dual_mu->mesh()) {
-      throw NormalError("diffusion coefficient is not defined on the dual mesh!");
-    }
+    Assert(DualMeshManager::instance().getDiamondDualMesh(*mesh) == dual_mu->mesh(),
+           "diffusion coefficient is not defined on the dual mesh!");
 
     using MeshDataType = MeshData<Dimension>;
 
@@ -1244,9 +1243,8 @@ class EnergyComputerHandler::EnergyComputer : public EnergyComputerHandler::IEne
     Assert(dual_lambdab->mesh() == dual_U->mesh());
     Assert(U->mesh() == source->mesh());
     Assert(dual_lambdab->mesh() == dual_mub->mesh());
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != dual_mub->mesh()) {
-      throw NormalError("diffusion coefficient is not defined on the dual mesh!");
-    }
+    Assert(DualMeshManager::instance().getDiamondDualMesh(*mesh) == dual_mub->mesh(),
+           "diffusion coefficient is not defined on the dual mesh!");
 
     using MeshDataType = MeshData<Dimension>;
 
@@ -1865,8 +1863,14 @@ VectorDiamondSchemeHandler::VectorDiamondSchemeHandler(
   const std::shared_ptr<const IDiscreteFunction>& f,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
 {
-  const std::shared_ptr i_mesh      = getCommonMesh({alpha, f});
-  const std::shared_ptr i_dual_mesh = getCommonMesh({dual_lambda, dual_mu});
+  const std::shared_ptr i_mesh = getCommonMesh({alpha, f});
+  if (not i_mesh) {
+    throw NormalError("primal discrete functions are not defined on the same mesh");
+  }
+  const std::shared_ptr i_dual_mesh = getCommonMesh({dual_lambda, dual_lambdab, dual_mu, dual_mub});
+  if (not i_dual_mesh) {
+    throw NormalError("dual discrete functions are not defined on the same mesh");
+  }
   checkDiscretizationType({alpha, dual_lambdab, dual_mub, dual_lambda, dual_mu, f}, DiscreteFunctionType::P0);
 
   switch (i_mesh->dimension()) {
@@ -1877,10 +1881,6 @@ VectorDiamondSchemeHandler::VectorDiamondSchemeHandler(
 
     std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
 
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != dual_mu->mesh()) {
-      throw NormalError("mu_dual is not defined on the dual mesh");
-    }
-
     m_scheme =
       std::make_unique<VectorDiamondScheme<1>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
                                                std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(
@@ -1899,8 +1899,8 @@ VectorDiamondSchemeHandler::VectorDiamondSchemeHandler(
 
     std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
 
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != dual_mu->mesh()) {
-      throw NormalError("mu_dual is not defined on the dual mesh");
+    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
+      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
     }
 
     m_scheme =
@@ -1921,8 +1921,8 @@ VectorDiamondSchemeHandler::VectorDiamondSchemeHandler(
 
     std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
 
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != dual_mu->mesh()) {
-      throw NormalError("mu_dual is not defined on the dual mesh");
+    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
+      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
     }
 
     m_scheme =
@@ -1952,9 +1952,15 @@ EnergyComputerHandler::EnergyComputerHandler(
   const std::shared_ptr<const IDiscreteFunction>& source,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
 {
-  const std::shared_ptr i_mesh      = getCommonMesh({U, source});
-  const std::shared_ptr i_dual_mesh = getCommonMesh({dual_lambdab, dual_mub});
-  checkDiscretizationType({dual_lambdab, dual_mub, U, source}, DiscreteFunctionType::P0);
+  const std::shared_ptr i_mesh = getCommonMesh({U, source});
+  if (not i_mesh) {
+    throw NormalError("primal discrete functions are not defined on the same mesh");
+  }
+  const std::shared_ptr i_dual_mesh = getCommonMesh({dual_lambdab, dual_mub, dual_U});
+  if (not i_dual_mesh) {
+    throw NormalError("dual discrete functions are not defined on the same mesh");
+  }
+  checkDiscretizationType({dual_lambdab, dual_mub, dual_U, U, source}, DiscreteFunctionType::P0);
 
   switch (i_mesh->dimension()) {
   case 1: {
@@ -1964,8 +1970,8 @@ EnergyComputerHandler::EnergyComputerHandler(
 
     std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
 
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != dual_mub->mesh()) {
-      throw NormalError("mu_dual is not defined on the dual mesh");
+    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
+      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
     }
 
     m_energy_computer =
@@ -1985,8 +1991,8 @@ EnergyComputerHandler::EnergyComputerHandler(
 
     std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
 
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != dual_mub->mesh()) {
-      throw NormalError("mu_dual is not defined on the dual mesh");
+    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
+      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
     }
 
     m_energy_computer =
@@ -2006,8 +2012,8 @@ EnergyComputerHandler::EnergyComputerHandler(
 
     std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
 
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != dual_mub->mesh()) {
-      throw NormalError("mu_dual is not defined on the dual mesh");
+    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
+      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
     }
 
     m_energy_computer =