diff --git a/src/mesh/DiamondDualConnectivityBuilder.cpp b/src/mesh/DiamondDualConnectivityBuilder.cpp
index 8bcbe9bad8e15cbef1dfa449c0d6defc115f126f..e20f972c4adecd1e2d449e82044357bac03c9ef4 100644
--- a/src/mesh/DiamondDualConnectivityBuilder.cpp
+++ b/src/mesh/DiamondDualConnectivityBuilder.cpp
@@ -278,7 +278,6 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
           node_array[i] = diamond_node_list[i];
         }
         diamond_descriptor.addRefItemList(RefNodeList{primal_ref_node_list.refId(), node_array});
-        std::cout << "stored " << primal_ref_node_list.refId() << '\n';
       }
     }
   }
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index 3865414c1dcf51e42dbf19aa69e4eb63474e2ce1..c959899f0a20fc97e37a01b493d18b9b6e61e5d2 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -11,7 +11,96 @@
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/RefId.hpp>
 #include <utils/Array.hpp>
+#include <utils/CastArray.hpp>
 #include <utils/Messenger.hpp>
+#include <utils/PugsAssert.hpp>
+
+template <size_t Dimension>
+class MeshToDualDataMapper
+{
+ private:
+  const IConnectivity* m_primal_connectivity;
+  const IConnectivity* m_dual_connectivity;
+
+  using NodeIdToNodeIdMap = Array<std::pair<NodeId, NodeId>>;
+  NodeIdToNodeIdMap m_primal_node_to_dual_node_map;
+
+  using CellIdToNodeIdMap = Array<std::pair<CellId, NodeId>>;
+  CellIdToNodeIdMap m_primal_cell_to_dual_node_map;
+
+ public:
+  template <typename DataType>
+  void
+  toDualNode(const NodeValue<const DataType>& primal_node_value,
+             const CellValue<const DataType>& primal_cell_value,
+             NodeValue<DataType>& dual_node_value)
+  {
+    Assert(m_primal_connectivity == primal_cell_value.connectivity_ptr().get());
+    Assert(m_primal_connectivity == primal_node_value.connectivity_ptr().get());
+    Assert(m_dual_connectivity == dual_node_value.connectivity_ptr().get());
+
+    parallel_for(
+      m_primal_node_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) {
+        const auto [primal_node_id, dual_node_id] = m_primal_node_to_dual_node_map[i];
+
+        dual_node_value[dual_node_id] = primal_node_value[primal_node_id];
+      });
+
+    parallel_for(
+      m_primal_cell_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) {
+        const auto [primal_cell_id, dual_node_id] = m_primal_cell_to_dual_node_map[i];
+        dual_node_value[dual_node_id]             = primal_cell_value[primal_cell_id];
+      });
+  }
+
+  MeshToDualDataMapper(const Connectivity<Dimension>& primal_connectivity,
+                       const Connectivity<Dimension>& dual_connectivity)
+    : m_primal_connectivity{&primal_connectivity}, m_dual_connectivity{&dual_connectivity}
+  {
+    if constexpr (Dimension == 1) {
+      const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
+
+      throw UnexpectedError("Looks like diamond mesh construction is buggy in 1D");
+
+      NodeId dual_node_id            = 0;
+      m_primal_node_to_dual_node_map = [&]() {
+        std::vector<std::pair<NodeId, NodeId>> primal_node_to_dual_node_vector;
+        for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
+          if (node_to_cell_matrix[primal_node_id].size() == 1) {
+            primal_node_to_dual_node_vector.push_back(std::make_pair(primal_node_id, dual_node_id++));
+          }
+        }
+        return convert_to_array(primal_node_to_dual_node_vector);
+      }();
+
+      m_primal_cell_to_dual_node_map = [&]() {
+        CellIdToNodeIdMap primal_cell_to_dual_node_map{primal_connectivity.numberOfCells()};
+        for (CellId cell_id = 0; cell_id < primal_cell_to_dual_node_map.size(); ++cell_id) {
+          primal_cell_to_dual_node_map[cell_id] = std::make_pair(cell_id, dual_node_id++);
+        }
+        return primal_cell_to_dual_node_map;
+      }();
+
+    } else {
+      m_primal_node_to_dual_node_map = [&]() {
+        NodeIdToNodeIdMap primal_node_to_dual_node_map{primal_connectivity.numberOfNodes()};
+        for (NodeId node_id = 0; node_id < primal_node_to_dual_node_map.size(); ++node_id) {
+          primal_node_to_dual_node_map[node_id] = std::make_pair(node_id, node_id);
+        }
+        return primal_node_to_dual_node_map;
+      }();
+
+      m_primal_cell_to_dual_node_map = [&]() {
+        CellIdToNodeIdMap primal_cell_to_dual_node_map{primal_connectivity.numberOfCells()};
+        NodeId dual_node_id = m_primal_node_to_dual_node_map.size();
+        for (CellId cell_id = 0; cell_id < primal_cell_to_dual_node_map.size(); ++cell_id) {
+          primal_cell_to_dual_node_map[cell_id] = std::make_pair(cell_id, dual_node_id++);
+        }
+        return primal_cell_to_dual_node_map;
+      }();
+    }
+  }
+};
 
 template <size_t Dimension>
 void
