diff --git a/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp
index 5d000813aeff1e670a2f547a69827a679cbfa699..892037219c389e6ff7ebd28e8b5afd586845c14f 100644
--- a/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp
@@ -234,10 +234,14 @@ dot(const std::shared_ptr<const DiscreteFunctionVariant>& f_v,
             throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
           }
         } else if constexpr (is_discrete_function_P0_vector_v<TypeOfF>) {
-          if (f.size() == g.size()) {
-            return std::make_shared<DiscreteFunctionVariant>(dot(f, g));
+          if constexpr (std::is_arithmetic_v<DataType>) {
+            if (f.size() == g.size()) {
+              return std::make_shared<DiscreteFunctionVariant>(dot(f, g));
+            } else {
+              throw NormalError("operands have different dimension");
+            }
           } else {
-            throw NormalError("operands have different dimension");
+            throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
           }
         } else {
           throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
@@ -691,8 +695,6 @@ sum_of_Vh_components(const std::shared_ptr<const DiscreteFunctionVariant>& f)
     [&](auto&& discrete_function) -> std::shared_ptr<const DiscreteFunctionVariant> {
       using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
       if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
-        using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
-        static_assert(std::is_same_v<DataType, double>);
         return std::make_shared<DiscreteFunctionVariant>(sumOfComponents(discrete_function));
       } else {
         throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
diff --git a/src/output/OutputNamedItemValueSet.hpp b/src/output/OutputNamedItemValueSet.hpp
index aa2c5e2ed67140d47c7d4ddf1f57caecc2f28360..8c43330420c91417e0c339228b4bfb41f9dc908b 100644
--- a/src/output/OutputNamedItemValueSet.hpp
+++ b/src/output/OutputNamedItemValueSet.hpp
@@ -38,7 +38,7 @@ class NamedItemData
   }
 
   NamedItemData& operator=(const NamedItemData&) = default;
-  NamedItemData& operator=(NamedItemData&&) = default;
+  NamedItemData& operator=(NamedItemData&&)      = default;
 
   NamedItemData(const std::string& name, const ItemDataT<DataType, item_type, ConnectivityPtr>& item_data)
     : m_name(name), m_item_data(item_data)
@@ -113,24 +113,48 @@ class OutputNamedItemDataSet
                                        NodeArray<const long int>,
                                        NodeArray<const unsigned long int>,
                                        NodeArray<const double>,
+                                       NodeArray<const TinyVector<1, double>>,
+                                       NodeArray<const TinyVector<2, double>>,
+                                       NodeArray<const TinyVector<3, double>>,
+                                       NodeArray<const TinyMatrix<1, 1, double>>,
+                                       NodeArray<const TinyMatrix<2, 2, double>>,
+                                       NodeArray<const TinyMatrix<3, 3, double>>,
 
                                        EdgeArray<const bool>,
                                        EdgeArray<const int>,
                                        EdgeArray<const long int>,
                                        EdgeArray<const unsigned long int>,
                                        EdgeArray<const double>,
+                                       EdgeArray<const TinyVector<1, double>>,
+                                       EdgeArray<const TinyVector<2, double>>,
+                                       EdgeArray<const TinyVector<3, double>>,
+                                       EdgeArray<const TinyMatrix<1, 1, double>>,
+                                       EdgeArray<const TinyMatrix<2, 2, double>>,
+                                       EdgeArray<const TinyMatrix<3, 3, double>>,
 
                                        FaceArray<const bool>,
                                        FaceArray<const int>,
                                        FaceArray<const long int>,
                                        FaceArray<const unsigned long int>,
                                        FaceArray<const double>,
+                                       FaceArray<const TinyVector<1, double>>,
+                                       FaceArray<const TinyVector<2, double>>,
+                                       FaceArray<const TinyVector<3, double>>,
+                                       FaceArray<const TinyMatrix<1, 1, double>>,
+                                       FaceArray<const TinyMatrix<2, 2, double>>,
+                                       FaceArray<const TinyMatrix<3, 3, double>>,
 
                                        CellArray<const bool>,
                                        CellArray<const int>,
                                        CellArray<const long int>,
                                        CellArray<const unsigned long int>,
-                                       CellArray<const double>>;
+                                       CellArray<const double>,
+                                       CellArray<const TinyVector<1, double>>,
+                                       CellArray<const TinyVector<2, double>>,
+                                       CellArray<const TinyVector<3, double>>,
+                                       CellArray<const TinyMatrix<1, 1, double>>,
+                                       CellArray<const TinyMatrix<2, 2, double>>,
+                                       CellArray<const TinyMatrix<3, 3, double>>>;
 
  private:
   // We do not use a map, we want variables to be written in the
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index 64d35434462e5980a11b69ea2257ce779b3ec581..2aedce3b1fe264bf3a11a94592077f4631e666aa 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -396,10 +396,17 @@ VTKWriter::_write(const MeshType& mesh,
          << "\">\n";
     fout << "<CellData>\n";
     for (const auto& [name, item_value_variant] : output_named_item_data_set) {
-      std::visit([&, var_name = name](
-                   auto&&
-                     item_value) { return this->_write_cell_data(fout, var_name, item_value, serialize_data_list); },
-                 item_value_variant);
+      std::visit(
+        [&, var_name = name](auto&& item_value) {
+          using IVType   = std::decay_t<decltype(item_value)>;
+          using DataType = typename IVType::data_type;
+          if constexpr (is_item_array_v<IVType> and not std::is_arithmetic_v<DataType>) {
+            throw NotImplementedError("DiscreteFunctionP0Vector of non arithmetic type");
+          } else {
+            return this->_write_cell_data(fout, var_name, item_value, serialize_data_list);
+          }
+        },
+        item_value_variant);
     }
     if (parallel::size() > 1) {
       CellValue<uint8_t> vtk_ghost_type{mesh.connectivity()};
@@ -413,10 +420,17 @@ VTKWriter::_write(const MeshType& mesh,
     fout << "</CellData>\n";
     fout << "<PointData>\n";
     for (const auto& [name, item_value_variant] : output_named_item_data_set) {
-      std::visit([&, var_name = name](
-                   auto&&
-                     item_value) { return this->_write_node_data(fout, var_name, item_value, serialize_data_list); },
-                 item_value_variant);
+      std::visit(
+        [&, var_name = name](auto&& item_value) {
+          using IVType   = std::decay_t<decltype(item_value)>;
+          using DataType = typename IVType::data_type;
+          if constexpr (is_item_array_v<IVType> and not std::is_arithmetic_v<DataType>) {
+            throw NotImplementedError("DiscreteFunctionP0Vector of non arithmetic type");
+          } else {
+            return this->_write_node_data(fout, var_name, item_value, serialize_data_list);
+          }
+        },
+        item_value_variant);
     }
     fout << "</PointData>\n";
     fout << "<Points>\n";
@@ -644,15 +658,33 @@ VTKWriter::_write(const MeshType& mesh,
 
     fout << "<PPointData>\n";
     for (const auto& [name, item_value_variant] : output_named_item_data_set) {
-      std::visit([&, var_name = name](auto&& item_value) { return this->_write_node_pvtu(fout, var_name, item_value); },
-                 item_value_variant);
+      std::visit(
+        [&, var_name = name](auto&& item_value) {
+          using IVType   = std::decay_t<decltype(item_value)>;
+          using DataType = typename IVType::data_type;
+          if constexpr (is_item_array_v<IVType> and not std::is_arithmetic_v<DataType>) {
+            throw NotImplementedError("DiscreteFunctionP0Vector of non arithmetic type");
+          } else {
+            return this->_write_node_pvtu(fout, var_name, item_value);
+          }
+        },
+        item_value_variant);
     }
     fout << "</PPointData>\n";
 
     fout << "<PCellData>\n";
     for (const auto& [name, item_value_variant] : output_named_item_data_set) {
-      std::visit([&, var_name = name](auto&& item_value) { return this->_write_cell_pvtu(fout, var_name, item_value); },
-                 item_value_variant);
+      std::visit(
+        [&, var_name = name](auto&& item_value) {
+          using IVType   = std::decay_t<decltype(item_value)>;
+          using DataType = typename IVType::data_type;
+          if constexpr (is_item_array_v<IVType> and not std::is_arithmetic_v<DataType>) {
+            throw NotImplementedError("DiscreteFunctionP0Vector of non arithmetic type");
+          } else {
+            return this->_write_cell_pvtu(fout, var_name, item_value);
+          }
+        },
+        item_value_variant);
     }
     if (parallel::size() > 1) {
       fout << "<PDataArray type=\"UInt8\" Name=\"vtkGhostType\" NumberOfComponents=\"1\"/>\n";
diff --git a/src/scheme/DiscreteFunctionDPkVariant.hpp b/src/scheme/DiscreteFunctionDPkVariant.hpp
index 0d683c1783c9aad299a81dcc4d7ef3817470d54a..2942f141f38452faed9e95a5f95ca675bc7d8c29 100644
--- a/src/scheme/DiscreteFunctionDPkVariant.hpp
+++ b/src/scheme/DiscreteFunctionDPkVariant.hpp
@@ -35,8 +35,28 @@ class DiscreteFunctionDPkVariant
                                DiscreteFunctionDPk<3, const TinyMatrix<3>>,
 
                                DiscreteFunctionDPkVector<1, const double>,
+                               DiscreteFunctionDPkVector<1, const TinyVector<1>>,
+                               DiscreteFunctionDPkVector<1, const TinyVector<2>>,
+                               DiscreteFunctionDPkVector<1, const TinyVector<3>>,
+                               DiscreteFunctionDPkVector<1, const TinyMatrix<1>>,
+                               DiscreteFunctionDPkVector<1, const TinyMatrix<2>>,
+                               DiscreteFunctionDPkVector<1, const TinyMatrix<3>>,
+
                                DiscreteFunctionDPkVector<2, const double>,
-                               DiscreteFunctionDPkVector<3, const double>>;
+                               DiscreteFunctionDPkVector<2, const TinyVector<1>>,
+                               DiscreteFunctionDPkVector<2, const TinyVector<2>>,
+                               DiscreteFunctionDPkVector<2, const TinyVector<3>>,
+                               DiscreteFunctionDPkVector<2, const TinyMatrix<1>>,
+                               DiscreteFunctionDPkVector<2, const TinyMatrix<2>>,
+                               DiscreteFunctionDPkVector<2, const TinyMatrix<3>>,
+
+                               DiscreteFunctionDPkVector<3, const double>,
+                               DiscreteFunctionDPkVector<3, const TinyVector<1>>,
+                               DiscreteFunctionDPkVector<3, const TinyVector<2>>,
+                               DiscreteFunctionDPkVector<3, const TinyVector<3>>,
+                               DiscreteFunctionDPkVector<3, const TinyMatrix<1>>,
+                               DiscreteFunctionDPkVector<3, const TinyMatrix<2>>,
+                               DiscreteFunctionDPkVector<3, const TinyMatrix<3>>>;
 
  private:
   Variant m_discrete_function_dpk;
@@ -91,7 +111,13 @@ class DiscreteFunctionDPkVariant
   DiscreteFunctionDPkVariant(const DiscreteFunctionDPkVector<Dimension, DataType>& discrete_function_dpk)
     : m_discrete_function_dpk{DiscreteFunctionDPkVector<Dimension, const DataType>{discrete_function_dpk}}
   {
-    static_assert(std::is_same_v<std::remove_const_t<DataType>, double>,
+    static_assert(std::is_same_v<std::remove_const_t<DataType>, double> or                       //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<1, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<2, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<3, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<1, 1, double>> or   //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<2, 2, double>> or   //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<3, 3, double>>,
                   "DiscreteFunctionDPkVector with this DataType is not allowed in variant");
   }
 
diff --git a/src/scheme/DiscreteFunctionP0Vector.hpp b/src/scheme/DiscreteFunctionP0Vector.hpp
index da3904b2ad9853a05a521920c3a3a9434aac14a9..bc08cb2a7f27286b9bafdd8bb15d0c761b892d1b 100644
--- a/src/scheme/DiscreteFunctionP0Vector.hpp
+++ b/src/scheme/DiscreteFunctionP0Vector.hpp
@@ -19,8 +19,6 @@ class DiscreteFunctionP0Vector
   friend class DiscreteFunctionP0Vector<std::add_const_t<DataType>>;
   friend class DiscreteFunctionP0Vector<std::remove_const_t<DataType>>;
 
-  static_assert(std::is_arithmetic_v<DataType>, "DiscreteFunctionP0Vector are only defined for arithmetic data type");
-
  private:
   std::shared_ptr<const MeshVariant> m_mesh;
   CellArray<DataType> m_cell_arrays;
@@ -188,16 +186,23 @@ class DiscreteFunctionP0Vector
     return product;
   }
 
-  PUGS_INLINE friend DiscreteFunctionP0<double>
+  PUGS_INLINE friend DiscreteFunctionP0<std::remove_const_t<DataType>>
   sumOfComponents(const DiscreteFunctionP0Vector& f)
   {
-    DiscreteFunctionP0<double> result{f.m_mesh};
+    DiscreteFunctionP0<std::remove_const_t<DataType>> result{f.m_mesh};
 
     parallel_for(
       f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
         const auto& f_cell_id = f[cell_id];
 
-        double sum = 0;
+        std::remove_const_t<DataType> sum = [] {
+          if constexpr (std::is_arithmetic_v<DataType>) {
+            return 0;
+          } else {
+            return zero;
+          }
+        }();
+
         for (size_t i = 0; i < f.size(); ++i) {
           sum += f_cell_id[i];
         }
@@ -213,6 +218,7 @@ class DiscreteFunctionP0Vector
   {
     Assert(f.meshVariant()->id() == g.meshVariant()->id(), "functions are nor defined on the same mesh");
     Assert(f.size() == g.size());
+    static_assert(std::is_arithmetic_v<std::decay_t<DataType>>);
     DiscreteFunctionP0<double> result{f.m_mesh};
     parallel_for(
       f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
diff --git a/src/scheme/DiscreteFunctionVariant.hpp b/src/scheme/DiscreteFunctionVariant.hpp
index 901b7989084249ca0e437f9eef24e0371dd723a8..616184298a212c10fe4865b451d0a7879a18b2ea 100644
--- a/src/scheme/DiscreteFunctionVariant.hpp
+++ b/src/scheme/DiscreteFunctionVariant.hpp
@@ -19,7 +19,13 @@ class DiscreteFunctionVariant
                                DiscreteFunctionP0<const TinyMatrix<2>>,
                                DiscreteFunctionP0<const TinyMatrix<3>>,
 
-                               DiscreteFunctionP0Vector<const double>>;
+                               DiscreteFunctionP0Vector<const double>,
+                               DiscreteFunctionP0Vector<const TinyVector<1>>,
+                               DiscreteFunctionP0Vector<const TinyVector<2>>,
+                               DiscreteFunctionP0Vector<const TinyVector<3>>,
+                               DiscreteFunctionP0Vector<const TinyMatrix<1>>,
+                               DiscreteFunctionP0Vector<const TinyMatrix<2>>,
+                               DiscreteFunctionP0Vector<const TinyMatrix<3>>>;
 
   Variant m_discrete_function;
 
@@ -70,7 +76,13 @@ class DiscreteFunctionVariant
   DiscreteFunctionVariant(const DiscreteFunctionP0Vector<DataType>& discrete_function)
     : m_discrete_function{DiscreteFunctionP0Vector<const DataType>{discrete_function}}
   {
-    static_assert(std::is_same_v<std::remove_const_t<DataType>, double>,
+    static_assert(std::is_same_v<std::remove_const_t<DataType>, double> or                       //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<1, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<2, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<3, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<1, 1, double>> or   //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<2, 2, double>> or   //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<3, 3, double>>,
                   "DiscreteFunctionP0Vector with this DataType is not allowed in variant");
   }
 
diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp
index d57eae2c45d93968df66fac6ffb3503d3aa74815..6f85610e0f1f87d2bf77ba516ddc9211f0b93fe6 100644
--- a/src/scheme/FluxingAdvectionSolver.cpp
+++ b/src/scheme/FluxingAdvectionSolver.cpp
@@ -86,10 +86,15 @@ class FluxingAdvectionSolver
     m_remapped_list.emplace_back(copy(old_q.cellValues()));
   }
 
+  template <typename DataType>
   void
-  _storeValues(const DiscreteFunctionP0Vector<const double>& old_q)
+  _storeValues(const DiscreteFunctionP0Vector<const DataType>& old_q)
   {
-    m_remapped_list.emplace_back(copy(old_q.cellArrays()));
+    if constexpr (std::is_arithmetic_v<DataType>) {
+      m_remapped_list.emplace_back(copy(old_q.cellArrays()));
+    } else {
+      throw NormalError("remapping DiscreteFunctionP0Vector of non arithmetic data type is not supported");
+    }
   }
 
   template <typename DataType>
@@ -741,8 +746,12 @@ FluxingAdvectionSolver<MeshType>::remap(
           new_variables.push_back(std::make_shared<DiscreteFunctionVariant>(
             DiscreteFunctionT(m_new_mesh, std::get<CellValue<DataType>>(m_remapped_list[i]))));
         } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
-          new_variables.push_back(std::make_shared<DiscreteFunctionVariant>(
-            DiscreteFunctionT(m_new_mesh, std::get<CellArray<DataType>>(m_remapped_list[i]))));
+          if constexpr (std::is_arithmetic_v<DataType>) {
+            new_variables.push_back(std::make_shared<DiscreteFunctionVariant>(
+              DiscreteFunctionT(m_new_mesh, std::get<CellArray<DataType>>(m_remapped_list[i]))));
+          } else {
+            throw NormalError("remapping DiscreteFunctionP0Vector of non arithmetic data type is not supported");
+          }
         } else {
           throw UnexpectedError("invalid discrete function type");
         }
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 2af32e157475b0fe99b9f2cc9cc6c445e05681b3..0e0b2a5afe28330dad7ad0a82b7dc9b277f027c8 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -69,8 +69,28 @@ class PolynomialReconstruction::MutableDiscreteFunctionDPkVariant
                                DiscreteFunctionDPk<3, TinyMatrix<3>>,
 
                                DiscreteFunctionDPkVector<1, double>,
+                               DiscreteFunctionDPkVector<1, TinyVector<1>>,
+                               DiscreteFunctionDPkVector<1, TinyVector<2>>,
+                               DiscreteFunctionDPkVector<1, TinyVector<3>>,
+                               DiscreteFunctionDPkVector<1, TinyMatrix<1>>,
+                               DiscreteFunctionDPkVector<1, TinyMatrix<2>>,
+                               DiscreteFunctionDPkVector<1, TinyMatrix<3>>,
+
                                DiscreteFunctionDPkVector<2, double>,
-                               DiscreteFunctionDPkVector<3, double>>;
+                               DiscreteFunctionDPkVector<2, TinyVector<1>>,
+                               DiscreteFunctionDPkVector<2, TinyVector<2>>,
+                               DiscreteFunctionDPkVector<2, TinyVector<3>>,
+                               DiscreteFunctionDPkVector<2, TinyMatrix<1>>,
+                               DiscreteFunctionDPkVector<2, TinyMatrix<2>>,
+                               DiscreteFunctionDPkVector<2, TinyMatrix<3>>,
+
+                               DiscreteFunctionDPkVector<3, double>,
+                               DiscreteFunctionDPkVector<3, TinyVector<1>>,
+                               DiscreteFunctionDPkVector<3, TinyVector<2>>,
+                               DiscreteFunctionDPkVector<3, TinyVector<3>>,
+                               DiscreteFunctionDPkVector<3, TinyMatrix<1>>,
+                               DiscreteFunctionDPkVector<3, TinyMatrix<2>>,
+                               DiscreteFunctionDPkVector<3, TinyMatrix<3>>>;
 
  private:
   Variant m_mutable_discrete_function_dpk;
@@ -123,7 +143,13 @@ class PolynomialReconstruction::MutableDiscreteFunctionDPkVariant
   MutableDiscreteFunctionDPkVariant(const DiscreteFunctionDPkVector<Dimension, DataType>& discrete_function_dpk)
     : m_mutable_discrete_function_dpk{discrete_function_dpk}
   {
-    static_assert(std::is_same_v<DataType, double>,
+    static_assert(std::is_same_v<DataType, double> or                       //
+                    std::is_same_v<DataType, TinyVector<1, double>> or      //
+                    std::is_same_v<DataType, TinyVector<2, double>> or      //
+                    std::is_same_v<DataType, TinyVector<3, double>> or      //
+                    std::is_same_v<DataType, TinyMatrix<1, 1, double>> or   //
+                    std::is_same_v<DataType, TinyMatrix<2, 2, double>> or   //
+                    std::is_same_v<DataType, TinyMatrix<3, 3, double>>,
                   "DiscreteFunctionDPkVector with this DataType is not allowed in variant");
   }
 
@@ -698,7 +724,7 @@ PolynomialReconstruction::_getNumberOfColumns(
         using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
         if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
           using data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
-          if constexpr (std::is_same_v<data_type, double>) {
+          if constexpr (std::is_arithmetic_v<data_type>) {
             return 1;
           } else if constexpr (is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>) {
             return data_type::Dimension;
@@ -708,7 +734,16 @@ PolynomialReconstruction::_getNumberOfColumns(
             // LCOV_EXCL_STOP
           }
         } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
-          return discrete_function.size();
+          using data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
+          if constexpr (std::is_arithmetic_v<data_type>) {
+            return discrete_function.size();
+          } else if constexpr (is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>) {
+            return discrete_function.size() * data_type::Dimension;
+          } else {
+            // LCOV_EXCL_START
+            throw UnexpectedError("unexpected data type " + demangle<data_type>());
+            // LCOV_EXCL_STOP
+          }
         } else {
           // LCOV_EXCL_START
           throw UnexpectedError("unexpected discrete function type");
@@ -927,7 +962,7 @@ PolynomialReconstruction::_build(
                     } else if constexpr (is_tiny_matrix_v<DataType>) {
                       if constexpr ((DataType::NumberOfColumns == DataType::NumberOfRows) and
                                     (DataType::NumberOfColumns == MeshType::Dimension)) {
-                        throw NotImplementedError("TinyMatrix symmetries for reconstruction");
+                        throw NotImplementedError("TinyMatrix symmetries for reconstruction of DiscreteFunctionP0");
                       }
                       const DataType& qi_qj = discrete_function[cell_i_id] - qj;
                       for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
@@ -945,32 +980,75 @@ PolynomialReconstruction::_build(
                   column_begin += DataType::Dimension;
                 }
               } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
-                using DataType       = std::decay_t<typename DiscreteFunctionT::data_type>;
+                using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
+
                 const auto qj_vector = discrete_function[cell_j_id];
-                size_t index         = 0;
-                for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
-                  const CellId cell_i_id = stencil_cell_list[i];
-                  for (size_t l = 0; l < qj_vector.size(); ++l) {
-                    const DataType& qj         = qj_vector[l];
-                    const DataType& qi_qj      = discrete_function[cell_i_id][l] - qj;
-                    B(index, column_begin + l) = qi_qj;
-                  }
-                }
 
-                for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
-                     ++i_symmetry) {
-                  auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
-                  auto ghost_cell_list = ghost_stencil[cell_j_id];
-                  for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
-                    const CellId cell_i_id = ghost_cell_list[i];
+                if constexpr (std::is_arithmetic_v<DataType>) {
+                  size_t index = 0;
+                  for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
+                    const CellId cell_i_id = stencil_cell_list[i];
                     for (size_t l = 0; l < qj_vector.size(); ++l) {
                       const DataType& qj         = qj_vector[l];
                       const DataType& qi_qj      = discrete_function[cell_i_id][l] - qj;
                       B(index, column_begin + l) = qi_qj;
                     }
                   }
+
+                  for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+                       ++i_symmetry) {
+                    auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+                    auto ghost_cell_list = ghost_stencil[cell_j_id];
+                    for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+                      const CellId cell_i_id = ghost_cell_list[i];
+                      for (size_t l = 0; l < qj_vector.size(); ++l) {
+                        const DataType& qj         = qj_vector[l];
+                        const DataType& qi_qj      = discrete_function[cell_i_id][l] - qj;
+                        B(index, column_begin + l) = qi_qj;
+                      }
+                    }
+                  }
+                } else if constexpr (is_tiny_vector_v<DataType>) {
+                  size_t index = 0;
+                  for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
+                    const CellId cell_i_id = stencil_cell_list[i];
+                    for (size_t l = 0; l < qj_vector.size(); ++l) {
+                      const DataType& qj    = qj_vector[l];
+                      const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
+                      for (size_t kB = column_begin + l * DataType::Dimension, k = 0; k < DataType::Dimension;
+                           ++k, ++kB) {
+                        B(index, kB) = qi_qj[k];
+                      }
+                    }
+                  }
+
+                  for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+                       ++i_symmetry) {
+                    auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+                    auto ghost_cell_list = ghost_stencil[cell_j_id];
+                    if (ghost_cell_list.size() > 0) {
+                      throw NotImplementedError("TinyVector symmetries for reconstruction of DiscreteFunctionP0Vector");
+                    }
+                  }
+                } else if constexpr (is_tiny_matrix_v<DataType>) {
+                  throw NotImplementedError("NIY TinyMatrix data");
+
+                  for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+                       ++i_symmetry) {
+                    auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+                    auto ghost_cell_list = ghost_stencil[cell_j_id];
+                    if (ghost_cell_list.size() > 0) {
+                      throw NotImplementedError("TinyMatrix symmetries for reconstruction of DiscreteFunctionP0Vector");
+                    }
+                  }
+                }
+
+                if constexpr (std::is_arithmetic_v<DataType>) {
+                  column_begin += qj_vector.size();
+                } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
+                  column_begin += qj_vector.size() * DataType::Dimension;
                 }
-                column_begin += qj_vector.size();
+
               } else {
                 // LCOV_EXCL_START
                 throw UnexpectedError("invalid discrete function type");
@@ -1289,6 +1367,24 @@ PolynomialReconstruction::_build(
                           dpk_j_ip1       = X(i, column_begin);
                         }
                         ++column_begin;
+                      } else if constexpr (is_tiny_vector_v<DataType>) {
+                        if (m_descriptor.degree() > 1) {
+                          auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
+                          for (size_t i = 0; i < basis_dimension - 1; ++i) {
+                            auto& dpk_j_0 = dpk_j[0];
+                            for (size_t k = 0; k < DataType::Dimension; ++k) {
+                              dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
+                            }
+                          }
+                        }
+
+                        for (size_t i = 0; i < basis_dimension - 1; ++i) {
+                          auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
+                          for (size_t k = 0; k < DataType::Dimension; ++k) {
+                            dpk_j_ip1[k] = X(i, column_begin + k);
+                          }
+                        }
+                        column_begin += DataType::Dimension;
                       } else {
                         // LCOV_EXCL_START
                         throw UnexpectedError("unexpected data type");
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp
index cdfb532f2cd6a30327f09c5f0b1b8e54e6917a29..3fd31d9527da6954cdbcd74bebc86871c01599cc 100644
--- a/tests/test_PolynomialReconstruction.cpp
+++ b/tests/test_PolynomialReconstruction.cpp
@@ -242,6 +242,82 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
             }
           }
 
+          SECTION("R3 vector data")
+          {
+            using R3 = TinyVector<3>;
+
+            for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+              SECTION(named_mesh.name())
+              {
+                auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+                auto& mesh  = *p_mesh;
+
+                auto vector_affine0 = [](const R1& x) -> R3 {
+                  return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]};
+                };
+
+                auto vector_affine1 = [](const R1& x) -> R3 {
+                  return R3{+1.6 + 0.7 * x[0], -2.1 + 1.2 * x[0], +1.1 - 0.3 * x[0]};
+                };
+
+                auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+                DiscreteFunctionP0Vector<R3> Vh{p_mesh, 2};
+
+                parallel_for(
+                  mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+                    Vh[cell_id][0] = vector_affine0(xj[cell_id]);
+                    Vh[cell_id][1] = vector_affine1(xj[cell_id]);
+                  });
+
+                auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+                auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const R3>>();
+
+                {
+                  double max_mean_error = 0;
+                  for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                    max_mean_error =
+                      std::max(max_mean_error, l2Norm(dpk_Vh(cell_id, 0)(xj[cell_id]) - vector_affine0(xj[cell_id])));
+                    max_mean_error =
+                      std::max(max_mean_error, l2Norm(dpk_Vh(cell_id, 1)(xj[cell_id]) - vector_affine1(xj[cell_id])));
+                  }
+                  REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+                }
+
+                {
+                  double max_slope_error = 0;
+                  {
+                    const TinyVector<3> slope0{+1.7, +2.1, -0.6};
+
+                    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                      for (size_t i = 0; i < R3::Dimension; ++i) {
+                        const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R1{0.1} + xj[cell_id])[i] -
+                                                                        dpk_Vh(cell_id, 0)(xj[cell_id] - R1{0.1})[i]);
+
+                        max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope0[i]));
+                      }
+                    }
+                  }
+
+                  {
+                    const TinyVector<3> slope1{+0.7, +1.2, -0.3};
+
+                    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                      for (size_t i = 0; i < R3::Dimension; ++i) {
+                        const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R1{0.1} + xj[cell_id])[i] -
+                                                                        dpk_Vh(cell_id, 1)(xj[cell_id] - R1{0.1})[i]);
+
+                        max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope1[i]));
+                      }
+                    }
+                  }
+                  REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
+                }
+              }
+            }
+          }
+
           SECTION("list of various types")
           {
             using R3x3 = TinyMatrix<3>;