diff --git a/src/main.cpp b/src/main.cpp
index cbaf6277e5ee0c1c7d689403e8d78cae54c3ed5f..ed341a0a4e37c41523f8cb761e83f48ff2a48308 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -7,6 +7,7 @@
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/StencilManager.hpp>
 #include <mesh/SynchronizerManager.hpp>
+#include <scheme/PolynomialBasisDataManager.hpp>
 #include <utils/ExecutionStatManager.hpp>
 #include <utils/GlobalVariableManager.hpp>
 #include <utils/PugsUtils.hpp>
@@ -32,6 +33,7 @@ main(int argc, char* argv[])
   RandomEngine::create();
   QuadratureManager::create();
   MeshDataManager::create();
+  PolynomialBasisDataManager::create();
   DualConnectivityManager::create();
   DualMeshManager::create();
   StencilManager::create();
@@ -44,6 +46,7 @@ main(int argc, char* argv[])
   StencilManager::destroy();
   DualMeshManager::destroy();
   DualConnectivityManager::destroy();
+  PolynomialBasisDataManager::destroy();
   MeshDataManager::destroy();
   QuadratureManager::destroy();
   RandomEngine::destroy();
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 68c6de46daaac99375fcb6b1d3d4e61bffdd83c8..be5b069773407c9010b0de5d81edd33e55cac5d4 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -54,3 +54,9 @@ target_link_libraries(
   ${PARMETIS_TARGET}
   ${HIGHFIVE_TARGET}
 )
+
+# Additional dependencies
+add_dependencies(
+  PugsMesh
+  PugsScheme
+)
diff --git a/src/mesh/Mesh.cpp b/src/mesh/Mesh.cpp
index baa0d00d71e7371514aad465870b7fafb0d3f460..6e368af538ba3edb5fe0b5f222743ce6c2ebfacf 100644
--- a/src/mesh/Mesh.cpp
+++ b/src/mesh/Mesh.cpp
@@ -3,11 +3,13 @@
 #include <mesh/Connectivity.hpp>
 #include <mesh/DualMeshManager.hpp>
 #include <mesh/MeshDataManager.hpp>
+#include <scheme/PolynomialBasisDataManager.hpp>
 
 template <size_t Dimension>
 Mesh<Dimension>::~Mesh()
 {
   MeshDataManager::instance().deleteMeshData(this->id());
+  PolynomialBasisDataManager::instance().deletePolynomialBasisData(this->id());
   DualMeshManager::instance().deleteMesh(this->id());
 }
 
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 852a6cd751d092a7ab1c59437c1e549ae2f8094f..7e5e0bf6226e8f17aa5ae38368b4f1a096815ced 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -13,6 +13,7 @@ add_library(
   FluxingAdvectionSolver.cpp
   HyperelasticSolver.cpp
   LoadBalancer.cpp
+  PolynomialBasisDataManager.cpp
   PolynomialReconstruction.cpp
 )
 
