diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 696d4c02d88b2427521be9ce5f87f1dcc2e9e788..0d8a2dde7bcffcff41e91698fdab7cbfc54fc90e 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -432,7 +432,7 @@ SchemeModule::SchemeModule()
 			                                          std::shared_ptr<const IMesh>,
                                                                   std::shared_ptr<const DiscreteFunctionVariant>,
                                                                   std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+			                                          std::shared_ptr<const DiscreteFunctionVariant>> {
                                 return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),getCommonMesh({rho2, u2, E2, c2, p2})}
                                   .solver()
                                   .apply(LocalDtAcousticSolverHandler::SolverType::Glace, dt1, q, rho1, rho2, u1, u2, E1, E2, c1, c2, p1, p2,
@@ -466,7 +466,7 @@ SchemeModule::SchemeModule()
 			                                          std::shared_ptr<const IMesh>,
                                                                   std::shared_ptr<const DiscreteFunctionVariant>,
                                                                   std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+			                                          std::shared_ptr<const DiscreteFunctionVariant>> {
                                 return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),getCommonMesh({rho2, u2, E2, c2, p2})}
                                   .solver()
                                   .apply(LocalDtAcousticSolverHandler::SolverType::Eucclhyd, dt1, q, rho1, rho2, u1, u2, E1, E2, c1, c2, p1, p2,
diff --git a/src/language/utils/OFStream.cpp b/src/language/utils/OFStream.cpp
index 39b7cc8ecedee3ff76e3d6d54790d25f90d01569..7fbae5890c2e300eb6eb8d244452457149c5fd70 100644
--- a/src/language/utils/OFStream.cpp
+++ b/src/language/utils/OFStream.cpp
@@ -8,6 +8,7 @@ OFStream::OFStream(const std::string& filename)
     m_fstream.open(filename);
     if (m_fstream.is_open()) {
       m_ostream = &m_fstream;
+      m_ostream->precision(15);
     } else {
       std::stringstream error_msg;
       error_msg << "cannot create file " << rang::fgB::yellow << filename << rang::style::reset;
diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp
index b7de42715eff88510119b54507ed32323365c817..77777a6981214be5975ce43ec0f791f09b855bff 100644
--- a/src/mesh/MeshFlatNodeBoundary.cpp
+++ b/src/mesh/MeshFlatNodeBoundary.cpp
@@ -14,7 +14,7 @@ MeshFlatNodeBoundary<Dimension>::_checkBoundaryIsFlat(const TinyVector<Dimension
   const NodeValue<const Rd>& xr = mesh.xr();
 
   bool is_bad = false;
-
+  
   auto node_list = this->m_ref_node_list.list();
   parallel_for(node_list.size(), [=, &is_bad](int r) {
     const Rd& x = xr[node_list[r]];
diff --git a/src/scheme/LocalDtAcousticSolver.cpp b/src/scheme/LocalDtAcousticSolver.cpp
index 0e865976f56f718d89d9f746db16f71e316e6dc7..256cd28506ced0e5a09049a561bfb03d148f3548 100644
--- a/src/scheme/LocalDtAcousticSolver.cpp
+++ b/src/scheme/LocalDtAcousticSolver.cpp
@@ -274,8 +274,8 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
       case IBoundaryConditionDescriptor::Type::coupling: {
 	const CouplingBoundaryConditionDescriptor& coupling_bc_descriptor =
 	  dynamic_cast<const CouplingBoundaryConditionDescriptor&>(*bc_descriptor);
-        bc_list.emplace_back(
-			     CouplingBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),coupling_bc_descriptor.label()));
+	bc_list.emplace_back(
+			     CouplingBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),coupling_bc_descriptor.label()));
         break;
       }
       default: {
@@ -299,20 +299,18 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
                                 const DiscreteScalarFunction& p,
                                 NodeValue<Rdxd>& Ar,
                                 NodeValue<Rd>& br) const;
-  void _applyCouplingBC(const BoundaryConditionList& bc_list1,
-			const BoundaryConditionList& bc_list2,
-			const MeshType& mesh1,
-			const MeshType& mesh2,
-			NodeValue<Rdxd>& Ar1,
+  void _applyCouplingBC(NodeValue<Rdxd>& Ar1,
 			NodeValue<Rdxd>& Ar2,
 			NodeValue<Rd>& br1,
-			NodeValue<Rd>& br2) const;
-  void _applyCouplingBC(const BoundaryConditionList& bc_list,
-			const MeshType& mesh,
+			NodeValue<Rd>& br2,
+			const std::vector<NodeId>& map1,
+			const std::vector<NodeId>& map2) const;
+  void _applyCouplingBC(const MeshType& mesh,
 			NodeValue<Rd>& ur,
 			NodeValue<Rd>& CR_ur,
 			NodeValuePerCell<Rd>& Fjr,
-			NodeValuePerCell<Rd>& CR_Fjr) const;
+			NodeValuePerCell<Rd>& CR_Fjr,
+			const std::vector<NodeId>& map2) const;
 
   void
   _applyBoundaryConditions(const BoundaryConditionList& bc_list,
@@ -361,6 +359,60 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
     return F;
   }
 
+  std::tuple<std::vector<NodeId>,
+	     std::vector<NodeId>>
+  _computeMapping(std::shared_ptr<const IMesh>& i_mesh1,
+		  std::shared_ptr<const IMesh>& i_mesh2,
+		  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+		  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
+  {
+    const MeshType& mesh1 = dynamic_cast<const MeshType&>(*i_mesh1);
+    const MeshType& mesh2 = dynamic_cast<const MeshType&>(*i_mesh2);
+    
+    const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
+    const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
+    
+    std::vector<NodeId> map1;
+    std::vector<NodeId> map2;
+    
+    NodeValue<Rd> xr1 = copy(mesh1.xr());
+    NodeValue<Rd> xr2 = copy(mesh2.xr());
+    for(const auto& boundary_condition1 : bc_list1){
+      std::visit([&](auto&& bc1) {
+	using T1 = std::decay_t<decltype(bc1)>;
+	if constexpr (std::is_same_v<CouplingBoundaryCondition,T1>){
+	  const auto& node_list1 = bc1.nodeList();
+	  for(const auto& boundary_condition2 : bc_list2){
+	    std::visit([&](auto&& bc2) {
+	      using T2 = std::decay_t<decltype(bc2)>;
+	      if constexpr (std::is_same_v<CouplingBoundaryCondition,T2>){
+		if(bc1.label() == bc2.label()){
+		  const auto& node_list2 = bc2.nodeList();
+		  for(size_t i = 0; i<node_list1.size(); i++){
+		    NodeId node_id1 = node_list1[i];
+		    NodeId node_id2;
+		    for(size_t j = 0; j<node_list2.size(); j++){
+		      node_id2 = node_list2[j];
+		      double err = 0;
+		      for(size_t j = 0; j<Dimension; j++){
+			err += (xr1[node_id1][j] - xr2[node_id2][j])*(xr1[node_id1][j] - xr2[node_id2][j]);
+		      }
+		      if(sqrt(err) < 1e-14){
+			map1.push_back(node_id1);
+			map2.push_back(node_id2);
+		      }
+		    }
+		  }
+		}
+	      }
+	    },boundary_condition2);
+	  }
+	}
+      },boundary_condition1);
+    }
+    return std::make_tuple(map1,map2);
+  }
+
  public:
   std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
   compute_fluxes(const SolverType& solver_type,
@@ -420,7 +472,9 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
                  const std::shared_ptr<const DiscreteFunctionVariant>& c2_v,
                  const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
                  const std::shared_ptr<const DiscreteFunctionVariant>& p2_v,
-                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+		 const std::vector<NodeId>& map1,
+		 const std::vector<NodeId>& map2) const
   {
     std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, c1_v, u1_v, p1_v});
     std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, c2_v, u2_v, p2_v});
