diff --git a/src/language/modules/CouplageModule.cpp b/src/language/modules/CouplageModule.cpp
index 1b818561336b70541a4ac76162bf00d535ade720..fc0553f8820f435bb84ed4543350e51bd521dbdc 100644
--- a/src/language/modules/CouplageModule.cpp
+++ b/src/language/modules/CouplageModule.cpp
@@ -17,6 +17,7 @@
 #include <scheme/FixedBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <utils/CommunicatorManager.hpp>
+#include <utils/CouplingData.hpp>
 #include <utils/Messenger.hpp>
 #include <utils/Serializer.hpp>
 #ifdef PUGS_HAS_MPI
@@ -27,16 +28,42 @@ CouplageModule::CouplageModule()
 {
   this->_addBuiltinFunction("serialize",
                             std::function(
-                              [](std::shared_ptr<const IMesh> i_mesh,
-                                 const std::vector<std::shared_ptr<const IInterfaceDescriptor>>& i_interface_list)
-                                -> void {
-                                const int tag = 20;
-                                return Serializer().apply(i_mesh, i_interface_list, tag);
+                              [](std::shared_ptr<const IMesh> mesh,
+                                 const std::vector<std::shared_ptr<const IInterfaceDescriptor>>& i_interface_list,
+                                 const int64_t dest, const int64_t tag) -> void {
+                                /* const int tag = 20; */
+                                switch (mesh->dimension()) {
+                                case 1: {
+                                  return Serializer<1>().interfaces(mesh, i_interface_list, dest, tag);
+                                }
+                                case 2: {
+                                  return Serializer<2>().interfaces(mesh, i_interface_list, dest, tag);
+                                }
+                                case 3: {
+                                  return Serializer<3>().interfaces(mesh, i_interface_list, dest, tag);
+                                }
+                                default: {
+                                  throw UnexpectedError("invalid mesh dimension");
+                                }
+                                }
                               }));
 
   this->_addBuiltinFunction("serialize", std::function([](std::shared_ptr<const IMesh> mesh, const int64_t dest,
                                                           const int64_t tag) -> void {
-                              return Serializer().apply(mesh, dest, tag);
+                              switch (mesh->dimension()) {
+                              case 1: {
+                                return Serializer<1>().mesh(mesh, dest, tag);
+                              }
+                              case 2: {
+                                return Serializer<2>().mesh(mesh, dest, tag);
+                              }
+                              case 3: {
+                                return Serializer<3>().mesh(mesh, dest, tag);
+                              }
+                              default: {
+                                throw UnexpectedError("invalid mesh dimension");
+                              }
+                              }
                             }));
 
   this->_addBuiltinFunction("serialize", std::function([](const std::shared_ptr<const IBoundaryDescriptor>& boundary,
@@ -44,34 +71,160 @@ CouplageModule::CouplageModule()
                                                           const int64_t dest, const int64_t tag) -> void {
                               std::cout << "\033[01;32m"
                                         << "dest: " << dest << "\ntag: " << tag << "\033[00;00m" << std::endl;
-                              return Serializer().apply(boundary, mesh, location, dest, tag);
+                              switch (mesh->dimension()) {
+                              case 1: {
+                                return Serializer<1>().BC(boundary, mesh, location, dest, tag);
+                              }
+                              case 2: {
+                                return Serializer<2>().BC(boundary, mesh, location, dest, tag);
+                              }
+                              case 3: {
+                                return Serializer<3>().BC(boundary, mesh, location, dest, tag);
+                              }
+                              default: {
+                                throw UnexpectedError("invalid mesh dimension");
+                              }
+                              }
                             }));
 
   this->_addBuiltinFunction("serialize_field",
-                            std::function([](const std::shared_ptr<const DiscreteFunctionVariant>& field) -> void {
-                              const int tag = 300;
-                              return Serializer().apply(field, tag);
+                            std::function([](const std::shared_ptr<const DiscreteFunctionVariant>& field,
+                                             const int64_t dest, const int64_t tag) -> void {
+                              /* const int tag                = 300; */
+                              const std::shared_ptr i_mesh = getCommonMesh({field});
+                              if (not i_mesh) {
+                                throw NormalError("primal discrete functions are not defined on the same mesh");
+                              }
+
+                              switch (i_mesh->dimension()) {
+                              case 1: {
+                                return Serializer<1>().field(field, dest, tag);
+                              }
+                              case 2: {
+                                return Serializer<2>().field(field, dest, tag);
+                              }
+                              case 3: {
+                                return Serializer<3>().field(field, dest, tag);
+                              }
+                              default: {
+                                throw UnexpectedError("invalid mesh dimension");
+                              }
+                              }
                             }));
 
   this->_addBuiltinFunction("serialize_field",
                             std::function([](std::shared_ptr<const IMesh> mesh,
                                              const std::shared_ptr<const IBoundaryDescriptor>& boundary,
-                                             const std::shared_ptr<const ItemValueVariant>& field) -> void {
-                              const int tag = 300;
-                              return Serializer().apply(mesh, boundary, field, tag);
+                                             const std::shared_ptr<const ItemValueVariant>& field, const int64_t dest,
+                                             const int64_t tag) -> void {
+                              /* const int tag = 300; */
+                              /* const std::shared_ptr i_mesh = getCommonMesh({mesh, field}); */
+                              /* if (not i_mesh) { */
+                              /*   throw NormalError("primal discrete functions are not defined on the same mesh"); */
+                              /* } */
+
+                              switch (mesh->dimension()) {
+                              case 1: {
+                                return Serializer<1>().BC_field(mesh, boundary, field, dest, tag);
+                              }
+                              case 2: {
+                                return Serializer<2>().BC_field(mesh, boundary, field, dest, tag);
+                              }
+                              case 3: {
+                                return Serializer<3>().BC_field(mesh, boundary, field, dest, tag);
+                              }
+                              default: {
+                                throw UnexpectedError("invalid mesh dimension");
+                              }
+                              }
                             }));
 
   this->_addBuiltinFunction("cpl_change_mesh_position",
-                            std::function([](std::shared_ptr<const IMesh> mesh) -> std::shared_ptr<const IMesh> {
-                              const int tag = 400;
-                              return Serializer().change_mesh_position(mesh, tag);
+                            std::function([](std::shared_ptr<const IMesh> mesh, const int64_t src,
+                                             const int64_t tag) -> std::shared_ptr<const IMesh> {
+                              /* tag == 400; */
+                              /* src = 1 */
+                              switch (mesh->dimension()) {
+                              case 1: {
+                                return Serializer<1>().change_mesh_position(mesh, src, tag);
+                              }
+                              case 2: {
+                                return Serializer<2>().change_mesh_position(mesh, src, tag);
+                              }
+                              case 3: {
+                                return Serializer<3>().change_mesh_position(mesh, src, tag);
+                              }
+                              default: {
+                                throw UnexpectedError("invalid mesh dimension");
+                              }
+                              }
                             }));
 
+  this->_addBuiltinFunction("get_flux_and_velocity",
+                            std::function(
+
+                              [](const std::shared_ptr<const IMesh>& i_mesh,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const std::shared_ptr<const ItemValueVariant>& ur,
+                                 const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
+                                 const std::shared_ptr<const ItemValueVariant>& ur_saved,
+                                 const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr_saved)
+                                -> std::tuple<std::shared_ptr<const ItemValueVariant>,
+                                              std::shared_ptr<const SubItemValuePerItemVariant>,
+                                              std::shared_ptr<const ItemValueVariant>,
+                                              std::shared_ptr<const SubItemValuePerItemVariant>> {
+                                switch (i_mesh->dimension()) {
+                                case 1: {
+                                  auto myCouplingData = CouplingData<1>::getInstanceOrBuildIt();
+                                  return myCouplingData.get_flux_and_velocity(i_mesh, bc_descriptor_list, ur, Fjr,
+                                                                              ur_saved, Fjr_saved);
+                                }
+                                case 2: {
+                                  auto myCouplingData = CouplingData<2>::getInstanceOrBuildIt();
+                                  return myCouplingData.get_flux_and_velocity(i_mesh, bc_descriptor_list, ur, Fjr,
+                                                                              ur_saved, Fjr_saved);
+                                }
+                                case 3: {
+                                  auto myCouplingData = CouplingData<3>::getInstanceOrBuildIt();
+                                  return myCouplingData.get_flux_and_velocity(i_mesh, bc_descriptor_list, ur, Fjr,
+                                                                              ur_saved, Fjr_saved);
+                                }
+                                default: {
+                                  throw UnexpectedError("invalid mesh dimension");
+                                }
+                                }
+                                /* } */
+                                /* } */
+                              }));
+
+  this->_addBuiltinFunction("update_mesh", std::function(
+
+                                             [](const std::shared_ptr<const IMesh>& i_mesh,
+                                                const std::shared_ptr<const IBoundaryDescriptor>& boundary,
+                                                const int64_t src, const int64_t tag) -> std::shared_ptr<const IMesh> {
+                                               switch (i_mesh->dimension()) {
+                                               case 1: {
+                                                 auto myCouplingData = CouplingData<1>::getInstanceOrBuildIt();
+                                                 return myCouplingData.update_mesh(i_mesh, boundary, src, tag);
+                                               }
+                                               case 2: {
+                                                 auto myCouplingData = CouplingData<2>::getInstanceOrBuildIt();
+                                                 return myCouplingData.update_mesh(i_mesh, boundary, src, tag);
+                                               }
+                                               case 3: {
+                                                 auto myCouplingData = CouplingData<3>::getInstanceOrBuildIt();
+                                                 return myCouplingData.update_mesh(i_mesh, boundary, src, tag);
+                                               }
+                                               default: {
+                                                 throw UnexpectedError("invalid mesh dimension");
+                                               }
+                                               }
+                                             }));
+
   this->_addBuiltinFunction("test", std::function(
 
                                       [](const std::shared_ptr<const IMesh>& i_mesh) -> void {
-                                        /* this->test(); */
-
                                         auto dim = i_mesh->dimension();
                                         std::cout << "in test dim: " << dim << std::endl;
                                         using MeshType = Mesh<Connectivity<2>>;
@@ -105,10 +258,10 @@ void
 CouplageModule::registerOperators() const
 {}
 
-void
-CouplageModule::test()
-{
-  std::cout << "in test" << std::endl;
-}
+/* void */
+/* CouplageModule::_test() */
+/* { */
+/*   std::cout << "in test" << std::endl; */
+/* } */
 
 #endif   // PUGS_HAS_COSTO
diff --git a/src/language/modules/CouplageModule.hpp b/src/language/modules/CouplageModule.hpp
index 810476441303c5f0a0444df243042a01a507ff9b..d074641f831432f0b005344f18de611888c72af3 100644
--- a/src/language/modules/CouplageModule.hpp
+++ b/src/language/modules/CouplageModule.hpp
@@ -18,8 +18,8 @@ class CouplageModule : public BuiltinModule
 
   CouplageModule();
   ~CouplageModule() = default;
-  void test();
   void registerOperators() const final;
+  /* void _test(); */
 };
 
 #endif   // PUGS_HAS_COSTO
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index f2fc4ec5dd6e8593dcad16eea02c405eaa751f8d..1c06883aa9f16f25634ca6dfecb9d79150b8f910 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -47,6 +47,7 @@
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <scheme/VariableBCDescriptor.hpp>
 #include <scheme/VectorDiamondScheme.hpp>
+#include <utils/CouplingData.hpp>
 #include <utils/Socket.hpp>
 
 SchemeModule::SchemeModule()
diff --git a/src/scheme/FixedBoundaryConditionDescriptor.hpp b/src/scheme/FixedBoundaryConditionDescriptor.hpp
index b69648e38364f80ef34552f10f7de41ceebf8f57..995c9483c4f65ab1cfaeefd9f0a7505bbd73d217 100644
--- a/src/scheme/FixedBoundaryConditionDescriptor.hpp
+++ b/src/scheme/FixedBoundaryConditionDescriptor.hpp
@@ -1,6 +1,5 @@
 #ifndef FIXED_BOUNDARY_CONDITION_DESCRIPTOR_HPP
 #define FIXED_BOUNDARY_CONDITION_DESCRIPTOR_HPP
-
 #include <language/utils/FunctionSymbolId.hpp>
 #include <mesh/IBoundaryDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
@@ -43,5 +42,4 @@ class FixedBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
 
   ~FixedBoundaryConditionDescriptor() = default;
 };
-
 #endif   // FIXED_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/HyperelasticSolver.cpp b/src/scheme/HyperelasticSolver.cpp
index dab77ca4cb74d692f3e63bd03cf85b18729271a6..1f073bfca607344ce4c785892ac7cfd7ab3137fa 100644
--- a/src/scheme/HyperelasticSolver.cpp
+++ b/src/scheme/HyperelasticSolver.cpp
@@ -17,6 +17,7 @@
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/CouplingData.hpp>
 #include <utils/Messenger.hpp>
 #include <utils/Socket.hpp>
 
@@ -412,7 +413,8 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
 
     return bc_list;
   }
-
+  /* NodeValue<Rd> CR_ur{mesh.connectivity()}; */
+  /* NodeValuePerCell<Rd> CR_Fjr{mesh.connectivity()}; */
   void _applyPressureBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
   void _applyNormalStressBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
 #ifdef PUGS_HAS_COSTO
@@ -421,6 +423,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
   void _applySymmetryBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
   void _applyVelocityBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
   void _applyCouplingBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
+
   void _applyCouplingBC2(const MeshType& mesh,
                          const BoundaryConditionList& bc_list,
                          NodeValue<Rd>& ur,
@@ -517,7 +520,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
     NodeValue<Rd> ur         = this->_computeUr(mesh, Ar, br);
     NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);
 
-    this->_applyCouplingBC2(mesh, bc_list, ur, Fjr);
+    /* this->_applyCouplingBC2(mesh, bc_list, ur, Fjr); */
 
     return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr));
@@ -935,7 +938,10 @@ HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyCouplingBC(const
               }
             }
           }
+          costo->m = static_cast<int>(recv[shape[0] * shape[1] + shape[2]]);
 
+          std::cout << "\033[01;32m"
+                    << "q recv: " << costo->m << "\033[00;00m" << std::endl;
           parallel_for(
             node_list.size(), PUGS_LAMBDA(size_t i_node) {
               NodeId node_id = node_list[i_node];
@@ -945,7 +951,7 @@ HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyCouplingBC(const
               }
             });
 
-#endif PUGS_HAS_COSTO
+#endif   // PUGS_HAS_COSTO
         }
       },
       boundary_condition);