@@ -24,56 +113,41 @@ DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(
 
   const ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
 
-  NodeValue<TinyVector<Dimension>> diamond_xr{diamond_connectivity};
-
   const NodeValue<const TinyVector<Dimension>> primal_xr = primal_mesh.xr();
+
   MeshData<Dimension>& primal_mesh_data                  = MeshDataManager::instance().getMeshData(primal_mesh);
   const CellValue<const TinyVector<Dimension>> primal_xj = primal_mesh_data.xj();
 
-  NodeId i_node = 0;
-  for (; i_node < primal_mesh.numberOfNodes(); ++i_node) {
-    diamond_xr[i_node] = primal_xr[i_node];
-  }
+  MeshToDualDataMapper<Dimension> mesh_to_dual_data_mapper{primal_mesh.connectivity(), diamond_connectivity};
 
-  for (CellId i_cell = 0; i_cell < primal_mesh.numberOfCells(); ++i_cell) {
-    diamond_xr[i_node++] = primal_xj[i_cell];
-  }
+  NodeValue<TinyVector<Dimension>> diamond_xr{diamond_connectivity};
+  mesh_to_dual_data_mapper.toDualNode(primal_xr, primal_xj, diamond_xr);
 
   m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
 }
 
-template <>
-void
-DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(const Mesh<Connectivity<1>>& primal_mesh,
-                                                  const std::shared_ptr<const Connectivity<1>>& p_diamond_connectivity)
-{
-  using ConnectivityType = Connectivity<1>;
-  using MeshType         = Mesh<Connectivity<1>>;
+// template <>
+// void
+// DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(const Mesh<Connectivity<1>>& primal_mesh,
+//                                                   const std::shared_ptr<const Connectivity<1>>&
+//                                                   p_diamond_connectivity)
+// {
+//   using ConnectivityType = Connectivity<1>;
+//   using MeshType         = Mesh<Connectivity<1>>;
 
-  const ConnectivityType& primal_connectivity = primal_mesh.connectivity();
-  const auto& node_to_cell_matrix             = primal_connectivity.nodeToCellMatrix();
+//   const ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
 
-  const ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
+//   const NodeValue<const TinyVector<1>> primal_xr = primal_mesh.xr();
+//   MeshData<1>& primal_mesh_data                  = MeshDataManager::instance().getMeshData(primal_mesh);
+//   const CellValue<const TinyVector<1>> primal_xj = primal_mesh_data.xj();
 
-  NodeValue<TinyVector<1>> diamond_xr{diamond_connectivity};
+//   MeshToDualDataMapper<1> mesh_to_dual_data_mapper{primal_mesh.connectivity(), diamond_connectivity};
 
-  const NodeValue<const TinyVector<1>> primal_xr = primal_mesh.xr();
-  MeshData<1>& primal_mesh_data                  = MeshDataManager::instance().getMeshData(primal_mesh);
-  const CellValue<const TinyVector<1>> primal_xj = primal_mesh_data.xj();
+//   NodeValue<TinyVector<1>> diamond_xr{diamond_connectivity};
+//   mesh_to_dual_data_mapper.toDualNode(primal_xr, primal_xj, diamond_xr);
 
-  NodeId next_node_id = 0;
-  for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
-    if (node_to_cell_matrix[primal_node_id].size() == 1) {
-      diamond_xr[next_node_id++] = primal_xr[primal_node_id];
-    }
-  }
-
-  for (CellId i_cell = 0; i_cell < primal_mesh.numberOfCells(); ++i_cell) {
-    diamond_xr[next_node_id++] = primal_xj[i_cell];
-  }
-
-  m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
-}
+//   m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
+// }
 
 DiamondDualMeshBuilder::DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>& p_mesh)
 {