@@ -440,6 +494,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
       throw NormalError("acoustic solver expects P0 functions");
     }
 
+
     const MeshType& mesh1              = dynamic_cast<const MeshType&>(*i_mesh1);
     const DiscreteScalarFunction& rho1 = rho1_v->get<DiscreteScalarFunction>();
     const DiscreteScalarFunction& c1   = c1_v->get<DiscreteScalarFunction>();
@@ -466,7 +521,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
     this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);
     this->_applyExternalVelocityBC(bc_list1, p1, Ar1, br1);
     this->_applyExternalVelocityBC(bc_list2, p2, Ar2, br2);
-    this->_applyCouplingBC(bc_list1,bc_list2,mesh1,mesh2,Ar1,Ar2,br1,br2);
+    this->_applyCouplingBC(Ar1,Ar2,br1,br2,map1,map2);
 
     synchronize(Ar1);
     synchronize(br1);
@@ -496,7 +551,8 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
                  const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
 		 NodeValue<Rd> CR_ur,
-		 NodeValuePerCell<Rd> CR_Fjr) const 
+		 NodeValuePerCell<Rd> CR_Fjr,
+		 const std::vector<NodeId> map2) const 
   {
     std::shared_ptr i_mesh = getCommonMesh({rho_v, c_v, u_v, p_v});
     if (not i_mesh) {
@@ -517,7 +573,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
 
     NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
     NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u, p);
-
+    
     const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
     this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
     this->_applyExternalVelocityBC(bc_list, p, Ar, br);
@@ -528,7 +584,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
     NodeValue<Rd> ur         = this->_computeUr(mesh, Ar, br);
     NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, p);
 