@@ -973,114 +979,104 @@ HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyCouplingBC2(cons
         /*           << "node_list.size()" << node_list.size() << "\033[00;00m" << std::endl; */
 
 #ifdef PUGS_HAS_COSTO
-          auto costo = parallel::Messenger::getInstance().myCoupling;
-          auto rank  = costo->myGlobalRank();
-          if (rank == 0) {
-            const size_t& n = node_list.size();
 
-            const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
-            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+          auto costo          = parallel::Messenger::getInstance().myCoupling;
+          auto rank           = costo->myGlobalRank();
+          auto myCouplingData = CouplingData<Dimension>::getInstanceOrBuildIt();
 
-            for (size_t i = 0; i < n; i++) {
-              const NodeId node_id                       = node_list[i];
-              const auto& node_to_cell                   = node_to_cell_matrix[node_id];
-              const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
-              /* std::cout << "\033[01;33m" */
-              /*           << "rank_0: " << node_id << " " << node_to_cell.size() << "\033[00;00m" << std::endl; */
-              for (size_t j = 0; j < node_to_cell.size(); j++) {
-                const CellId J       = node_to_cell[j];
-                const unsigned int R = node_local_number_in_its_cells[j];
-                /* Fjr(J, R)            = CR_Fjr(J, R); */
+          if (costo->m == 1) {
+            /* std::cout << "\033[01;33m" */
+            /*           << "(" << rank << ")" */
+            /*           << "NOTHING TO DO " */
+            /*           << "\033[00;00m" << std::endl; */
+          } else {
+            costo->ctm += 1;
+
+            // first call
+            if (costo->ctm == 1) {
+              std::cout << "\033[01;34m"
+                        << "(" << rank << ")"
+                        << "(" << costo->ctm << "|" << costo->m << ")"
+                        << "applied BC2 (SAVE QTT)\n";
+              /* myCouplingData.save_flux_and_velocity(mesh, node_list, ur, Fjr); */
+              myCouplingData.CR_ur  = copy(ur);
+              myCouplingData.CR_Fjr = copy(Fjr);
+              const size_t& n       = node_list.size();
+
+              const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+              const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+              for (size_t i = 0; i < n; i++) {
+                const NodeId node_id                       = node_list[i];
+                const auto& node_to_cell                   = node_to_cell_matrix[node_id];
+                const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+                std::cout << "\033[01;32m"
+                          << "rank_1: " << node_id << " " << node_to_cell.size() << "\033[00;00m" << std::endl;
+
+                for (size_t j = 0; j < node_to_cell.size(); j++) {
+                  const CellId J       = node_to_cell[j];
+                  const unsigned int R = node_local_number_in_its_cells[j];
+                  std::cout << "\033[01;31m" << myCouplingData.CR_Fjr(J, R) << "\n";
+                }
+                std::cout << "\033[01;31m" << myCouplingData.CR_ur[node_id] << "\n";
               }
-              /* ur[node_id] = CR_ur[node_id]; */
-            }
-          }
-          if (rank == 1) {
-            /* std::vector<double> recv; */
-            /* tag = 310; */
-            /* costo->recvData(recv, dest, tag); */
 
-            /* CR_Fjr = ... */
-            /*   CR_ur =  */
-            const size_t& n = node_list.size();
+            } else {
+              std::cout << "\033[01;34m"
+                        << "(" << rank << ")"
+                        << "(" << costo->ctm << "|" << costo->m << ")"
+                        << "applied BC2 (APPLIED QTT FLUX AND UR)\n";
+              /* auto CR_ur      = myCouplingData.getUr(); */
+              /* auto CR_Fjr     = myCouplingData.getFjr(); */
+              const size_t& n = node_list.size();
+
+              const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+              const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+              for (size_t i = 0; i < n; i++) {
+                const NodeId node_id                       = node_list[i];
+                const auto& node_to_cell                   = node_to_cell_matrix[node_id];
+                const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+                std::cout << "\033[01;32m"
+                          << "rank_1: " << node_id << " " << node_to_cell.size() << "\033[00;00m" << std::endl;
+
+                for (size_t j = 0; j < node_to_cell.size(); j++) {
+                  const CellId J       = node_to_cell[j];
+                  const unsigned int R = node_local_number_in_its_cells[j];
+                  std::cout << "\033[01;31m" << myCouplingData.CR_Fjr(J, R) << "\n";
+                }
+                std::cout << "\033[01;31m" << myCouplingData.CR_ur[node_id] << "\n";
+              }
 
-            const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
-            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+              /* const size_t& n = node_list.size(); */
 
-            for (size_t i = 0; i < n; i++) {
-              const NodeId node_id                       = node_list[i];
-              const auto& node_to_cell                   = node_to_cell_matrix[node_id];
-              const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
-              /* std::cout << "\033[01;32m" */
-              /*           << "rank_1: " << node_id << " " << node_to_cell.size() << "\033[00;00m" << std::endl; */
-
-              for (size_t j = 0; j < node_to_cell.size(); j++) {
-                const CellId J       = node_to_cell[j];
-                const unsigned int R = node_local_number_in_its_cells[j];
-                /* Fjr(J, R)            = CR_Fjr(J, R); */
+              /* const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix(); */
+              /* const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); */
+
+              for (size_t i = 0; i < n; i++) {
+                const NodeId node_id                       = node_list[i];
+                const auto& node_to_cell                   = node_to_cell_matrix[node_id];
+                const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+                std::cout << "\033[01;32m"
+                          << "rank_1: " << node_id << " " << node_to_cell.size() << "\033[00;00m" << std::endl;
+
+                for (size_t j = 0; j < node_to_cell.size(); j++) {
+                  const CellId J       = node_to_cell[j];
+                  const unsigned int R = node_local_number_in_its_cells[j];
+                  Fjr(J, R)            = myCouplingData.CR_Fjr(J, R);
+                }
+                ur[node_id] = myCouplingData.CR_ur[node_id];
               }
-              /* ur[node_id] = CR_ur[node_id]; */
+            }
+            if (costo->ctm == costo->m) {
+              costo->ctm = 0;
+              /* std::cout << "\033[01;31m" */
+              /*           << "reset CTM" */
+              /*           << "\033[00;00m" << std::endl; */
             }
           }
 
-        /* const int dest = costo->myGlobalSize() - 1; */
-        /* int tag        = 300; */
-        /* std::vector<int> shape; */
-        /* shape.resize(3); */
-        /* shape[0] = node_list.size(); */
-        /* shape[1] = Dimension + Dimension * Dimension; */
-        /* /\* shape[1] = Dimension; *\/ */
-        /* shape[2] = 0; */
-
-        /* std::vector<double> data; */
-        /* data.resize(shape[0] * shape[1] + shape[2]); */
-        /* Array<TinyVector<Dimension>> br_extracted(node_list.size()); */
-        /* Array<TinyMatrix<Dimension>> Ar_extracted(node_list.size()); */
-        /* parallel_for( */
-        /*   node_list.size(), PUGS_LAMBDA(size_t i_node) { */
-        /*     NodeId node_id = node_list[i_node]; */
-        /*     for (size_t i_dim = 0; i_dim < Dimension; i_dim++) { */
-        /*       br_extracted[i_node] = br[node_id]; */
-        /*       Ar_extracted[i_node] = Ar[node_id]; */
-        /*     } */
-        /*   }); */
-        /* /\*TODO: serializer directement Ar et br pour eviter une copie*\/ */
-        /* for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { */
-        /*   for (unsigned short i_dim = 0; i_dim < Dimension; ++i_dim) { */
-        /*     data[i_node * shape[1] + i_dim] = br_extracted[i_node][i_dim]; */
-        /*     for (size_t j_dim = 0; j_dim < Dimension; j_dim++) { */
-        /*       data[i_node * shape[1] + Dimension + i_dim * Dimension + j_dim] = Ar_extracted[i_node](i_dim,
-j_dim);
-         */
-        /*     } */
-        /*   } */
-        /* } */
-        /* costo->sendData(shape, data, dest, tag); */
-        /* std::vector<double> recv; */
-
-        /* tag = 310; */
-        /* costo->recvData(recv, dest, tag); */
-        /* for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { */
-        /*   for (unsigned short i_dim = 0; i_dim < Dimension; ++i_dim) { */
-        /*     br_extracted[i_node][i_dim] = recv[i_node * shape[1] + i_dim]; */
-        /*     for (size_t j_dim = 0; j_dim < Dimension; j_dim++) { */
-        /*       Ar_extracted[i_node](i_dim, j_dim) = recv[i_node * shape[1] + Dimension + i_dim * Dimension +
-j_dim];
-         */
-        /*     } */
-        /*   } */
-        /* } */
-
-        /* parallel_for( */
-        /*   node_list.size(), PUGS_LAMBDA(size_t i_node) { */
-        /*     NodeId node_id = node_list[i_node]; */
-        /*     for (size_t i_dim = 0; i_dim < Dimension; i_dim++) { */
-        /*       br[node_id] = br_extracted[i_node]; */
-        /*       Ar[node_id] = Ar_extracted[i_node]; */
-        /*     } */
-        /*   }); */
-
-#endif PUGS_HAS_COSTO
+#endif   // PUGS_HAS_COSTO
         }
       },
       boundary_condition);
diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt
index 7f3290a6237a9be9ac3c025652d7fcf95d6cceee..17932cb1e22fdac7dbd415ca8817800a98736e7f 100644
--- a/src/utils/CMakeLists.txt
+++ b/src/utils/CMakeLists.txt
@@ -6,6 +6,7 @@ add_library(
   BacktraceManager.cpp
   CommunicatorManager.cpp
   ConsoleManager.cpp
+  CouplingData.cpp
   Demangle.cpp
   Exceptions.cpp
   FPEManager.cpp
diff --git a/src/utils/CouplingData.cpp b/src/utils/CouplingData.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c44fec1c650acd5694f7f47be9733b2f22ccce0a
--- /dev/null
+++ b/src/utils/CouplingData.cpp
@@ -0,0 +1,6 @@
+#include <utils/CouplingData.hpp>
+#include <utils/Messenger.hpp>
+
+#ifdef PUGS_HAS_COSTO
+
+#endif   // PUGS_HAS_COSTO
diff --git a/src/utils/CouplingData.hpp b/src/utils/CouplingData.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..485fe17b6a8b291def84512cfc8e40d84cb6e7d1
--- /dev/null
+++ b/src/utils/CouplingData.hpp
@@ -0,0 +1,294 @@
+#pragma once
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/ItemValueVariant.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/SubItemValuePerItemVariant.hpp>
+#include <scheme/CouplingBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <utils/Serializer.hpp>
+#include <utils/pugs_config.hpp>
+#include <vector>
+#ifdef PUGS_HAS_COSTO
+
+class IMesh;
+
+template <size_t Dimension>
+class CouplingData
+{
+ private:
+  using MeshType = Mesh<Connectivity<Dimension>>;
+  using Rd       = TinyVector<Dimension>;
+  class CouplingBoundaryCondition;
+
+  using BoundaryCondition1     = std::variant<CouplingBoundaryCondition>;
+  using BoundaryConditionList1 = std::vector<CouplingBoundaryCondition>;
+
+ public:
+  NodeValue<Rd> CR_ur;
+  NodeValuePerCell<Rd> CR_Fjr;
+
+  ~CouplingData() = default;
+
+  static auto&
+  getInstanceOrBuildIt()
+  {
+    static CouplingData<Dimension> cdata;
+    return cdata;
+  }
+  std::tuple<const std::shared_ptr<const ItemValueVariant>,
+             const std::shared_ptr<const SubItemValuePerItemVariant>,
+             const std::shared_ptr<const ItemValueVariant>,
+             const std::shared_ptr<const SubItemValuePerItemVariant>>
+  get_flux_and_velocity(const std::shared_ptr<const IMesh>& i_mesh,
+                        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+                        const std::shared_ptr<const ItemValueVariant>& ur,
+                        const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
+                        const std::shared_ptr<const ItemValueVariant>& ur_saved,
+                        const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr_saved);
+
+  std::shared_ptr<const IMesh> update_mesh(const std::shared_ptr<const IMesh>& i_mesh,
+                                           const std::shared_ptr<const IBoundaryDescriptor>& boundary,
+                                           const int64_t src,
+                                           const int64_t tag);
+
+  BoundaryConditionList1
+  _getBCList(const MeshType& mesh,
+             const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+  {
+    BoundaryConditionList1 coupling_bc_list;
+    for (const auto& bc_descriptor : bc_descriptor_list) {
+      switch (bc_descriptor->type()) {
+      case IBoundaryConditionDescriptor::Type::coupling: {
+        const CouplingBoundaryConditionDescriptor& coupling_bc_descriptor =
+          dynamic_cast<const CouplingBoundaryConditionDescriptor&>(*bc_descriptor);
+        coupling_bc_list.emplace_back(
+          CouplingBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+        break;
+      }
+      default: {
+        break;
+      }
+      }
+    }
+    return coupling_bc_list;
+  }
+
+ private:
+  CouplingData();
+};
+
+template <size_t Dimension>
+class CouplingData<Dimension>::CouplingBoundaryCondition
+{
+ private:
+  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+
+ public:
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_node_boundary.nodeList();
+  }
+
+  CouplingBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary)
+    : m_mesh_node_boundary{mesh_node_boundary}
+  {
+    ;
+  }
+
+  ~CouplingBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+std::shared_ptr<const IMesh>
+CouplingData<Dimension>::update_mesh(const std::shared_ptr<const IMesh>& i_mesh,
+                                     const std::shared_ptr<const IBoundaryDescriptor>& boundary,
+                                     const int64_t src,
+                                     const int64_t tag)
+{
+  auto costo = parallel::Messenger::getInstance().myCoupling;
+  auto rank  = costo->myGlobalRank();
+
+  /* const int src = costo->myGlobalSize() - 1; */
+  std::vector<double> recv;
+
+  costo->recvData(recv, src, tag);
+
+  /* const std::shared_ptr given_mesh               = std::dynamic_pointer_cast<const MeshType>(*i_mesh); */
+  const MeshType& mesh                           = dynamic_cast<const MeshType&>(*i_mesh);
+  MeshNodeBoundary<Dimension> mesh_node_boundary = getMeshNodeBoundary(mesh, *boundary);
+  const auto node_list                           = mesh_node_boundary.nodeList();
+
+  NodeValue<Rd> new_xr = copy(mesh.xr());
+  const size_t& n      = node_list.size();
+
+  const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+
+  for (size_t i = 0; i < n; i++) {
+    const NodeId node_id     = node_list[i];
+    const auto& node_to_cell = node_to_cell_matrix[node_id];
+    /* std::cout << "\033[01;32m" */
+    /*           << "rank_1: " << node_id << " " << node_to_cell.size() << "\033[00;00m" << std::endl; */
+    for (unsigned short j = 0; j < Dimension; ++j) {
+      new_xr[node_id][j] = recv[Dimension * i + j];
+    }
+  }
+  return std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+}
+
+template <size_t Dimension>
+std::tuple<const std::shared_ptr<const ItemValueVariant>,
+           const std::shared_ptr<const SubItemValuePerItemVariant>,
+           const std::shared_ptr<const ItemValueVariant>,
+           const std::shared_ptr<const SubItemValuePerItemVariant>>
+CouplingData<Dimension>::get_flux_and_velocity(
+  const std::shared_ptr<const IMesh>& i_mesh,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const std::shared_ptr<const ItemValueVariant>& ur,
+  const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
+  const std::shared_ptr<const ItemValueVariant>& ur_saved,
+  const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr_saved)
+{
+  auto costo = parallel::Messenger::getInstance().myCoupling;
+  auto rank  = costo->myGlobalRank();
+
+  const MeshType& mesh = dynamic_cast<const MeshType&>(*i_mesh);
+
+  auto coupling_bcs = this->_getBCList(mesh, bc_descriptor_list);
+  auto ur_ok        = copy(ur->get<NodeValue<const Rd>>());
+  auto Fjr_ok       = copy(Fjr->get<NodeValuePerCell<const Rd>>());
+
+  for (const auto& boundary_condition : coupling_bcs) {
+    const auto& node_list = boundary_condition.nodeList();
+    if (costo->m == 1) {
+      std::cout << "\033[01;33m"
+                << "m=1"
+                << "\033[00;00m" << std::endl;
+      return std::make_tuple(std::make_shared<const ItemValueVariant>(ur_ok),
+                             std::make_shared<const SubItemValuePerItemVariant>(Fjr_ok),
+                             std::make_shared<const ItemValueVariant>(ur_ok),
+                             std::make_shared<const SubItemValuePerItemVariant>(Fjr_ok));
+    } else {
+      costo->ctm += 1;
+      if (costo->ctm == 1) {
+        auto ur_saved_ok  = copy(ur->get<NodeValue<const Rd>>());
+        auto Fjr_saved_ok = copy(Fjr->get<NodeValuePerCell<const Rd>>());
+        std::cout << "\033[01;34m"
+                  << "(" << rank << ")"
+                  << "(" << costo->ctm << "|" << costo->m << ")"
+                  << "(SAVE QTT)\n"
+                  << "\033[00;00m" << std::endl;
+
+        if (costo->ctm == costo->m) {
+          costo->ctm = 0;
+        }
+        return std::make_tuple(std::make_shared<const ItemValueVariant>(ur_ok),
+                               std::make_shared<const SubItemValuePerItemVariant>(Fjr_ok),
+                               std::make_shared<const ItemValueVariant>(ur_saved_ok),
+                               std::make_shared<const SubItemValuePerItemVariant>(Fjr_saved_ok));
+      } else {
+        std::cout << "\033[01;34m"
+                  << "(" << rank << ")"
+                  << "(" << costo->ctm << "|" << costo->m << ")"
+                  << "(APPLIED QTT FLUX AND UR)\n"
+                  << "\033[00;00m" << std::endl;
+
+        auto ur_saved_ok  = copy(ur_saved->get<NodeValue<const Rd>>());
+        auto Fjr_saved_ok = copy(Fjr_saved->get<NodeValuePerCell<const Rd>>());
+
+        const size_t& n = node_list.size();
+
+        const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+        const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+        for (size_t i = 0; i < n; i++) {
+          const NodeId node_id                       = node_list[i];
+          const auto& node_to_cell                   = node_to_cell_matrix[node_id];
+          const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+          /* std::cout << "\033[01;32m" */
+          /*           << "rank_1: " << node_id << " " << node_to_cell.size() << "\033[00;00m" << std::endl; */
+
+          for (size_t j = 0; j < node_to_cell.size(); j++) {
+            const CellId J       = node_to_cell[j];
+            const unsigned int R = node_local_number_in_its_cells[j];
+            Fjr_ok(J, R)         = Fjr_saved_ok(J, R);
+          }
+          ur_ok[node_id] = ur_saved_ok[node_id];
+        }
+
+        if (costo->ctm == costo->m) {
+          /* Serializer<Dimension>().BC(std::make_shared<const IBoundaryDescriptor>( */
+          /*                              bc_descriptor_list[0]->boundaryDescriptor()), */
+          /*                            i_mesh, 0, 2, 800); */
+          /* auto ptr = std::make_shared<const IBoundaryDescriptor>(bc_descriptor_list[0]->boundaryDescriptor()); */
+          /* Serializer<Dimension>().BC(ptr, i_mesh, 0, 2, 800); */
+
+          auto costo = parallel::Messenger::getInstance().myCoupling;
+          /* std::vector<int> shape; */
+          /* std::vector<double> pts; */
+          /* const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh); */
+
+          /* /\* const ItemValue<const Rd>& xl = mesh_data.xl(); *\/ */
+          /* const auto xr = mesh->xr(); */
+
+          /* Array<TinyVector<Dimension>> positions(mesh->numberOfNodes()); */
+          /* parallel_for( */
+          /*   mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { */
+          /*     for (unsigned short i = 0; i < Dimension; ++i) { */
+          /*       positions[r][i] = xr[r][i]; */
+          /*     } */
+          /*     /\* for (unsigned short i = Dimension; i < 3; ++i) { *\/ */
+          /*     /\*   positions[r][i] = 0; *\/ */
+          /*     /\* } *\/ */
+          /*   }); */
+
+          /* /\*TODO: serializer directement position pour eviter une copie*\/ */   //
+          /* MeshNodeBoundary<Dimension> mesh_node_boundary = */
+          /*   getMeshNodeBoundary(*mesh, bc_descriptor_list[0]->boundaryDescriptor()); */
+          /* const auto node_list = mesh_node_boundary.nodeList(); */
+          /* pts.resize(node_list.size() * Dimension); */
+
+          /* for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { */
+          /*   const NodeId node_id = node_list[i_node]; */
+          /*   for (unsigned short j = 0; j < Dimension; ++j) { */
+          /*     pts[Dimension * i_node + j] = positions[node_id][j]; */
+          /*   } */
+          /*   /\* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl;
+           * *\/
+           */
+          /* } */
+          /* shape.resize(3); */
+          /* shape[0] = node_list.size(); */
+          /* shape[1] = Dimension; */
+          /* shape[2] = 0; */
+          /* costo->sendData(shape, pts, 2, 800); */
+          /* auto ptr = std::make_shared<const IBoundaryDescriptor>(bc_descriptor_list[0]->boundaryDescriptor()); */
+
+          Serializer<Dimension>().BC(bc_descriptor_list[0]->boundaryDescriptor(), i_mesh, 0, 2, 800);
+          costo->ctm = 0;
+        }
+        return std::make_tuple(std::make_shared<const ItemValueVariant>(ur_ok),
+                               std::make_shared<const SubItemValuePerItemVariant>(Fjr_ok),
+                               std::make_shared<const ItemValueVariant>(ur_saved_ok),
+                               std::make_shared<const SubItemValuePerItemVariant>(Fjr_saved_ok));
+      }
+    }
+  }
+}
+
+template <size_t Dimension>
+CouplingData<Dimension>::CouplingData()
+{
+  std::cout << "\033[01;31m"
+            << "COUPLING DATA CREATION" << Dimension << "\033[00;00m" << std::endl;
+}
+
+using CouplingData3D = CouplingData<3>;
+using CouplingData2D = CouplingData<2>;
+using CouplingData1D = CouplingData<1>;
+
+#endif   // PUGS_HAS_COSTO
diff --git a/src/utils/Serializer.cpp b/src/utils/Serializer.cpp
index dec86016c6573b68a937daed8194b91ba45e708c..1bf1e080ba8a5bfe5c8c31268ae05353e68c4f2e 100644
--- a/src/utils/Serializer.cpp
+++ b/src/utils/Serializer.cpp
@@ -1,635 +1,5 @@
-#include <mesh/MeshNodeBoundary.hpp>
-#include <utils/FPEManager.hpp>
-#include <utils/Messenger.hpp>
 #include <utils/Serializer.hpp>
 
 #ifdef PUGS_HAS_COSTO
 
-Serializer::Serializer() {}
-
-Serializer::~Serializer() = default;
-
-// MESH
-void
-Serializer::apply(std::shared_ptr<const IMesh> i_mesh, const int dest, const int tag)
-{
-  auto costo = parallel::Messenger::getInstance().myCoupling;
-  std::vector<int> shape;
-  std::vector<double> pts;
-  switch (i_mesh->dimension()) {
-  case 1: {
-    using MeshType             = Mesh<Connectivity<1>>;
-    const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    using Rd                      = TinyVector<MeshType::Dimension>;
-    const NodeValue<const Rd>& xr = mesh->xr();
-    Array<TinyVector<1>> positions(mesh->numberOfNodes());
-    parallel_for(
-      mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-        for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
-          positions[r][i] = xr[r][i];
-        }
-      });
-
-    /*TODO: serializer directement position pour eviter une copie*/
-    pts.resize(mesh->numberOfNodes() * MeshType::Dimension);
-    for (unsigned short r = 0; r < mesh->numberOfNodes(); ++r) {
-      for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
-        pts[MeshType::Dimension * r + j] = positions[r][j];
-      }
-    }
-
-    shape.resize(3);
-    shape[0] = mesh->numberOfNodes();
-    shape[1] = i_mesh->dimension();
-    shape[2] = 0;
-    break;
-  }
-  case 2: {
-    using MeshType             = Mesh<Connectivity<2>>;
-    const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    using Rd                      = TinyVector<MeshType::Dimension>;
-    const NodeValue<const Rd>& xr = mesh->xr();
-    Array<TinyVector<2>> positions(mesh->numberOfNodes());
-    parallel_for(
-      mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-        for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
-          positions[r][i] = xr[r][i];
-        }
-      });
-
-    /*TODO: serializer directement position pour eviter une copie*/
-    pts.resize(mesh->numberOfNodes() * MeshType::Dimension);
-    for (unsigned short r = 0; r < mesh->numberOfNodes(); ++r) {
-      for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
-        pts[MeshType::Dimension * r + j] = positions[r][j];
-      }
-    }
-
-    shape.resize(3);
-    shape[0] = mesh->numberOfNodes();
-    shape[1] = i_mesh->dimension();
-    shape[2] = 0;
-    break;
-  }
-  case 3: {
-    using MeshType             = Mesh<Connectivity<3>>;
-    const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    using Rd                      = TinyVector<MeshType::Dimension>;
-    const NodeValue<const Rd>& xr = mesh->xr();
-    Array<TinyVector<2>> positions(mesh->numberOfNodes());
-    parallel_for(
-      mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-        for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
-          positions[r][i] = xr[r][i];
-        }
-      });
-
-    /*TODO: serializer directement position pour eviter une copie*/
-    pts.resize(mesh->numberOfNodes() * MeshType::Dimension);
-    for (unsigned short r = 0; r < mesh->numberOfNodes(); ++r) {
-      for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
-        pts[MeshType::Dimension * r + j] = positions[r][j];
-      }
-    }
-
-    shape.resize(3);
-    shape[0] = mesh->numberOfNodes();
-    shape[1] = i_mesh->dimension();
-    shape[2] = 0;
-    break;
-  }
-  }
-  costo->sendData(shape, pts, dest, tag);
-  /* const int src  = 1; */
-  /* std::vector<int> recv; */
-  /* parallel::Messenger::getInstance().myCoupling->recvData(recv, src, tag + 1); */
-}
-
-// MESH BC (faces)
-void
-Serializer::apply(const std::shared_ptr<const IBoundaryDescriptor>& boundary,
-                  std::shared_ptr<const IMesh> i_mesh,
-                  const int location,
-                  const int dest,
-                  const int tag)
-{
-  auto costo = parallel::Messenger::getInstance().myCoupling;
-  std::vector<int> shape;
-  std::vector<double> pts;
-  std::cout << "\033[01;31m"
-            << "L" << location << " ff " << at_face << "\033[00;00m" << std::endl;
-
-  if (location == at_face) {
-    std::cout << "\033[01;31m"
-              << "ICI AT FACE"
-              << "\033[00;00m" << std::endl;
-    switch (i_mesh->dimension()) {
-    case 1: {
-      using MeshType             = Mesh<Connectivity<1>>;
-      const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-      using MeshDataType = MeshData<1>;
-
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-      /* const ItemValue<const Rd>& xl = mesh_data.xl(); */
-      const auto xl = mesh_data.xl();
-
-      Array<TinyVector<1>> fpositions(mesh->numberOfFaces());
-      parallel_for(
-        mesh->numberOfFaces(), PUGS_LAMBDA(FaceId l) {
-          for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
-            fpositions[l][i] = xl[l][i];
-          }
-          /* for (unsigned short i = MeshType::Dimension; i < 3; ++i) { */
-          /*   positions[r][i] = 0; */
-          /* } */
-        });
-
-      /*TODO: serializer directement fposition pour eviter une copie*/
-      MeshFaceBoundary<1> mesh_face_boundary = getMeshFaceBoundary(*mesh, *boundary);
-      /* mesh_face_boundary.faceList() */
-      const auto face_list = mesh_face_boundary.faceList();
-      pts.resize(face_list.size() * MeshType::Dimension);
-
-      for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-        const FaceId face_id = face_list[i_face];
-        for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
-          pts[MeshType::Dimension * i_face + j] = fpositions[face_id][j];
-        }
-        /* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl; */
-      }
-      shape.resize(3);
-      shape[0] = face_list.size();
-      shape[1] = i_mesh->dimension();
-      shape[2] = 0;
-      break;
-    }
-    case 2: {
-      using MeshType             = Mesh<Connectivity<2>>;
-      const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-      using MeshDataType = MeshData<2>;
-
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-      /* const ItemValue<const Rd>& xl = mesh_data.xl(); */
-      const auto xl = mesh_data.xl();
-
-      Array<TinyVector<2>> fpositions(mesh->numberOfFaces());
-      parallel_for(
-        mesh->numberOfFaces(), PUGS_LAMBDA(FaceId l) {
-          for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
-            fpositions[l][i] = xl[l][i];
-          }
-          /* for (unsigned short i = MeshType::Dimension; i < 3; ++i) { */
-          /*   positions[r][i] = 0; */
-          /* } */
-        });
-
-      /*TODO: serializer directement fposition pour eviter une copie*/
-      MeshFaceBoundary<2> mesh_face_boundary = getMeshFaceBoundary(*mesh, *boundary);
-      /* mesh_face_boundary.faceList() */
-      const auto face_list = mesh_face_boundary.faceList();
-      pts.resize(face_list.size() * MeshType::Dimension);
-
-      for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-        const FaceId face_id = face_list[i_face];
-        for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
-          pts[MeshType::Dimension * i_face + j] = fpositions[face_id][j];
-        }
-        /* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl; */
-      }
-      shape.resize(3);
-      shape[0] = face_list.size();
-      shape[1] = i_mesh->dimension();
-      shape[2] = 0;
-      break;
-    }
-    case 3: {
-      using MeshType             = Mesh<Connectivity<3>>;
-      const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-      using MeshDataType = MeshData<3>;
-
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-      /* const ItemValue<const Rd>& xl = mesh_data.xl(); */
-      const auto xl = mesh_data.xl();
-
-      Array<TinyVector<3>> fpositions(mesh->numberOfFaces());
-      parallel_for(
-        mesh->numberOfFaces(), PUGS_LAMBDA(FaceId l) {
-          for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
-            fpositions[l][i] = xl[l][i];
-          }
-          /* for (unsigned short i = MeshType::Dimension; i < 3; ++i) { */
-          /*   positions[r][i] = 0; */
-          /* } */
-        });
-
-      /*TODO: serializer directement fposition pour eviter une copie*/
-      MeshFaceBoundary<3> mesh_face_boundary = getMeshFaceBoundary(*mesh, *boundary);
-      /* mesh_face_boundary.faceList() */
-      const auto face_list = mesh_face_boundary.faceList();
-      pts.resize(face_list.size() * MeshType::Dimension);
-
-      for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-        const FaceId face_id = face_list[i_face];
-        for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
-          pts[MeshType::Dimension * i_face + j] = fpositions[face_id][j];
-        }
-        /* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl; */
-      }
-      shape.resize(3);
-      shape[0] = face_list.size();
-      shape[1] = i_mesh->dimension();
-      shape[2] = 0;
-
-      break;
-    }
-    }
-    costo->sendData(shape, pts, dest, tag);
-  } else if (location == at_node) {
-    std::cout << "\033[01;31m"
-              << "ICI AT node"
-              << "\033[00;00m" << std::endl;
-    switch (i_mesh->dimension()) {
-    case 1: {
-      using MeshType             = Mesh<Connectivity<1>>;
-      const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-      /* const ItemValue<const Rd>& xl = mesh_data.xl(); */
-      const auto xr = mesh->xr();
-
-      Array<TinyVector<1>> positions(mesh->numberOfNodes());
-      parallel_for(
-        mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-          for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
-            positions[r][i] = xr[r][i];
-          }
-          /* for (unsigned short i = MeshType::Dimension; i < 3; ++i) { */
-          /*   positions[r][i] = 0; */
-          /* } */
-        });
-
-      /*TODO: serializer directement position pour eviter une copie*/
-      MeshNodeBoundary<1> mesh_node_boundary = getMeshNodeBoundary(*mesh, *boundary);
-      /* mesh_face_boundary.faceList() */
-      const auto node_list = mesh_node_boundary.nodeList();
-      pts.resize(node_list.size() * MeshType::Dimension);
-
-      for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
-        const NodeId node_id = node_list[i_node];
-        for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
-          pts[MeshType::Dimension * i_node + j] = positions[node_id][j];
-        }
-        /* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl; */
-      }
-      shape.resize(3);
-      shape[0] = node_list.size();
-      shape[1] = i_mesh->dimension();
-      shape[2] = 0;
-      break;
-    }
-    case 2: {
-      using MeshType             = Mesh<Connectivity<2>>;
-      const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-      /* const ItemValue<const Rd>& xl = mesh_data.xl(); */
-      const auto xr = mesh->xr();
-
-      Array<TinyVector<2>> positions(mesh->numberOfNodes());
-      parallel_for(
-        mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-          for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
-            positions[r][i] = xr[r][i];
-          }
-          /* for (unsigned short i = MeshType::Dimension; i < 3; ++i) { */
-          /*   positions[r][i] = 0; */
-          /* } */
-        });
-
-      /*TODO: serializer directement position pour eviter une copie*/
-      MeshNodeBoundary<2> mesh_node_boundary = getMeshNodeBoundary(*mesh, *boundary);
-      /* mesh_face_boundary.faceList() */
-      const auto node_list = mesh_node_boundary.nodeList();
-      pts.resize(node_list.size() * MeshType::Dimension);
-
-      for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
-        const NodeId node_id = node_list[i_node];
-        for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
-          pts[MeshType::Dimension * i_node + j] = positions[node_id][j];
-        }
-        /* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl; */
-      }
-      shape.resize(3);
-      shape[0] = node_list.size();
-      shape[1] = i_mesh->dimension();
-      shape[2] = 0;
-      break;
-    }
-    case 3: {
-      using MeshType             = Mesh<Connectivity<3>>;
-      const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-      const auto xr = mesh->xr();
-
-      Array<TinyVector<3>> positions(mesh->numberOfNodes());
-      parallel_for(
-        mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-          for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
-            positions[r][i] = xr[r][i];
-          }
-          /* for (unsigned short i = MeshType::Dimension; i < 3; ++i) { */
-          /*   positions[r][i] = 0; */
-          /* } */
-        });
-
-      /*TODO: serializer directement position pour eviter une copie*/
-      MeshNodeBoundary<3> mesh_node_boundary = getMeshNodeBoundary(*mesh, *boundary);
-      /* mesh_face_boundary.faceList() */
-      const auto node_list = mesh_node_boundary.nodeList();
-      pts.resize(node_list.size() * MeshType::Dimension);
-
-      for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
-        const NodeId node_id = node_list[i_node];
-        for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
-          pts[MeshType::Dimension * i_node + j] = positions[node_id][j];
-        }
-        /* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl; */
-      }
-      shape.resize(3);
-      shape[0] = node_list.size();
-      shape[1] = i_mesh->dimension();
-      shape[2] = 0;
-
-      break;
-    }
-    }
-    costo->sendData(shape, pts, dest, tag);
-  } else {
-    throw UnexpectedError(" location must be 0 or 1");
-  }
-
-  /* const int src  = 1; */
-  /* std::vector<int> recv; */
-  /* parallel::Messenger::getInstance().myCoupling->recvData(recv, src, tag + 1); */
-}
-
-void
-Serializer::apply(std::shared_ptr<const IMesh> i_mesh,
-                  const std::vector<std::shared_ptr<const IInterfaceDescriptor>>& i_interface_list,
-                  const int tag)
-{
-  auto costo     = parallel::Messenger::getInstance().myCoupling;
-  const int dest = costo->myGlobalSize() - 1;
-  std::vector<int> shape;
-  std::vector<double> data;
-
-  switch (i_mesh->dimension()) {
-  case 1: {
-    throw NotImplementedError("1D not implemented");
-    break;
-  }
-  case 2: {
-    using MeshType               = Mesh<Connectivity<2>>;
-    const std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    for (const auto& i_interface : i_interface_list) {
-      auto node_interface = getMeshNodeInterface(*p_mesh, *i_interface);
-      auto node_list      = node_interface.nodeList();
-
-      data.resize(node_list.size() * MeshType::Dimension);
-
-      const auto xr = p_mesh->xr();
-
-      Array<TinyVector<2>> positions(p_mesh->numberOfNodes());
-      parallel_for(
-        p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-          for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
-            positions[r][i] = xr[r][i];
-          }
-          /* for (unsigned short i = MeshType::Dimension; i < 3; ++i) { */
-          /*   positions[r][i] = 0; */
-          /* } */
-        });
-
-      /*TODO: serializer directement positions pour eviter une copie */
-      for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
-        const NodeId node_id = node_list[i_node];
-        for (unsigned short i_dim = 0; i_dim < MeshType::Dimension; ++i_dim) {
-          data[i_node * MeshType::Dimension + i_dim] = positions[node_id][i_dim];
-        }
-      }
-
-      shape.resize(3);
-      shape[0] = node_list.size();
-      shape[1] = 2;
-      shape[2] = 0;
-      costo->sendData(shape, data, dest, tag);
-    }
-
-    break;
-  }
-  case 3: {
-    throw NotImplementedError("3D not implemented");
-    break;
-  }
-  }
-}
-// Field at node
-void
-Serializer::apply(std::shared_ptr<const IMesh> i_mesh,
-                  const std::shared_ptr<const IBoundaryDescriptor>& boundary,
-                  const std::shared_ptr<const ItemValueVariant>& field,
-                  const int tag)
-{
-  auto costo     = parallel::Messenger::getInstance().myCoupling;
-  const int dest = costo->myGlobalSize() - 1;
-  std::vector<int> shape;
-  std::vector<double> data;
-
-  /* std::cout << "\033[01;31m" */
-  /*           << "ival" << ival.type() << "\033[00;00m" << std::endl; */
-  /* std::cout << "\033[01;31m" << ival << "\033[00;00m" << std::endl; */
-  /* const auto test = field->mesh(); */
-  /* const std::shared_ptr i_mesh           = getCommonMesh({field, field}); */
-  using MeshType = Mesh<Connectivity<2>>;
-
-  const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-  if (not i_mesh) {
-    throw NormalError("primal discrete functions are not defined on the same mesh");
-  }
-  switch (i_mesh->dimension()) {
-  case 1: {
-    break;
-  }
-  case 2: {
-    const NodeValue<const TinyVector<2>> nField = field->get<NodeValue<const TinyVector<2>>>();
-    using MeshType                              = Mesh<Connectivity<2>>;
-    /* assert(MeshType(i_mesh) == field->mesh()); */
-
-    MeshNodeBoundary<2> mesh_node_boundary = getMeshNodeBoundary(*mesh, *boundary);
-
-    const auto node_list = mesh_node_boundary.nodeList();
-    data.resize(node_list.size() * MeshType::Dimension);
-
-    for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
-      const NodeId node_id = node_list[i_node];
-      for (unsigned short i_dim = 0; i_dim < MeshType::Dimension; ++i_dim) {
-        data[i_node * MeshType::Dimension + i_dim] = nField[node_id][i_dim];
-      }
-    }
-
-    shape.resize(3);
-    shape[0] = node_list.size();
-    shape[1] = 2;
-    shape[2] = 0;
-    costo->sendData(shape, data, dest, tag);
-
-    break;
-  }
-  }
-}
-
-// Field (faces)
-void
-Serializer::apply(const std::shared_ptr<const DiscreteFunctionVariant>& field, const int tag)
-{
-  const std::shared_ptr i_mesh = getCommonMesh({field});
-  if (not i_mesh) {
-    throw NormalError("primal discrete functions are not defined on the same mesh");
-  }
-
-  switch (i_mesh->dimension()) {
-  case 1: {
-    break;
-  }
-  case 2: {
-    using MeshType = Mesh<Connectivity<2>>;
-    /* assert(MeshType(i_mesh) == field->mesh()); */
-    const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-    /* const IBoundaryDescriptor& bd = *boundary; */
-    /* const Mesh<Connectivity<2>> mt = <const MeshType>(i_mesh); */
-
-    /* const auto& n_list = MeshNodeBoundary<2>(mt, bd); */
-
-    using Rd                      = TinyVector<MeshType::Dimension>;
-    const NodeValue<const Rd>& xr = mesh->xr();
-    Array<TinyVector<2>> positions(mesh->numberOfNodes());
-    parallel_for(
-      mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-        for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
-          positions[r][i] = xr[r][i];
-        }
-        /* for (unsigned short i = MeshType::Dimension; i < 3; ++i) { */
-        /*   positions[r][i] = 0; */
-        /* } */
-      });
-
-    std::vector<double> pts;
-    pts.resize(mesh->numberOfNodes() * MeshType::Dimension);
-    for (unsigned short r = 0; r < mesh->numberOfNodes(); ++r) {
-      for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
-        pts[MeshType::Dimension * r + j] = positions[r][j];
-      }
-      /* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl; */
-    }
-
-    auto costo     = parallel::Messenger::getInstance().myCoupling;
-    const int dest = costo->myGlobalSize() - 1;
-    std::vector<int> shape;
-    shape.resize(3);
-    shape[0] = mesh->numberOfNodes();
-    shape[1] = 2;
-    shape[2] = 0;
-    costo->sendData(shape, pts, dest, tag);
-
-    std::vector<int> recv;
-    const int src = costo->myGlobalSize() - 1;
-
-    costo->recvData(recv, src, tag + 1);
-    /* mess.myCoupling->sendData(shape, pts, dest, tag); */
-
-    /* for (size_t i_node = 0; i_node < mesh->numberOfNodes(); i_node++) { */
-    /*   std::cout << "\033[01;31m" << coords [[i_node]] << "\033[00;00m" << std::endl; */
-    /* } */
-    /* const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix(); */
-    /* const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); */
-    /* const auto& coords              = mesh->xr(); */
-    /* for (size_t i_node = 0; i_node < mesh->numberOfCells(); i_node++) { */
-    /*   std::cout << "\033[01;31m" << coords[cell_to_node_matrix[i_node]] << "\033[00;00m" << std::endl; */
-    /* } */
-    break;
-  }
-  }
-}
-
-std::shared_ptr<const IMesh>
-Serializer::change_mesh_position(std::shared_ptr<const IMesh> i_mesh, const int tag)
-{
-  const int src = 1;
-  std::vector<double> recv_xr;
-  parallel::Messenger::getInstance().myCoupling->recvData(recv_xr, src, tag);
-
-  switch (i_mesh->dimension()) {
-  case 1: {
-    using MeshType                   = Mesh<Connectivity<1>>;
-    const std::shared_ptr given_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-    using Rd                         = TinyVector<MeshType::Dimension>;
-    /* const NodeValue<const Rd>& xr    = given_mesh->xr(); */
-
-    NodeValue<Rd> new_xr{given_mesh->connectivity()};
-    parallel_for(
-      given_mesh->numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
-        for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
-          /* new_xr[node_id][j] = xr[node_id][j]; */
-          new_xr[node_id][j] = recv_xr[MeshType::Dimension * node_id + j];
-        }
-      });
-
-    return std::make_shared<MeshType>(given_mesh->shared_connectivity(), new_xr);
-  }
-  case 2: {
-    using MeshType                   = Mesh<Connectivity<2>>;
-    const std::shared_ptr given_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-    using Rd                         = TinyVector<MeshType::Dimension>;
-    /* const NodeValue<const Rd>& xr    = given_mesh->xr(); */
-
-    NodeValue<Rd> new_xr{given_mesh->connectivity()};
-    parallel_for(
-      given_mesh->numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
-        for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
-          /* new_xr[node_id][j] = xr[node_id][j]; */
-          new_xr[node_id][j] = recv_xr[MeshType::Dimension * node_id + j];
-        }
-      });
-
-    return std::make_shared<MeshType>(given_mesh->shared_connectivity(), new_xr);
-  }
-  case 3: {
-    using MeshType                   = Mesh<Connectivity<3>>;
-    const std::shared_ptr given_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-    using Rd                         = TinyVector<MeshType::Dimension>;
-    /* const NodeValue<const Rd>& xr    = given_mesh->xr(); */
-
-    NodeValue<Rd> new_xr{given_mesh->connectivity()};
-    parallel_for(
-      given_mesh->numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
-        for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
-          /* new_xr[node_id][j] = xr[node_id][j]; */
-          new_xr[node_id][j] = recv_xr[MeshType::Dimension * node_id + j];
-        }
-      });
-
-    return std::make_shared<MeshType>(given_mesh->shared_connectivity(), new_xr);
-  }
-  default: {
-    throw NormalError("Must not pass here");
-  }
-  }
-}
-
 #endif   // PUGS_HAS_COSTO