diff --git a/src/scheme/PolynomialBasisData.hpp b/src/scheme/PolynomialBasisData.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..43640b00eeb38464120cd113a6ab24b6aa3d5853
--- /dev/null
+++ b/src/scheme/PolynomialBasisData.hpp
@@ -0,0 +1,46 @@
+#ifndef POLYNOMIAL_BASIS_DATA_HPP
+#define POLYNOMIAL_BASIS_DATA_HPP
+
+#include <mesh/MeshTraits.hpp>
+
+template <size_t Dimension>
+class Mesh;
+
+template <MeshConcept MeshType>
+class PolynomialBasisData;
+
+template <size_t Dimension>
+class PolynomialBasisData<Mesh<Dimension>>
+{
+ public:
+  using MeshType = Mesh<Dimension>;
+
+  static_assert(Dimension > 0, "dimension must be strictly positive");
+  static_assert((Dimension <= 3), "only 1d, 2d and 3d are implemented");
+
+ private:
+  const MeshType& m_mesh;
+
+ public:
+  PUGS_INLINE
+  const MeshType&
+  mesh() const
+  {
+    return m_mesh;
+  }
+
+ private:
+  // PolynomialBasisData **must** be constructed through PolynomialBasisDataManager
+  friend class PolynomialBasisDataManager;
+  PolynomialBasisData(const MeshType& mesh) : m_mesh(mesh) {}
+
+ public:
+  PolynomialBasisData() = delete;
+
+  PolynomialBasisData(const PolynomialBasisData&) = delete;
+  PolynomialBasisData(PolynomialBasisData&&)      = delete;
+
+  ~PolynomialBasisData() {}
+};
+
+#endif   // POLYNOMIAL_BASIS_DATA_HPP
diff --git a/src/scheme/PolynomialBasisDataManager.cpp b/src/scheme/PolynomialBasisDataManager.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8f659e9230dce54c7c607342e87110c58692a971
--- /dev/null
+++ b/src/scheme/PolynomialBasisDataManager.cpp
@@ -0,0 +1,63 @@
+#include <utils/PugsAssert.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <scheme/PolynomialBasisData.hpp>
+#include <scheme/PolynomialBasisDataManager.hpp>
+#include <scheme/PolynomialBasisDataVariant.hpp>
+#include <utils/Exceptions.hpp>
+
+#include <sstream>
+
+PolynomialBasisDataManager* PolynomialBasisDataManager::m_instance{nullptr};
+
+void
+PolynomialBasisDataManager::create()
+{
+  Assert(m_instance == nullptr, "PolynomialBasisDataManager is already created");
+  m_instance = new PolynomialBasisDataManager;
+}
+
+void
+PolynomialBasisDataManager::destroy()
+{
+  Assert(m_instance != nullptr, "PolynomialBasisDataManager was not created!");
+
+  if (m_instance->m_mesh_id_mesh_polynomial_data_map.size() > 0) {
+    std::stringstream error;
+    error << ": some mesh data is still registered\n";
+    for (const auto& i_mesh_data : m_instance->m_mesh_id_mesh_polynomial_data_map) {
+      error << " - mesh data " << rang::fgB::magenta << i_mesh_data.first << rang::style::reset << '\n';
+    }
+    throw UnexpectedError(error.str());
+  }
+  delete m_instance;
+  m_instance = nullptr;
+}
+
+void
+PolynomialBasisDataManager::deletePolynomialBasisData(const size_t mesh_id)
+{
+  m_mesh_id_mesh_polynomial_data_map.erase(mesh_id);
+}
+
+template <MeshConcept MeshType>
+PolynomialBasisData<MeshType>&
+PolynomialBasisDataManager::getPolynomialBasisData(const MeshType& mesh)
+{
+  if (auto i_mesh_data = m_mesh_id_mesh_polynomial_data_map.find(mesh.id());
+      i_mesh_data != m_mesh_id_mesh_polynomial_data_map.end()) {
+    const auto& mesh_data_v = *i_mesh_data->second;
+    return *mesh_data_v.template get<MeshType>();
+  } else {
+    // **cannot** use make_shared since PolynomialBasisData constructor is **private**
+    std::shared_ptr<PolynomialBasisData<MeshType>> mesh_data{new PolynomialBasisData<MeshType>(mesh)};
+
+    m_mesh_id_mesh_polynomial_data_map[mesh.id()] = std::make_shared<PolynomialBasisDataVariant>(mesh_data);
+    return *mesh_data;
+  }
+}
+
+template PolynomialBasisData<Mesh<1>>& PolynomialBasisDataManager::getPolynomialBasisData(const Mesh<1>&);
+template PolynomialBasisData<Mesh<2>>& PolynomialBasisDataManager::getPolynomialBasisData(const Mesh<2>&);
+template PolynomialBasisData<Mesh<3>>& PolynomialBasisDataManager::getPolynomialBasisData(const Mesh<3>&);
diff --git a/src/scheme/PolynomialBasisDataManager.hpp b/src/scheme/PolynomialBasisDataManager.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..43b96dce8a1ab7ff1248316d9d061fadb6d2253a
--- /dev/null
+++ b/src/scheme/PolynomialBasisDataManager.hpp
@@ -0,0 +1,47 @@
+#ifndef POLYNOMIAL_BASIS_DATA_MANAGER_HPP
+#define POLYNOMIAL_BASIS_DATA_MANAGER_HPP
+
+#include <mesh/MeshTraits.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+
+#include <memory>
+#include <unordered_map>
+
+class PolynomialBasisDataVariant;
+
+template <MeshConcept MeshType>
+class PolynomialBasisData;
+
+class PolynomialBasisDataManager
+{
+ private:
+  std::unordered_map<size_t, std::shared_ptr<PolynomialBasisDataVariant>> m_mesh_id_mesh_polynomial_data_map;
+
+  static PolynomialBasisDataManager* m_instance;
+
+  PolynomialBasisDataManager(const PolynomialBasisDataManager&) = delete;
+  PolynomialBasisDataManager(PolynomialBasisDataManager&&)      = delete;
+
+  PolynomialBasisDataManager()  = default;
+  ~PolynomialBasisDataManager() = default;
+
+ public:
+  static void create();
+  static void destroy();
+
+  PUGS_INLINE
+  static PolynomialBasisDataManager&
+  instance()
+  {
+    Assert(m_instance != nullptr, "PolynomialBasisDataManager was not created!");
+    return *m_instance;
+  }
+
+  void deletePolynomialBasisData(const size_t mesh_id);
+
+  template <MeshConcept MeshType>
+  PolynomialBasisData<MeshType>& getPolynomialBasisData(const MeshType&);
+};
+
+#endif   // POLYNOMIAL_BASIS_DATA_MANAGER_HPP
diff --git a/src/scheme/PolynomialBasisDataVariant.hpp b/src/scheme/PolynomialBasisDataVariant.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b713350b6fd2a352a49c578143f0c00c8caa8167
--- /dev/null
+++ b/src/scheme/PolynomialBasisDataVariant.hpp
@@ -0,0 +1,64 @@
+#ifndef POLYNOMIAL_BASIS_DATA_VARIANT_HPP
+#define POLYNOMIAL_BASIS_DATA_VARIANT_HPP
+
+#include <mesh/MeshTraits.hpp>
+#include <scheme/PolynomialBasisData.hpp>
+#include <utils/Demangle.hpp>
+#include <utils/Exceptions.hpp>
+#include <utils/PugsMacros.hpp>
+
+#include <rang.hpp>
+
+#include <memory>
+#include <sstream>
+#include <variant>
+
+template <size_t Dimension>
+class Mesh;
+
+class PolynomialBasisDataVariant
+{
+ private:
+  using Variant = std::variant<std::shared_ptr<PolynomialBasisData<Mesh<1>>>,   //
+                               std::shared_ptr<PolynomialBasisData<Mesh<2>>>,   //
+                               std::shared_ptr<PolynomialBasisData<Mesh<3>>>>;
+
+  Variant m_p_mesh_data_variant;
+
+ public:
+  template <MeshConcept MeshType>
+  PUGS_INLINE std::shared_ptr<PolynomialBasisData<MeshType>>
+  get() const
+  {
+    if (not std::holds_alternative<std::shared_ptr<PolynomialBasisData<MeshType>>>(this->m_p_mesh_data_variant)) {
+      std::ostringstream error_msg;
+      error_msg << "invalid mesh type type\n";
+      error_msg << "- required " << rang::fgB::red << demangle<PolynomialBasisData<MeshType>>() << rang::fg::reset
+                << '\n';
+      error_msg << "- contains " << rang::fgB::yellow
+                << std::visit(
+                     [](auto&& p_mesh_data) -> std::string {
+                       using FoundPolynomialBasisDataType = typename std::decay_t<decltype(p_mesh_data)>::element_type;
+                       return demangle<FoundPolynomialBasisDataType>();
+                     },
+                     this->m_p_mesh_data_variant)
+                << rang::fg::reset;
+      throw NormalError(error_msg.str());
+    }
+    return std::get<std::shared_ptr<PolynomialBasisData<MeshType>>>(m_p_mesh_data_variant);
+  }
+
+  PolynomialBasisDataVariant() = delete;
+
+  template <MeshConcept MeshType>
+  PolynomialBasisDataVariant(const std::shared_ptr<PolynomialBasisData<MeshType>>& p_mesh_data)
+    : m_p_mesh_data_variant{p_mesh_data}
+  {}
+
+  PolynomialBasisDataVariant(const PolynomialBasisDataVariant&) = default;
+  PolynomialBasisDataVariant(PolynomialBasisDataVariant&&)      = default;
+
+  ~PolynomialBasisDataVariant() = default;
+};
+
+#endif   // POLYNOMIAL_BASIS_DATA_VARIANT_HPP
diff --git a/tests/mpi_test_main.cpp b/tests/mpi_test_main.cpp
index cabb10386d4e8e99ab0dd9b1946bc7dba818502f..36a28b131ec0f578744400a75e631eff6d6c0e35 100644
--- a/tests/mpi_test_main.cpp
+++ b/tests/mpi_test_main.cpp
@@ -9,6 +9,7 @@
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/StencilManager.hpp>
 #include <mesh/SynchronizerManager.hpp>