-    this->_applyCouplingBC(bc_list,mesh,ur,CR_ur,Fjr,CR_Fjr);
+    this->_applyCouplingBC(mesh,ur,CR_ur,Fjr,CR_Fjr,map2);
 
     return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr));
@@ -621,6 +677,62 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
                               Fjr->get<NodeValuePerCell<const Rd>>());
   }
 
+  std::tuple<std::shared_ptr<const IMesh>,
+	     std::shared_ptr<const DiscreteFunctionVariant>,
+	     std::shared_ptr<const DiscreteFunctionVariant>,
+	     std::shared_ptr<const DiscreteFunctionVariant>>
+  mesh_correction(const MeshType& mesh1,
+		  const MeshType& mesh2,
+		  const DiscreteScalarFunction& rho,
+		  const DiscreteVectorFunction& u,
+		  const DiscreteScalarFunction& E,
+		  std::vector<NodeId>& map1,
+		  std::vector<NodeId>& map2) const
+  { 
+    NodeValue<Rd> xr1     = copy(mesh1.xr());
+    NodeValue<Rd> new_xr2 = copy(mesh2.xr());
+
+    size_t n = map1.size();
+
+    for(size_t i = 0; i<n; i++){
+      for(size_t j = 0; j<Dimension; j++){
+	new_xr2[map2[i]][j] = xr1[map1[i]][j];
+      }
+    }
+
+    std::shared_ptr<const MeshType> new_mesh2 = std::make_shared<MeshType>(mesh2.shared_connectivity(), new_xr2);
+
+    CellValue<double> new_rho = copy(rho.cellValues());
+    CellValue<Rd> new_u       = copy(u.cellValues());
+    CellValue<double> new_E   = copy(E.cellValues());
+
+    return {new_mesh2,std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh2, new_rho)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh2, new_u)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh2, new_E))};
+  }
+
+   std::tuple<std::shared_ptr<const IMesh>,
+	     std::shared_ptr<const DiscreteFunctionVariant>,
+	     std::shared_ptr<const DiscreteFunctionVariant>,
+	     std::shared_ptr<const DiscreteFunctionVariant>>
+  mesh_correction(const std::shared_ptr<const IMesh>& i_mesh1,
+		  const std::shared_ptr<const IMesh>& i_mesh2,
+		  const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+		  const std::shared_ptr<const DiscreteFunctionVariant>& u,
+		  const std::shared_ptr<const DiscreteFunctionVariant>& E,
+		  std::vector<NodeId>& map1,
+		  std::vector<NodeId>& map2) const
+  {
+
+    return this->mesh_correction(dynamic_cast<const MeshType&>(*i_mesh1),
+				 dynamic_cast<const MeshType&>(*i_mesh2),
+				 rho->get<DiscreteScalarFunction>(),
+				 u->get<DiscreteVectorFunction>(),
+				 E->get<DiscreteScalarFunction>(),
+				 map1,
+				 map2);
+  }
+
   std::tuple<std::shared_ptr<const IMesh>,
              std::shared_ptr<const DiscreteFunctionVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>,