diff --git a/src/utils/Serializer.hpp b/src/utils/Serializer.hpp
index 26779f778c4ceae4b6d1328fa6e8db27064e39d1..7982581d1087be53abdc84abd96a28d2306d71bf 100644
--- a/src/utils/Serializer.hpp
+++ b/src/utils/Serializer.hpp
@@ -11,6 +11,7 @@
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/MeshFaceBoundary.hpp>
 #include <mesh/MeshFlatFaceBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
 #include <mesh/MeshNodeInterface.hpp>
 #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
@@ -18,6 +19,8 @@
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <string>
 #include <type_traits>
+#include <utils/FPEManager.hpp>
+#include <utils/Messenger.hpp>
 #include <utils/PugsAssert.hpp>
 #include <utils/PugsMacros.hpp>
 #include <utils/pugs_config.hpp>
@@ -27,8 +30,14 @@
 #include <mpi.h>
 #endif   // PUGS_HAS_MPI
 
+template <size_t Dimension>
 class Serializer
 {
+ private:
+  using MeshType     = Mesh<Connectivity<Dimension>>;
+  using Rd           = TinyVector<Dimension>;
+  using MeshDataType = MeshData<Dimension>;
+
  public:
   enum locator : int
   {
@@ -36,27 +45,354 @@ class Serializer
     at_face,
     at_cell_center,
   };
-  Serializer();
-  ~Serializer();
-  void apply(std::shared_ptr<const IMesh> i_mesh, const int dest, const int tag);
-  void apply(const std::shared_ptr<const IBoundaryDescriptor>& boundary,
-             std::shared_ptr<const IMesh> i_mesh,
-             const int location,
-             const int dest,
-             const int tag);
-
-  void apply(std::shared_ptr<const IMesh> i_mesh,
-             const std::shared_ptr<const IBoundaryDescriptor>& boundary,
-             const std::shared_ptr<const ItemValueVariant>& field,
-             const int tag);
-
-  void apply(std::shared_ptr<const IMesh> i_mesh,
-             const std::vector<std::shared_ptr<const IInterfaceDescriptor>>& i_interface_list,
-             const int tag);
-
-  void apply(const std::shared_ptr<const DiscreteFunctionVariant>& field, const int tag);
-
-  std::shared_ptr<const IMesh> change_mesh_position(std::shared_ptr<const IMesh> i_mesh, const int tag);
+  Serializer()  = default;
+  ~Serializer() = default;
+  void mesh(std::shared_ptr<const IMesh> i_mesh, const int64_t dest, const int64_t tag);
+  void
+  BC(const std::shared_ptr<const IBoundaryDescriptor>& boundary,
+     std::shared_ptr<const IMesh> i_mesh,
+     const int location,
+     const int64_t dest,
+     const int64_t tag)
+  {
+    this->BC(*boundary, i_mesh, location, dest, tag);
+  }
+
+  void BC(const IBoundaryDescriptor& boundary,
+          std::shared_ptr<const IMesh> i_mesh,
+          const int location,
+          const int64_t dest,
+          const int64_t tag);
+
+  void BC_field(std::shared_ptr<const IMesh> i_mesh,
+                const std::shared_ptr<const IBoundaryDescriptor>& boundary,
+                const std::shared_ptr<const ItemValueVariant>& field,
+                const int64_t dest,
+                const int64_t tag);
+
+  void interfaces(std::shared_ptr<const IMesh> i_mesh,
+                  const std::vector<std::shared_ptr<const IInterfaceDescriptor>>& i_interface_list,
+                  const int64_t dest,
+                  const int64_t tag);
+
+  void field(const std::shared_ptr<const DiscreteFunctionVariant>& field, const int64_t dest, const int64_t tag);
+
+  std::shared_ptr<const IMesh> change_mesh_position(std::shared_ptr<const IMesh> i_mesh,
+                                                    const int64_t src,
+                                                    const int64_t tag);
 };
 
