Add reconstructions for LocalDtHyperelasticSolver

src/scheme/Order2LocalDtHyperelasticSolver.cpp

@@ -450,6 +450,228 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     return std::make_tuple(map1, map2);
   }
 
+     void
+  _vector_limiter(const MeshType& mesh,
+                  const DiscreteFunctionP0<const Rd>& f,
+                  DiscreteFunctionDPk<Dimension, Rd>& DPk_fh) const
+  {
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
+    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor, symmetry_boundary_descriptor_list);
+
+    const auto xj = mesh_data.xj();
+
+    parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+      const Rd fj = f[cell_id];
+
+      Rd f_min = fj;
+      Rd f_max = fj;
+
+      const auto cell_stencil = stencil[cell_id];
+      for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+        for(size_t dim = 0; dim < Dimension; ++dim){
+          f_min[dim] = std::min(f_min[dim], f[cell_stencil[i_cell]][dim]);
+          f_max[dim] = std::max(f_max[dim], f[cell_stencil[i_cell]][dim]);
+        }
+      }
+
+      Rd f_bar_min = fj;
+      Rd f_bar_max = fj;
+
+      for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+        const CellId cell_k_id = cell_stencil[i_cell];
+        const Rd f_xk    = DPk_fh[cell_id](xj[cell_k_id]);
+
+        for(size_t dim = 0; dim < Dimension; ++dim){
+          f_bar_min[dim] = std::min(f_bar_min[dim], f_xk[dim]);
+          f_bar_max[dim] = std::max(f_bar_max[dim], f_xk[dim]);
+        }
+      }
+
+      const double eps = 1E-14;
+      Rd coef1;
+      for(size_t dim = 0; dim < Dimension; ++dim){
+        coef1[dim] = 1;
+        if (std::abs(f_bar_max[dim] - fj[dim]) > eps) {
+          coef1[dim] = (f_max[dim] - fj[dim]) / ((f_bar_max[dim] - fj[dim]));
+        }
+      }
+
+      Rd coef2;
+      for(size_t dim = 0; dim < Dimension; ++dim){
+        coef2[dim] = 1;
+        if (std::abs(f_bar_min[dim] - fj[dim]) > eps) {
+          coef2[dim] = (f_min[dim] - fj[dim]) / ((f_bar_min[dim] - fj[dim]));
+        }
+      }
+
+      double min_coef1 = coef1[0];
+      double min_coef2 = coef2[0];
+      for(size_t dim = 1; dim < Dimension; ++dim){
+        min_coef1 = std::min(min_coef1,coef1[dim]);
+        min_coef2 = std::min(min_coef2,coef2[dim]);
+      }
+
+      const double lambda = std::max(0., std::min(1., std::min(min_coef1, min_coef2)));
+
+      auto coefficients = DPk_fh.coefficients(cell_id);
+
+      coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
+      for (size_t i = 1; i < coefficients.size(); ++i) {
+        coefficients[i] *= lambda;
+      }
+    });
+  }
+
+  CellValue<double> 
+  _scalar_limiter_coefficient(const MeshType& mesh,
+                              const DiscreteFunctionP0<const double>& f,
+                              const DiscreteFunctionDPk<Dimension, const double>& DPk_fh) const
+  {
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
+    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor, symmetry_boundary_descriptor_list);
+
+    const auto xj = mesh_data.xj();
+
+    CellValue<double> lambda{mesh.connectivity()};
+
+    parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+      const double fj = f[cell_id];
+
+      double f_min = fj;
+      double f_max = fj;
+
+      const auto cell_stencil = stencil[cell_id];
+      for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+        f_min = std::min(f_min, f[cell_stencil[i_cell]]);
+        f_max = std::max(f_max, f[cell_stencil[i_cell]]);
+      }
+
+      double f_bar_min = fj;
+      double f_bar_max = fj;
+
+      for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+        const CellId cell_k_id = cell_stencil[i_cell];
+        const double f_xk    = DPk_fh[cell_id](xj[cell_k_id]);
+
+        f_bar_min = std::min(f_bar_min, f_xk);
+        f_bar_max = std::max(f_bar_max, f_xk);
+      }
+
+      const double eps = 1E-14;
+      double coef1     = 1;
+      if (std::abs(f_bar_max - fj) > eps) {
+        coef1 = (f_max - fj) / ((f_bar_max - fj));
+      }
+
+      double coef2 = 1.;
+      if (std::abs(f_bar_min - fj) > eps) {
+        coef2 = (f_min - fj) / ((f_bar_min - fj));
+      }
+
+      lambda[cell_id] = std::max(0., std::min(1., std::min(coef1, coef2)));
+    });
+
+    return lambda;
+  }
+
+  NodeValuePerCell<Rdxd>
+  _compute_tensor_reconstruction(const std::shared_ptr<const MeshVariant>& v_mesh,
+                                 DiscreteTensorFunction sigma,
+                                 PolynomialReconstructionDescriptor reconstruction_descriptor) const
+    {
+        const MeshType& mesh                = *v_mesh->get<MeshType>();
+
+        CellValue<double> limiter_j{mesh.connectivity()};
+        parallel_for(
+          mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
+            limiter_j[j] = 1;
+          }
+        );
+
+        const CellValue<const Rdxd>& sigma_j = copy(sigma.cellValues());
+        std::vector<DiscreteFunctionDPk<Dimension, double>> sigma_coef;
+        std::vector<size_t> row;    
+        std::vector<size_t> col;
+        for(size_t i = 0; i < Dimension; ++i){
+            for(size_t k = i; k < Dimension; ++k){
+
+                CellValue<double> coef{mesh.connectivity()};
+                for(CellId j = 0; j <mesh.numberOfCells(); ++j){
+                    coef[j] = sigma_j[j](i,k);
+                }
+
+                const auto& coef_scalar_function = DiscreteScalarFunction(v_mesh, coef);
+                auto coef_v = std::make_shared<DiscreteFunctionVariant>(coef_scalar_function);
+                auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(coef_v);
+                auto DPk_coef = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const double>>();
+
+                const CellValue<const double>& limiter_ik = _scalar_limiter_coefficient(mesh,coef_scalar_function,DPk_coef);
+                parallel_for(
+                  mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+                    limiter_j[j] = std::min(limiter_j[j],limiter_ik[j]);
+                });
+
+                sigma_coef.push_back(copy(DPk_coef));
+                row.push_back(i);
+                col.push_back(k);
+            }
+        }
+
+        NodeValuePerCell<Rdxd> sigma_lim{mesh.connectivity()} ;
+        const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+        parallel_for(
+          mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
+            const auto& cell_nodes = cell_to_node_matrix[j];
+            for(size_t i = 0; i < Dimension; ++i){
+              for(size_t k = 0; k < Dimension; ++k){
+                for(size_t R = 0; R < cell_nodes.size(); ++R){
+                  sigma_lim(j,R)(i,k) = 0;
+                }
+              }
+            }
+          }
+        );
+
+        const NodeValue<const Rd>& xr = copy(mesh.xr());
+        parallel_for(
+            mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+                const auto& cell_nodes = cell_to_node_matrix[j];
+
+                for(size_t R = 0; R < cell_nodes.size(); ++R){
+                    const NodeId r = cell_nodes[R];
+
+                    for(size_t l = 0; l < sigma_coef.size(); ++l){
+                        const size_t& i = row[l];
+                        const size_t& k = col[l];
+
+                        auto coefficients = sigma_coef[l].coefficients(j);
+
+                        coefficients[0] = (1-limiter_j[j])*sigma_j[j](i,k) + limiter_j[j] * coefficients[0];
+                        for (size_t indice = 1; indice < coefficients.size(); ++indice) {
+                          coefficients[indice] *= limiter_j[j];
+                        }
+
+                        sigma_lim(j,R)(i,k) = sigma_coef[l][j](xr[r]);
+                        if(i != k){
+                          sigma_lim(j,R)(i,k) *= 2;
+                        }
+                    }
+
+                    sigma_lim(j,R) += transpose(sigma_lim(j,R));
+                    sigma_lim(j,R) *= 0.5;
+                }
+
+            }
+        );
+
+        return sigma_lim;
+    }
+
  public:
   std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
   compute_fluxes(const SolverType& solver_type,
@@ -476,10 +698,28 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     const DiscreteScalarFunction& aT    = aT_v->get<DiscreteScalarFunction>();
     const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
 
+    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
+
+    for(auto&& bc_descriptor : bc_descriptor_list){
+      if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
+        symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
+      }
+    }
+
+    PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
+                                                                symmetry_boundary_descriptor_list);
+
+    NodeValuePerCell<Rdxd> sigma_lim = _compute_tensor_reconstruction(i_mesh,sigma,reconstruction_descriptor);
+    auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v);
+    auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+
+    DiscreteFunctionDPk<Dimension, Rd> u_lim  = copy(DPk_uh);
+    this->_vector_limiter(mesh,u,u_lim);
+
     NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
 
     NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