+#include <scheme/PolynomialBasisDataManager.hpp>
 #include <utils/GlobalVariableManager.hpp>
 #include <utils/Messenger.hpp>
 #include <utils/PETScWrapper.hpp>
@@ -111,6 +112,7 @@ main(int argc, char* argv[])
       RandomEngine::create();
       QuadratureManager::create();
       MeshDataManager::create();
+      PolynomialBasisDataManager::create();
       DualConnectivityManager::create();
       DualMeshManager::create();
       StencilManager::create();
@@ -132,6 +134,7 @@ main(int argc, char* argv[])
       GlobalVariableManager::destroy();
       DualMeshManager::destroy();
       DualConnectivityManager::destroy();
+      PolynomialBasisDataManager::destroy();
       MeshDataManager::destroy();
       QuadratureManager::destroy();
       RandomEngine::destroy();
diff --git a/tests/test_main.cpp b/tests/test_main.cpp
index 98c01f7879d00fbf33718ab679d05ef3d7c73fde..83ff999d5602c4267c68977178d0dc5e0fad4721 100644
--- a/tests/test_main.cpp
+++ b/tests/test_main.cpp
@@ -9,6 +9,7 @@
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/StencilManager.hpp>
 #include <mesh/SynchronizerManager.hpp>
+#include <scheme/PolynomialBasisDataManager.hpp>
 #include <utils/GlobalVariableManager.hpp>
 #include <utils/Messenger.hpp>
 #include <utils/PETScWrapper.hpp>
@@ -66,6 +67,7 @@ main(int argc, char* argv[])
       RandomEngine::create();
       QuadratureManager::create();
       MeshDataManager::create();
+      PolynomialBasisDataManager::create();
       DualConnectivityManager::create();
       DualMeshManager::create();
       StencilManager::create();
@@ -87,6 +89,7 @@ main(int argc, char* argv[])
       StencilManager::destroy();
       DualMeshManager::destroy();
       DualConnectivityManager::destroy();
+      PolynomialBasisDataManager::destroy();
       MeshDataManager::destroy();
       QuadratureManager::destroy();
       RandomEngine::destroy();