+// MESH node position
+template <size_t Dimension>
+void
+Serializer<Dimension>::mesh(std::shared_ptr<const IMesh> i_mesh, const int64_t dest, const int64_t tag)
+{
+  auto costo = parallel::Messenger::getInstance().myCoupling;
+  std::vector<int> shape;
+  std::vector<double> pts;
+
+  const std::shared_ptr mesh    = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+  const NodeValue<const Rd>& xr = mesh->xr();
+  Array<TinyVector<Dimension>> positions(mesh->numberOfNodes());
+  parallel_for(
+    mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+      for (unsigned short i = 0; i < Dimension; ++i) {
+        positions[r][i] = xr[r][i];
+      }
+    });
+
+  /*TODO: serializer directement position pour eviter une copie*/
+  pts.resize(mesh->numberOfNodes() * Dimension);
+  for (unsigned short r = 0; r < mesh->numberOfNodes(); ++r) {
+    for (unsigned short j = 0; j < Dimension; ++j) {
+      pts[Dimension * r + j] = positions[r][j];
+    }
+  }
+
+  shape.resize(3);
+  shape[0] = mesh->numberOfNodes();
+  shape[1] = Dimension;
+  shape[2] = 0;
+  costo->sendData(shape, pts, dest, tag);
+  /* const int src  = 1; */
+  /* std::vector<int> recv; */
+  /* parallel::Messenger::getInstance().myCoupling->recvData(recv, src, tag + 1); */
+}
+
+// BC nodes position
+template <size_t Dimension>
+void
+Serializer<Dimension>::BC(const IBoundaryDescriptor& boundary,
+                          std::shared_ptr<const IMesh> i_mesh,
+                          const int location,
+                          const int64_t dest,
+                          const int64_t tag)
+{
+  auto costo = parallel::Messenger::getInstance().myCoupling;
+  std::vector<int> shape;
+  std::vector<double> pts;
+  std::cout << "\033[01;31m"
+            << "L" << location << " ff " << at_face << "\033[00;00m" << std::endl;
+
+  if (location == at_face) {
+    std::cout << "\033[01;31m"
+              << "ICI AT FACE"
+              << "\033[00;00m" << std::endl;
+    const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+    /* const ItemValue<const Rd>& xl = mesh_data.xl(); */
+    const auto xl = mesh_data.xl();
+
+    Array<TinyVector<Dimension>> fpositions(mesh->numberOfFaces());
+    parallel_for(
+      mesh->numberOfFaces(), PUGS_LAMBDA(FaceId l) {
+        for (unsigned short i = 0; i < Dimension; ++i) {
+          fpositions[l][i] = xl[l][i];
+        }
+        /* for (unsigned short i = Dimension; i < 3; ++i) { */
+        /*   positions[r][i] = 0; */
+        /* } */
+      });
+
+    /*TODO: serializer directement fposition pour eviter une copie*/
+    MeshFaceBoundary<Dimension> mesh_face_boundary = getMeshFaceBoundary(*mesh, boundary);
+    /* mesh_face_boundary.faceList() */
+    const auto face_list = mesh_face_boundary.faceList();
+    pts.resize(face_list.size() * Dimension);
+
+    for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+      const FaceId face_id = face_list[i_face];
+      for (unsigned short j = 0; j < Dimension; ++j) {
+        pts[Dimension * i_face + j] = fpositions[face_id][j];
+      }
+      /* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl; */
+    }
+    shape.resize(3);
+    shape[0] = face_list.size();
+    shape[1] = Dimension;
+    shape[2] = 0;
+    costo->sendData(shape, pts, dest, tag);
+  } else if (location == at_node) {
+    std::cout << "\033[01;31m"
+              << "ICI AT node"
+              << "\033[00;00m" << std::endl;
+    const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+    /* const ItemValue<const Rd>& xl = mesh_data.xl(); */
+    const auto xr = mesh->xr();
+
+    Array<TinyVector<Dimension>> positions(mesh->numberOfNodes());
+    parallel_for(
+      mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+        for (unsigned short i = 0; i < Dimension; ++i) {
+          positions[r][i] = xr[r][i];
+        }
+        /* for (unsigned short i = Dimension; i < 3; ++i) { */
+        /*   positions[r][i] = 0; */
+        /* } */
+      });
+
+    /*TODO: serializer directement position pour eviter une copie*/
+    MeshNodeBoundary<Dimension> mesh_node_boundary = getMeshNodeBoundary(*mesh, boundary);
+    /* mesh_face_boundary.faceList() */
+    const auto node_list = mesh_node_boundary.nodeList();
+    pts.resize(node_list.size() * Dimension);
+
+    for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+      const NodeId node_id = node_list[i_node];
+      for (unsigned short j = 0; j < Dimension; ++j) {
+        pts[Dimension * i_node + j] = positions[node_id][j];
+      }
+      /* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl; */
+    }
+    shape.resize(3);
+    shape[0] = node_list.size();
+    shape[1] = Dimension;
+    shape[2] = 0;
+    costo->sendData(shape, pts, dest, tag);
+  } else {
+    throw UnexpectedError(" location must be 0 or 1");
+  }
+
+  /* const int src  = 1; */
+  /* std::vector<int> recv; */
+  /* parallel::Messenger::getInstance().myCoupling->recvData(recv, src, tag + 1); */
+}
+
+template <size_t Dimension>
+void
+Serializer<Dimension>::interfaces(std::shared_ptr<const IMesh> i_mesh,
+                                  const std::vector<std::shared_ptr<const IInterfaceDescriptor>>& i_interface_list,
+                                  const int64_t dest,
+                                  const int64_t tag)
+{
+  auto costo = parallel::Messenger::getInstance().myCoupling;
+  /* const int64_t dest = costo->myGlobalSize() - 1; */
+  std::vector<int> shape;
+  std::vector<double> data;
+
+  const std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+  for (const auto& i_interface : i_interface_list) {
+    auto node_interface = getMeshNodeInterface(*p_mesh, *i_interface);
+    auto node_list      = node_interface.nodeList();
+
+    data.resize(node_list.size() * Dimension);
+
+    const auto xr = p_mesh->xr();
+
+    Array<TinyVector<Dimension>> positions(p_mesh->numberOfNodes());
+    parallel_for(
+      p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+        for (unsigned short i = 0; i < Dimension; ++i) {
+          positions[r][i] = xr[r][i];
+        }
+        /* for (unsigned short i = Dimension; i < 3; ++i) { */
+        /*   positions[r][i] = 0; */
+        /* } */
+      });
+
+    /*TODO: serializer directement positions pour eviter une copie */
+    for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+      const NodeId node_id = node_list[i_node];
+      for (unsigned short i_dim = 0; i_dim < Dimension; ++i_dim) {
+        data[i_node * Dimension + i_dim] = positions[node_id][i_dim];
+      }
+    }
+
+    shape.resize(3);
+    shape[0] = node_list.size();
+    shape[1] = Dimension;
+    shape[2] = 0;
+    costo->sendData(shape, data, dest, tag);
+  }
+}
+
+// Field at node
+template <size_t Dimension>
+void
+Serializer<Dimension>::BC_field(std::shared_ptr<const IMesh> i_mesh,
+                                const std::shared_ptr<const IBoundaryDescriptor>& boundary,
+                                const std::shared_ptr<const ItemValueVariant>& field,
+                                const int64_t dest,
+                                const int64_t tag)
+{
+  auto costo = parallel::Messenger::getInstance().myCoupling;
+  /* const int64_t dest = costo->myGlobalSize() - 1; */
+  std::vector<int> shape;
+  std::vector<double> data;
+
+  /* std::cout << "\033[01;31m" */
+  /*           << "ival" << ival.type() << "\033[00;00m" << std::endl; */
+  /* std::cout << "\033[01;31m" << ival << "\033[00;00m" << std::endl; */
+  /* const auto test = field->mesh(); */
+  /* const std::shared_ptr i_mesh           = getCommonMesh({field, field}); */
+
+  const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+  if (not i_mesh) {
+    throw NormalError("primal discrete functions are not defined on the same mesh");
+  }
+  const NodeValue<const TinyVector<Dimension>> nField = field->get<NodeValue<const TinyVector<Dimension>>>();
+  /* assert(MeshType(i_mesh) == field->mesh()); */
+
+  MeshNodeBoundary<Dimension> mesh_node_boundary = getMeshNodeBoundary(*mesh, *boundary);
+
+  const auto node_list = mesh_node_boundary.nodeList();
+  data.resize(node_list.size() * Dimension);
+
+  for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+    const NodeId node_id = node_list[i_node];
+    for (unsigned short i_dim = 0; i_dim < Dimension; ++i_dim) {
+      data[i_node * Dimension + i_dim] = nField[node_id][i_dim];
+    }
+  }
+
+  shape.resize(3);
+  shape[0] = node_list.size();
+  shape[1] = Dimension;
+  shape[2] = 0;
+  costo->sendData(shape, data, dest, tag);
+}
+
+// Field (faces)
+template <size_t Dimension>
+void
+Serializer<Dimension>::field(const std::shared_ptr<const DiscreteFunctionVariant>& field,
+                             const int64_t dest,
+                             const int64_t tag)
+{
+  const std::shared_ptr i_mesh = getCommonMesh({field});
+  if (not i_mesh) {
+    throw NormalError("primal discrete functions are not defined on the same mesh");
+  }
+  /* assert(MeshType(i_mesh) == field->mesh()); */
+  const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+  /* const IBoundaryDescriptor& bd = *boundary; */
+  /* const Mesh<Connectivity<2>> mt = <const MeshType>(i_mesh); */
+
+  /* const auto& n_list = MeshNodeBoundary<2>(mt, bd); */
+
+  const NodeValue<const Rd>& xr = mesh->xr();
+  Array<TinyVector<Dimension>> positions(mesh->numberOfNodes());
+  parallel_for(
+    mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+      for (unsigned short i = 0; i < Dimension; ++i) {
+        positions[r][i] = xr[r][i];
+      }
+      /* for (unsigned short i = Dimension; i < 3; ++i) { */
+      /*   positions[r][i] = 0; */
+      /* } */
+    });
+
+  std::vector<double> pts;
+  pts.resize(mesh->numberOfNodes() * Dimension);
+  for (unsigned short r = 0; r < mesh->numberOfNodes(); ++r) {
+    for (unsigned short j = 0; j < Dimension; ++j) {
+      pts[Dimension * r + j] = positions[r][j];
+    }
+    /* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl; */
+  }
+
+  auto costo = parallel::Messenger::getInstance().myCoupling;
+  /* const int64_t dest = costo->myGlobalSize() - 1; */
+  std::vector<int> shape;
+  shape.resize(3);
+  shape[0] = mesh->numberOfNodes();
+  shape[1] = Dimension;
+  shape[2] = 0;
+  costo->sendData(shape, pts, dest, tag);
+
+  std::vector<int> recv;
+  const int src = costo->myGlobalSize() - 1;
+
+  costo->recvData(recv, src, tag + 1);
+}
+
+template <size_t Dimension>
+std::shared_ptr<const IMesh>
+Serializer<Dimension>::change_mesh_position(std::shared_ptr<const IMesh> i_mesh, const int64_t src, const int64_t tag)
+{
+  std::vector<double> recv_xr;
+  parallel::Messenger::getInstance().myCoupling->recvData(recv_xr, src, tag);
+
+  const std::shared_ptr given_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+  NodeValue<Rd> new_xr{given_mesh->connectivity()};
+  parallel_for(
+    given_mesh->numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+      for (unsigned short j = 0; j < Dimension; ++j) {
+        /* new_xr[node_id][j] = xr[node_id][j]; */
+        new_xr[node_id][j] = recv_xr[Dimension * node_id + j];
+      }
+    });
+
+  return std::make_shared<MeshType>(given_mesh->shared_connectivity(), new_xr);
+}
+
+using Serializer3D = Serializer<3>;
+using Serializer2D = Serializer<2>;
+using Serializer1D = Serializer<1>;
+
 #endif   // PUGS_HAS_COSTO