-    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u, sigma);
+    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u_lim, sigma_lim);
 
     const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
     this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
@@ -488,7 +728,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     synchronize(br);
 
     NodeValue<const Rd> ur         = this->_computeUr(mesh, Ar, br);
-    NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);
+    NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, sigma_lim);
 
     return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr));
@@ -547,13 +787,45 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     const DiscreteScalarFunction& aT2    = aT2_v->get<DiscreteScalarFunction>();
     const DiscreteTensorFunction& sigma2 = sigma2_v->get<DiscreteTensorFunction>();
 
+    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list1;
+    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list2;
+
+    for(auto&& bc_descriptor : bc_descriptor_list1){
+      if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
+        symmetry_boundary_descriptor_list1.push_back(bc_descriptor->boundaryDescriptor_shared());
+      }
+    }
+    for(auto&& bc_descriptor : bc_descriptor_list2){
+      if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
+        symmetry_boundary_descriptor_list2.push_back(bc_descriptor->boundaryDescriptor_shared());
+      }
+    }
+
+    PolynomialReconstructionDescriptor reconstruction_descriptor1(IntegrationMethodType::cell_center, 1,
+                                                                 symmetry_boundary_descriptor_list1);
+    PolynomialReconstructionDescriptor reconstruction_descriptor2(IntegrationMethodType::cell_center, 1,
+                                                                 symmetry_boundary_descriptor_list2);
+
+    NodeValuePerCell<Rdxd> sigma1_lim = _compute_tensor_reconstruction(i_mesh1,sigma1,reconstruction_descriptor1);
+    NodeValuePerCell<Rdxd> sigma2_lim = _compute_tensor_reconstruction(i_mesh2,sigma2,reconstruction_descriptor2);
+
+    auto reconstruction1 = PolynomialReconstruction{reconstruction_descriptor1}.build(u1_v);
+    auto DPk_uh1 = reconstruction1[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+    auto reconstruction2 = PolynomialReconstruction{reconstruction_descriptor2}.build(u2_v);
+    auto DPk_uh2 = reconstruction2[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+
+    DiscreteFunctionDPk<Dimension, Rd> u1_lim  = copy(DPk_uh1);
+    this->_vector_limiter(mesh1,u1,u1_lim);
+    DiscreteFunctionDPk<Dimension, Rd> u2_lim  = copy(DPk_uh2);
+    this->_vector_limiter(mesh2,u2,u2_lim);
+
     NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * aL1, rho1 * aT1);
     NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * aL2, rho2 * aT2);
 
     NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
