Add support for variable dirichlet boundary condition

kernel/include/BCs/BoundaryConditions.hpp

@@ -57,7 +57,7 @@ class BoundaryConditions {
  public:
   template <class... Args>
   BoundaryConditions(SpatialDiscretization<T, DIM>* spatial, const Args&... boundaries);
-  void SetBoundaryConditions(mfem::Vector& u);
+  void SetDirichletBoundaryConditions(mfem::Vector& u, Coefficients& coefficients, std::vector<mfem::Vector> auxvars_unk);
   mfem::Array<int> GetEssentialDofs();
   ~BoundaryConditions();
   mfem::Array<int> get_marker_array(const std::string& boundary_type);

kernel/include/BCs/BoundaryConditions.tpp

@@ -163,23 +163,86 @@ mfem::Array<int> BoundaryConditions<T, DIM>::get_marker_array(const std::string&
 }
 
 /**
- * @brief Set boundary conditions
+ * @brief Set dirichlet boundary conditions
  *
  * @param u unknown vector
+ * @param coefficients coefficients list for the variable
+ * @param auxvars_unk unknown vectors of the auxiliary variables
+ * 
  */
 template <class T, int DIM>
-void BoundaryConditions<T, DIM>::SetBoundaryConditions(mfem::Vector& u) {
+void BoundaryConditions<T, DIM>::SetDirichletBoundaryConditions(mfem::Vector& u, Coefficients& coefficients, std::vector<mfem::Vector> auxvars_unk) {
+
+  const int nb_bdr = this->Dirichlet_bdr_.Size();
   mfem::Array<int> tmp_array_bdr(this->Dirichlet_bdr_.Size());
-  for (auto i = 0; i < this->Dirichlet_bdr_.Size(); i++) {
-    tmp_array_bdr = 0;
-    mfem::Array<int> dof;
+  const int coef_size = coefficients.size();
+  Coefficient dirichlet_coef = Coefficient(Glossary::Default, 0.0);
+
+  std::vector<double> u_values;
+  std::vector<double> vaux_values;
+
+  // Inline function to compute Dirichlet coefficient
+  auto compute_dirichlet_coefficient = [&](Coefficient& coef, 
+    const std::span<const double>& values,
+    const std::span<const double>& aux_values) -> double {
+    if (coef.is_scalar()) {
+      return coef.compute();
+    } else {
+      return coef.compute(values, aux_values);
+    }
+  };
+
+  // Loop over boundaries
+  for (int i = 0; i < nb_bdr; i++) {
+
+    // If the boundary is dirichlet
     if (this->Dirichlet_bdr_[i] > 0) {
+
+      // Check if there is a coefficient given and store it
+      bool has_dirichlet_coef = false;
+      for (int l = 0; l < coef_size; l++) {
+        const auto& coef = coefficients[l];
+        if (coef.get_type() == GlossaryType::Dirichlet) {
+          auto bdr_ids = coef.get_bdr_index_coef();
+          if (std::find(bdr_ids.begin(), bdr_ids.end(), i) != bdr_ids.end()) {
+            dirichlet_coef = coef;
+            has_dirichlet_coef = true;
+            break;
+          }
+        }
+      }
+
+      // Get the list of essential true dofs
+      tmp_array_bdr = 0;
       tmp_array_bdr[i] = 1;
-      this->fespace_->GetEssentialTrueDofs(tmp_array_bdr, dof);
-      u.SetSubVector(dof, this->Dirichlet_value_[i]);
+      mfem::Array<int> dof_list;
+      this->fespace_->GetEssentialTrueDofs(tmp_array_bdr, dof_list);
+
+      if (has_dirichlet_coef) {
+        mfem::Vector dirichlet_at_dofs(dof_list.Size());
+        // Loop over essential dofs and compute the coefficient at the dofs
+        for (int j = 0; j < dof_list.Size(); j++) {
+          u_values.clear();
+          vaux_values.clear();
+
+          int dof = dof_list[j];
+          u_values.push_back(u(dof));
+          for (auto aux : auxvars_unk)
+              vaux_values.emplace_back(aux(dof));
+
+          dirichlet_at_dofs[j] = compute_dirichlet_coefficient(dirichlet_coef, std::span<const double>(u_values),
+          std::span<const double>(vaux_values));
+        }
+
+        // SetSubVector with the calculated dirichlet values
+        u.SetSubVector(dof_list, dirichlet_at_dofs);
+      } else {
+        // If no dirichlet coefficient is given, use the constant dirichlet value
+        u.SetSubVector(dof_list, this->Dirichlet_value_[i]);
+      }
     }
   }
-}  // end of SetBoundaryConditions
+}  // end of SetDirichletBoundaryConditions
 
 /**
  * @brief Destroy the Boundary Conditions:: Boundary Conditions object

kernel/include/Glossary/Glossary.hpp

@@ -120,7 +120,8 @@ enum class GlossaryType {
   Neumann,
   RobinA,
   RobinB,
-  ExplicitTime
+  ExplicitTime,
+  Dirichlet
 };
 
 struct GlossaryQuantity {
@@ -359,6 +360,13 @@ static const GlossaryQuantity Robin_a =
 static const GlossaryQuantity Robin_b =
     GlossaryQuantity(GlossaryType::RobinB, GlossaryUnit::None, "B-Robin boundary condition");
 
+/**
+ * @brief Quantity associated with the Dirichlet boundary condition
+ *
+ */
+static const GlossaryQuantity Dirichlet =
+    GlossaryQuantity(GlossaryType::Dirichlet, GlossaryUnit::None, "Dirichlet boundary condition");
+
 /**
  * @brief Quantity associated with the MPI rank
  *

kernel/include/Operators/OperatorBase.tpp

@@ -374,7 +374,15 @@ void OperatorBase<T, DIM>::initialize([[maybe_unused]] const double& initial_tim
 
     this->bcs_.emplace_back(vv.get_boundary_conditions());
     this->ess_tdof_list_.emplace_back(this->bcs_[iv]->GetEssentialDofs());
-    this->bcs_[iv]->SetBoundaryConditions(u);
+
+    // Get the unknown vectors of the auxiliary variables
+    std::vector<mfem::Vector> auxvars_unk;
+    for (const auto& auxvar_vec : this->auxvariables_) {
+      for (const auto& auxvar : auxvar_vec->getVariables()) {
+        auxvars_unk.emplace_back(auxvar.get_unknown());
+      }
+    }
+    this->bcs_[iv]->SetDirichletBoundaryConditions(u, this->coefficients_[iv], auxvars_unk);
     vv.update(u);
     u_vect.emplace_back(u);
   }

kernel/include/Operators/SteadyOperator.tpp

@@ -168,9 +168,16 @@ void SteadyOperator<T, DIM>::solve(std::vector<std::unique_ptr<mfem::Vector>>& v
   this->SetTransientParameters(dt, u_vect);
 
   /// Apply BCs: check if need to be uncomment
+  // Get the unknown vectors of the auxiliary variables
+  // std::vector<mfem::Vector> auxvars_unk;
+  // for (const auto& auxvar_vec : this->auxvariables_) {
+  //   for (const auto& auxvar : auxvar_vec->getVariables()) {
+  //     auxvars_unk.emplace_back(auxvar.get_unknown());
+  //   }
+  // }
   // for (size_t i = 0; i < unk_size; i++) {
   //   auto &unk_i = *(vect_unk[i]);
-  //   this->bcs_[i]->SetBoundaryConditions(unk_i);
+  //   this->bcs_[i]->SetDirichletBoundaryConditions(unk_i, this->coefficients_[i], auxvars_unk);
   // }
 
   // Source term

kernel/include/Operators/TransientOperator.tpp

@@ -534,11 +534,19 @@ void TransientOperator<T, DIM>::ImplicitSolve(const double dt, const mfem::Vecto
   {
     MATools::MATrace::start();
     Catch_Time_Section("ImplicitSolve::ApplyBCs");
+
+    // Get the unknown vectors of the auxiliary variables
+    std::vector<mfem::Vector> auxvars_unk;
+    for (const auto& auxvar_vec : this->auxvariables_) {
+      for (const auto& auxvar : auxvar_vec->getVariables()) {
+        auxvars_unk.emplace_back(auxvar.get_unknown());
+      }
+    }
     auto sc_1 = 0;
     auto sc_2 = sc / fes_size;
     for (int i = 0; i < fes_size; ++i) {
       mfem::Vector v_i(u.GetData() + sc_1, sc_2);
-      this->bcs_[i]->SetBoundaryConditions(v_i);
+      this->bcs_[i]->SetDirichletBoundaryConditions(v_i, this->coefficients_[i], auxvars_unk);
       sc_1 += sc_2;
     }
     reduced_oper->SetParameters(dt, &v);

kernel/include/Problems/Problem.tpp

@@ -394,7 +394,9 @@ void Problem<OPE, VAR, PST>::do_time_step(
     [[maybe_unused]] const std::vector<std::vector<std::string>>& unks_info) {
   const size_t unk_size = vect_unk.size();
 
-  this->set_time_coefficients(next_time);
+  // Set time for coefficients, /!\ this is correct only for implicit time scheme
+  // TODO: use effective time step according to the time scheme
+  this->set_time_coefficients(next_time+current_time_step);
 
   this->oper_.setGeometry(this->geometry_);
   this->oper_.solve(vect_unk, next_time, current_time, current_time_step, iter);

tests/HeatTransfer/2D/CMakeLists.txt

@@ -2,3 +2,4 @@ add_subdirectory(test1)
 add_subdirectory(test2)
 add_subdirectory(test3)
 add_subdirectory(test4)
+add_subdirectory(test5)

tests/HeatTransfer/2D/test5/CMakeLists.txt

@@ -0,0 +1,6 @@
+execute_process(COMMAND python3 ${BUILD_SCRIPT_DIR}/GenerateCoefficient.py -r  -f  ${CMAKE_CURRENT_SOURCE_DIR}/coefficient.json
+    WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
+
+create_test("HeatTransfer2Dtest5" "HeatTransfer2Dtest5" FALSE "2D EXE HEAT" 1)
+configure_file(${CMAKE_SOURCE_DIR}/tests/tools/convergence_study.py ${CMAKE_CURRENT_BINARY_DIR}/convergence_study.py COPYONLY)
+create_col_comparison_convergence("CompareHeat2Dtest5Convergence" "convergence_output_ref.csv" "convergence_output.csv" -1 absolute 1e-16 FALSE "HeatTransfer2Dtest5" "2D Heat" 1. 100000)

tests/HeatTransfer/2D/test5/Coefficient.hpp

@@ -0,0 +1,190 @@
+/**
+ *
+ * Copyright CEA (C) 2026
+ *
+ * This file is part of SLOTH.
+ *
+ * SLOTH is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * SLOTH is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#include <algorithm>
+#include <cmath>
+#include <functional>
+#include <numeric>
+#include <span>
+#include <vector>
+
+
+#include "Options/PhysicalPropertiesOptions.hpp"
+
+#include "Coefficients/FunctionCoefficient.hpp"
+
+
+#pragma once
+
+/**
+ *
+ * @brief C++ function of the analytical expression
+ *
+ *       F = -pi*t*cos(pi*x)
+ */
+class NeumannCoefficient : public FunctionCoefficient {
+ private:
+  double prefactor_;
+ protected:
+  std::function<double(const std::span<const double>&,const std::span<const double>&, const std::span<const double>&, const unsigned int dimension)> F() final;
+  std::function<std::vector<double>(const std::span<const double>&,const std::span<const double>&, const std::span<const double>&, const unsigned int dimension)> GradientF() final;
+  std::function<std::vector<double>(const std::span<const double>&,const std::span<const double>&,const std::span<const double>&, const unsigned int dimension)> HessianF() final;
+
+ public:
+  NeumannCoefficient() : prefactor_(1.0) {}
+ explicit NeumannCoefficient(const double prefactor): prefactor_(prefactor) {}
+  virtual ~NeumannCoefficient() = default;
+};
+
+/**
+ *
+ * @brief C++ function of the expression
+ *
+ *
+ * @return std::function<double(const std::span<const double>&,const std::span<const double>&, const std::span<const double>&, const unsigned int dimension)>
+ */
+ std::function<double(const std::span<const double>&,const std::span<const double>&, const std::span<const double>&, const unsigned int dimension)> NeumannCoefficient::F() {
+  auto func = [&](const std::span<const double>& input_vector, [[maybe_unused]] const std::span<const double>&,const std::span<const double>& auxiliary_vector, [[maybe_unused]] const unsigned int dimension) {
+    double T = input_vector[0];
+    double x = auxiliary_vector[0];
+    double y = auxiliary_vector[1];
+    double t = this->time_;
+ double   F = -M_PI*t*std::cos(M_PI*x);
+    return this->prefactor_ * F;
+  };
+  return func;
+}
+
+/**
+ *
+ * @brief Gradient
+ * 
+ * @return std::function<std::vector<double>(const std::span<const double>&,const std::span<const double>&, const unsigned int dimension)> 
+ */
+std::function<std::vector<double>(const std::span<const double>&,const std::span<const double>&, const std::span<const double>&, const unsigned int dimension)> NeumannCoefficient::GradientF() {
+  auto func = [&](const std::span<const double>& input_vector, [[maybe_unused]] const std::span<const double>&,const std::span<const double>& auxiliary_vector, [[maybe_unused]] const unsigned int dimension) {
+    double T = input_vector[0];
+    double x = auxiliary_vector[0];
+    double y = auxiliary_vector[1];
+    double t = this->time_;
+    std::vector<double> gradient(1);
+    gradient[0] = this->prefactor_ * (0);
+    return gradient;
+  };
+  return func;
+}
+
+/**
+ *
+ * @brief Hessian
+ * @remark Hessian matrix stored in vector : H(i,j)->H(i*n+j)
+ * 
+ * @return std::function<std::vector<double>(const std::span<const double>&,const std::span<const double>&, const std::span<const double>&, const unsigned int dimension)> 
+ */
+std::function<std::vector<double>(const std::span<const double>&,const std::span<const double>&, const std::span<const double>&, const unsigned int dimension)> NeumannCoefficient::HessianF() {
+  auto func = [&](const std::span<const double>& input_vector, [[maybe_unused]] const std::span<const double>&,const std::span<const double>& auxiliary_vector, [[maybe_unused]] const unsigned int dimension) {
+    double T = input_vector[0];
+    double x = auxiliary_vector[0];
+    double y = auxiliary_vector[1];
+    double t = this->time_;
+    std::vector<double> hessian(1);
+    hessian[0] = this->prefactor_ * (0);
+    return hessian;
+  };
+  return func;
+}
+/**
+ *
+ * @brief C++ function of the analytical expression
+ *
+ *       F = t*sin(pi*y)
+ */
+class DirichletCoefficient : public FunctionCoefficient {
+ private:
+  double prefactor_;
+ protected:
+  std::function<double(const std::span<const double>&,const std::span<const double>&, const std::span<const double>&, const unsigned int dimension)> F() final;
+  std::function<std::vector<double>(const std::span<const double>&,const std::span<const double>&, const std::span<const double>&, const unsigned int dimension)> GradientF() final;
+  std::function<std::vector<double>(const std::span<const double>&,const std::span<const double>&,const std::span<const double>&, const unsigned int dimension)> HessianF() final;
+
+ public:
+  DirichletCoefficient() : prefactor_(1.0) {}
+ explicit DirichletCoefficient(const double prefactor): prefactor_(prefactor) {}
+  virtual ~DirichletCoefficient() = default;
+};
+
+/**
+ *
+ * @brief C++ function of the expression
+ *
+ *
+ * @return std::function<double(const std::span<const double>&,const std::span<const double>&, const std::span<const double>&, const unsigned int dimension)>
+ */
+ std::function<double(const std::span<const double>&,const std::span<const double>&, const std::span<const double>&, const unsigned int dimension)> DirichletCoefficient::F() {
+  auto func = [&](const std::span<const double>& input_vector, [[maybe_unused]] const std::span<const double>&,const std::span<const double>& auxiliary_vector, [[maybe_unused]] const unsigned int dimension) {
+    double T = input_vector[0];
+    double x = auxiliary_vector[0];
+    double y = auxiliary_vector[1];
+    double t = this->time_;
+ double   F = t*std::sin(M_PI*y);
+    return this->prefactor_ * F;
+  };
+  return func;
+}
+
+/**
+ *
+ * @brief Gradient
+ * 
+ * @return std::function<std::vector<double>(const std::span<const double>&,const std::span<const double>&, const unsigned int dimension)> 
+ */
+std::function<std::vector<double>(const std::span<const double>&,const std::span<const double>&, const std::span<const double>&, const unsigned int dimension)> DirichletCoefficient::GradientF() {
+  auto func = [&](const std::span<const double>& input_vector, [[maybe_unused]] const std::span<const double>&,const std::span<const double>& auxiliary_vector, [[maybe_unused]] const unsigned int dimension) {
+    double T = input_vector[0];
+    double x = auxiliary_vector[0];
+    double y = auxiliary_vector[1];
+    double t = this->time_;
+    std::vector<double> gradient(1);
+    gradient[0] = this->prefactor_ * (0);
+    return gradient;
+  };
+  return func;
+}
+
+/**
+ *
+ * @brief Hessian
+ * @remark Hessian matrix stored in vector : H(i,j)->H(i*n+j)
+ * 
+ * @return std::function<std::vector<double>(const std::span<const double>&,const std::span<const double>&, const std::span<const double>&, const unsigned int dimension)> 
+ */
+std::function<std::vector<double>(const std::span<const double>&,const std::span<const double>&, const std::span<const double>&, const unsigned int dimension)> DirichletCoefficient::HessianF() {
+  auto func = [&](const std::span<const double>& input_vector, [[maybe_unused]] const std::span<const double>&,const std::span<const double>& auxiliary_vector, [[maybe_unused]] const unsigned int dimension) {
+    double T = input_vector[0];
+    double x = auxiliary_vector[0];
+    double y = auxiliary_vector[1];
+    double t = this->time_;
+    std::vector<double> hessian(1);
+    hessian[0] = this->prefactor_ * (0);
+    return hessian;
+  };
+  return func;
+}