@@ -646,17 +758,29 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
 	const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
   {
     std::shared_ptr<const IMesh> new_m2 = getCommonMesh({rho2, c2, u2, p2});
+    std::shared_ptr<const IMesh> m1     = getCommonMesh({rho1, c1, u1, p1});
+    
     std::shared_ptr<const DiscreteFunctionVariant> new_rho2 = rho2;
     std::shared_ptr<const DiscreteFunctionVariant> new_u2   = u2;
     std::shared_ptr<const DiscreteFunctionVariant> new_E2   = E2;
 
-    auto  [ur1, Fjr1, ur2, Fjr2,CR_ur,CR_Fjr] = compute_fluxes(solver_type, rho1, c1, u1, p1, bc_descriptor_list1, rho2, c2, u2, p2, bc_descriptor_list2);
-    auto [new_m1,new_rho1,new_u1,new_E1] = apply_fluxes(dt1, rho1, u1, E1, ur1, Fjr1);
+    auto [map1, map2] = _computeMapping(m1,new_m2,bc_descriptor_list1,bc_descriptor_list2);
+    
+    auto [ur1, Fjr1, ur2, Fjr2,CR_ur,CR_Fjr] = compute_fluxes(solver_type, rho1, c1, u1, p1, bc_descriptor_list1, rho2, c2, u2, p2, bc_descriptor_list2,map1,map2);
+    auto [new_m1,new_rho1,new_u1,new_E1]     = apply_fluxes(dt1, rho1, u1, E1, ur1, Fjr1);
     
-    const double& dt2 = dt1/q;
+    double dt2   = dt1/q;
     const double& gamma = 1.4;
+    double sum_dt = 0;
+
+    std::shared_ptr<const DiscreteFunctionVariant> new_p2;
+
+    while(sum_dt != dt1){
+
+      if(sum_dt + dt2 > dt1){
+	dt2 = dt1 - sum_dt;
+      }
 
-    for(size_t i = 0; i<q; i++){
       std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,ur2,Fjr2);
 
       DiscreteScalarFunction rho_d = new_rho2->get<DiscreteScalarFunction>();
@@ -666,12 +790,16 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
       DiscreteScalarFunction p_d   = (gamma-1) * rho_d * eps;
       DiscreteScalarFunction c_d   = sqrt(gamma * p_d / rho_d);
 
-      std::shared_ptr<const DiscreteFunctionVariant> new_p2 = std::make_shared<DiscreteFunctionVariant>(p_d);
+      new_p2 = std::make_shared<DiscreteFunctionVariant>(p_d);
       std::shared_ptr<const DiscreteFunctionVariant> new_c2 = std::make_shared<DiscreteFunctionVariant>(c_d);
 
-      auto [ur2,Fjr2] = compute_fluxes(solver_type,new_rho2,new_c2,new_u2,new_p2,bc_descriptor_list2,CR_ur,CR_Fjr);
+      auto [ur2,Fjr2] = compute_fluxes(solver_type,new_rho2,new_c2,new_u2,new_p2,bc_descriptor_list2,CR_ur,CR_Fjr,map2);
+
+      sum_dt += dt2;
     }
     
+    std::tie(new_m2,new_rho2,new_u2,new_E2) = mesh_correction(new_m1,new_m2,new_rho2,new_u2,new_E2,map1,map2);
+    
     return {new_m1,new_rho1,new_u1,new_E1,new_m2,new_rho2,new_u2,new_E2};
   }
 
@@ -844,128 +972,54 @@ LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyExternalVe
 }
 
 template <size_t Dimension>