-    NodeValue<Rd> br1   = this->_computeBr(mesh1, Ajr1, u1, sigma1);
+    NodeValue<Rd> br1   = this->_computeBr(mesh1, Ajr1, u1_lim, sigma1_lim);
     NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
-    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2, sigma2);
+    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2_lim, sigma2_lim);
     this->_applyCouplingBC(Ar1, Ar2, br1, br2, map1, map2);
 
     const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
@@ -567,11 +839,11 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     synchronize(br2);
 
     NodeValue<Rd> ur1           = this->_computeUr(mesh1, Ar1, br1);
-    NodeValuePerCell<Rd> Fjr1   = this->_computeFjr(mesh1, Ajr1, ur1, u1, sigma1);
+    NodeValuePerCell<Rd> Fjr1   = this->_computeFjr(mesh1, Ajr1, ur1, u1_lim, sigma1_lim);
     NodeValue<Rd> ur2           = this->_computeUr(mesh2, Ar2, br2);
-    NodeValuePerCell<Rd> Fjr2   = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2);
+    NodeValuePerCell<Rd> Fjr2   = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, sigma2_lim);
     NodeValue<Rd> CR_ur         = this->_computeUr(mesh2, Ar2, br2);
-    NodeValuePerCell<Rd> CR_Fjr = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2);
+    NodeValuePerCell<Rd> CR_Fjr = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, sigma2_lim);
 
     return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
@@ -607,10 +879,28 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     const DiscreteScalarFunction& aT    = aT_v->get<DiscreteScalarFunction>();
     const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
 
+    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
+
+    for(auto&& bc_descriptor : bc_descriptor_list){
+      if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
+        symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
+      }
+    }
+
+    PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
+                                                                symmetry_boundary_descriptor_list);
+
+    NodeValuePerCell<Rdxd> sigma_lim = _compute_tensor_reconstruction(i_mesh,sigma,reconstruction_descriptor);
+    auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v);
+    auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+
+    DiscreteFunctionDPk<Dimension, Rd> u_lim  = copy(DPk_uh);
+    this->_vector_limiter(mesh,u,u_lim);
+
     NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
 
     NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
-    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u, sigma);
+    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u_lim, sigma_lim);
 
     const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
     this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
@@ -619,7 +909,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     synchronize(br);
 
     NodeValue<Rd> ur         = this->_computeUr(mesh, Ar, br);
-    NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);
+    NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, sigma_lim);
 
     this->_applyCouplingBC(mesh, ur, CR_ur, Fjr, CR_Fjr, map2);
 
@@ -680,13 +970,45 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     const DiscreteScalarFunction& aT2    = aT2_v->get<DiscreteScalarFunction>();
     const DiscreteTensorFunction& sigma2 = sigma2_v->get<DiscreteTensorFunction>();
 
+    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list1;
+    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list2;
+
+    for(auto&& bc_descriptor : bc_descriptor_list1){
+      if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
+        symmetry_boundary_descriptor_list1.push_back(bc_descriptor->boundaryDescriptor_shared());
+      }
+    }
+    for(auto&& bc_descriptor : bc_descriptor_list2){
+      if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
+        symmetry_boundary_descriptor_list2.push_back(bc_descriptor->boundaryDescriptor_shared());
+      }
+    }
+
+    PolynomialReconstructionDescriptor reconstruction_descriptor1(IntegrationMethodType::cell_center, 1,
+                                                                 symmetry_boundary_descriptor_list1);
+    PolynomialReconstructionDescriptor reconstruction_descriptor2(IntegrationMethodType::cell_center, 1,
+                                                                 symmetry_boundary_descriptor_list2);
+
+    NodeValuePerCell<Rdxd> sigma1_lim = _compute_tensor_reconstruction(i_mesh1,sigma1,reconstruction_descriptor1);
+    NodeValuePerCell<Rdxd> sigma2_lim = _compute_tensor_reconstruction(i_mesh2,sigma2,reconstruction_descriptor2);
+
+    auto reconstruction1 = PolynomialReconstruction{reconstruction_descriptor1}.build(u1_v);
+    auto DPk_uh1 = reconstruction1[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+    auto reconstruction2 = PolynomialReconstruction{reconstruction_descriptor2}.build(u2_v);
+    auto DPk_uh2 = reconstruction2[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+
+    DiscreteFunctionDPk<Dimension, Rd> u1_lim  = copy(DPk_uh1);
+    this->_vector_limiter(mesh1,u1,u1_lim);
+    DiscreteFunctionDPk<Dimension, Rd> u2_lim  = copy(DPk_uh2);
+    this->_vector_limiter(mesh2,u2,u2_lim);
+
     NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * aL1, rho1 * aT1);
     NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * aL2, rho2 * aT2);
 
     NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
-    NodeValue<Rd> br1   = this->_computeBr(mesh1, Ajr1, u1, sigma1);
+    NodeValue<Rd> br1   = this->_computeBr(mesh1, Ajr1, u1_lim, sigma1_lim);
     NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
-    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2, sigma2);
+    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2_lim, sigma2_lim);
 
     const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
     const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
@@ -701,8 +1023,8 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
     NodeValue<Rd> ur1           = this->_computeUr(mesh1, Ar1, br1);
     NodeValue<Rd> ur2           = this->_computeUr(mesh2, Ar2, br2);
     this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);
-    NodeValuePerCell<Rd> Fjr1   = this->_computeFjr(mesh1, Ajr1, ur1, u1, sigma1);
-    NodeValuePerCell<Rd> Fjr2   = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2);
+    NodeValuePerCell<Rd> Fjr1   = this->_computeFjr(mesh1, Ajr1, ur1, u1_lim, sigma1_lim);
+    NodeValuePerCell<Rd> Fjr2   = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, sigma2_lim);
 
     return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr1),