-void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCouplingBC(const BoundaryConditionList& bc_list1,
-										      const BoundaryConditionList& bc_list2,
-										      const MeshType& mesh1,
-										      const MeshType& mesh2,
-										      NodeValue<Rdxd>& Ar1,
+void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCouplingBC(NodeValue<Rdxd>& Ar1,
 										      NodeValue<Rdxd>& Ar2,
 										      NodeValue<Rd>& br1,
-										      NodeValue<Rd>& br2) const
+										      NodeValue<Rd>& br2,
+										      const std::vector<NodeId>& map1,
+										      const std::vector<NodeId>& map2) const
 {
-  NodeValue<Rd> xr1 = copy(mesh1.xr());
-  NodeValue<Rd> xr2 = copy(mesh2.xr());
-  for(const auto& boundary_condition1 : bc_list1){
-    std::visit([&](auto&& bc1) {
-      using T1 = std::decay_t<decltype(bc1)>;
-      if constexpr (std::is_same_v<CouplingBoundaryCondition,T1>){
-	const auto& node_list1 = bc1.nodeList();
-	for(const auto& boundary_condition2 : bc_list2){
-	  std::visit([&](auto&& bc2) {
-	    using T2 = std::decay_t<decltype(bc2)>;
-	    if constexpr (std::is_same_v<CouplingBoundaryCondition,T2>){
-	      if(bc1.label() == bc2.label()){
-		const auto& node_list2 = bc2.nodeList();
-		if constexpr (Dimension == 1){
-		  for(size_t i = 0; i<node_list1.size(); i++){
-		    NodeId node_id1 = node_list1[i];
-		    NodeId node_id2;
-		    for(size_t j = 0; j<node_list2.size(); j++){
-		      node_id2 = node_list2[j];
-		      if(abs(xr1[node_id1][0] - xr2[node_id2][0]) < 1e-13){
-			Ar1[node_id1] += Ar2[node_id2];
-			Ar2[node_id2] = Ar1[node_id1];
+  size_t n = map1.size();
 
-			br1[node_id1] += br2[node_id2];
-			br2[node_id2] = br1[node_id1];
+  for(size_t i = 0; i<n; i++){
 
-			break;
-		      }
-		    }
-		  }
-		} else if constexpr (Dimension == 2){
-		  for(size_t i = 0; i<node_list1.size(); i++){
-		    NodeId node_id1 = node_list1[i];
-		    NodeId node_id2;
-		    for(size_t j = 0; j<node_list2.size(); j++){
-		      node_id2 = node_list2[j];
-		      double err = 0;
-		      for(size_t k = 0; k<Dimension; k++){
-			err += (xr1[node_id1][k] - xr2[node_id2][k])*(xr1[node_id1][k] - xr2[node_id2][k]);
-		      }
-		      if(sqrt(err) < 1e-14){
-			Ar1[node_id1] += Ar2[node_id2];
-			Ar2[node_id2] = Ar1[node_id1];
-
-			br1[node_id1] += br2[node_id2];
-			br2[node_id2] = br1[node_id1];
-
-			break;
-		      }
-		    }
-		  }
-		} else {
-		  for(size_t i = 0; i<node_list1.size(); i++){
-		    NodeId node_id1 = node_list1[i];
-		    NodeId node_id2;
-		    for(size_t j = 0; j<node_list2.size(); j++){
-		      node_id2 = node_list2[j];
-		      double err = 0;
-		      for(size_t k = 0; k<Dimension; k++){
-			err += (xr1[node_id1][k] - xr2[node_id2][k])*(xr1[node_id1][k] - xr2[node_id2][k]);
-		      }
-		      if(sqrt(err) < 1e-14){
-			Ar1[node_id1] += Ar2[node_id2];
-			Ar2[node_id2] = Ar1[node_id1];
+    NodeId node_id1 = map1[i];
+    NodeId node_id2 = map2[i];
 
-			br1[node_id1] += br2[node_id2];
-			br2[node_id2] = br1[node_id1];
+    Ar1[node_id1] += Ar2[node_id2];
+    Ar2[node_id2] = Ar1[node_id1];
 
-			break;
-		      }
-		    }
-		  }
-		}
-	      }
-	    }
-	  },boundary_condition2);
-	}
-      }
-    }, boundary_condition1);
+    br1[node_id1] += br2[node_id2];
+    br2[node_id2] = br1[node_id1];
+    
   }
 }
 
 template <size_t Dimension>
-void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCouplingBC(const BoundaryConditionList& bc_list,
- 										      const MeshType& mesh,
+void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCouplingBC(const MeshType& mesh,
  										      NodeValue<Rd>& ur,
  										      NodeValue<Rd>& CR_ur,
  										      NodeValuePerCell<Rd>& Fjr,
-										      NodeValuePerCell<Rd>& CR_Fjr) const
+										      NodeValuePerCell<Rd>& CR_Fjr,
+										      const std::vector<NodeId>& map2) const
 {
-  for(const auto& boundary_condition : bc_list){
-    std::visit([&](auto&& bc) {
-      using T = std::decay_t<decltype(bc)>;
-      if constexpr (std::is_same_v<CouplingBoundaryCondition,T>){
-	
-	const auto& node_list                         = bc.nodeList();
-	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 < node_list.size(); i++){
-	  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);
-	  
-	  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);
-            }
-	  ur[node_id] = CR_ur[node_id];
-	}
-      }
-    }, boundary_condition);
+  const size_t& n = map2.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                       = map2[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);
+
+    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);
+    }
+    ur[node_id] = CR_ur[node_id];
   }
 }